72 NumGlobalNonzeros_(0),
83 PetscObjectGetComm( (PetscObject)Amat, &comm);
91 MatGetType(Amat, &MatType_);
92 if ( strcmp(MatType_,MATSEQAIJ) != 0 && strcmp(MatType_,MATMPIAIJ) != 0 ) {
93 sprintf(errMsg,
"PETSc matrix must be either seqaij or mpiaij (but it is %s)",MatType_);
94 throw Comm_->ReportError(errMsg,-1);
98 if (strcmp(MatType_,MATMPIAIJ) == 0) {
100 aij = (Mat_MPIAIJ*)Amat->data;
102 else if (strcmp(MatType_,MATSEQAIJ) == 0) {
105 int numLocalRows, numLocalCols;
106 ierr = MatGetLocalSize(Amat,&numLocalRows,&numLocalCols);
108 sprintf(errMsg,
"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetLocalSize() returned error code %d",__LINE__,ierr);
109 throw Comm_->ReportError(errMsg,-1);
111 NumMyRows_ = numLocalRows;
112 NumMyCols_ = numLocalCols;
114 if (mt == PETSC_MPI_AIJ)
115 NumMyCols_ += aij->B->cmap->n;
117 ierr = MatGetInfo(Amat,MAT_LOCAL,&info);
119 sprintf(errMsg,
"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetInfo() returned error code %d",__LINE__,ierr);
120 throw Comm_->ReportError(errMsg,-1);
122 NumMyNonzeros_ = (int) info.nz_used;
123 Comm_->SumAll(&(info.nz_used), &NumGlobalNonzeros_, 1);
127 int rowStart, rowEnd;
128 ierr = MatGetOwnershipRange(Amat,&rowStart,&rowEnd);
130 sprintf(errMsg,
"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetOwnershipRange() returned error code %d",__LINE__,ierr);
131 throw Comm_->ReportError(errMsg,-1);
134 PetscRowStart_ = rowStart;
135 PetscRowEnd_ = rowEnd;
137 int* MyGlobalElements =
new int[rowEnd-rowStart];
138 for (
int i=0; i<rowEnd-rowStart; i++)
139 MyGlobalElements[i] = rowStart+i;
141 ierr = MatGetInfo(Amat,MAT_GLOBAL_SUM,&info);
143 sprintf(errMsg,
"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetInfo() returned error code %d",__LINE__,ierr);
144 throw Comm_->ReportError(errMsg,-1);
147 ierr = MatGetSize(Amat,&NumGlobalRows_,&tmp);
149 DomainMap_ =
new Epetra_Map(NumGlobalRows_, NumMyRows_, MyGlobalElements, 0, *Comm_);
154 int * ColGIDs =
new int[NumMyCols_];
155 for (
int i=0; i<numLocalCols; i++) ColGIDs[i] = MyGlobalElements[i];
156 for (
int i=numLocalCols; i<NumMyCols_; i++) ColGIDs[i] = aij->garray[i-numLocalCols];
158 ColMap_ =
new Epetra_Map(-1, NumMyCols_, ColGIDs, 0, *Comm_);
162 delete [] MyGlobalElements;
182 PetscErrorCode MatCreateColmap_MPIAIJ_Private(Mat);
189 PetscInt *gcols, *lcols, ierr;
196 ierr=MatGetRow(
Amat_, globalRow, &nz, (
const PetscInt**) &gcols, (
const PetscScalar **) &vals);CHKERRQ(ierr);
200 if (strcmp(
MatType_,MATMPIAIJ) == 0) {
201 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)
Amat_->data;
202 lcols = (PetscInt *) malloc(nz *
sizeof(
int));
205 ierr = MatCreateColmap_MPIAIJ_Private(
Amat_);CHKERRQ(ierr);
223 int offset =
Amat_->cmap->n-1;
225 for (
int i=0; i<nz; i++) {
226 if (gcols[i] >=
Amat_->cmap->rstart && gcols[i] <
Amat_->cmap->rend) {
227 lcols[i] = gcols[i] -
Amat_->cmap->rstart;
229# ifdef PETSC_USE_CTABLE
230 ierr = PetscTableFind(aij->colmap,gcols[i]+1,lcols+i);CHKERRQ(ierr);
231 lcols[i] = lcols[i] + offset;
233 lcols[i] = aij->colmap[gcols[i]] + offset;
242 if (NumEntries > Length)
return(-1);
243 for (
int i=0; i<NumEntries; i++) {
244 Indices[i] = lcols[i];
247 if (alloc) free(lcols);
248 MatRestoreRow(
Amat_,globalRow,&nz,(
const PetscInt**) &gcols, (
const PetscScalar **) &vals);
259 MatGetRow(
Amat_, globalRow, &NumEntries, PETSC_NULL, PETSC_NULL);
260 MatRestoreRow(
Amat_,globalRow,PETSC_NULL, PETSC_NULL, PETSC_NULL);
273 int ierr=VecCreate(
Comm_->Comm(),&petscDiag);CHKERRQ(ierr);
276 ierr = VecSetType(petscDiag,VECMPI);CHKERRQ(ierr);
278 VecSetType(petscDiag,VECSEQ);
281 MatGetDiagonal(
Amat_, petscDiag);
282 VecGetArray(petscDiag,&vals);
283 VecGetLocalSize(petscDiag,&length);
284 for (
int i=0; i<length; i++) Diagonal[i] = vals[i];
285 VecRestoreArray(petscDiag,&vals);
286 VecDestroy(&petscDiag);
317 for (
int i=0; i<NumVectors; i++) {
326 ierr = MatMult(
Amat_,petscX,petscY);CHKERRQ(ierr);
328 ierr = VecGetArray(petscY,&vals);CHKERRQ(ierr);
329 ierr = VecGetLocalSize(petscY,&length);CHKERRQ(ierr);
330 for (
int j=0; j<length; j++) yptrs[i][j] = vals[j];
331 ierr = VecRestoreArray(petscY,&vals);CHKERRQ(ierr);
334 VecDestroy(&petscX); VecDestroy(&petscY);
338 flops *= (double) NumVectors;
397 int NumEntries =
GetRow(i);
399 for (j=0; j < NumEntries; j++) scale += fabs(
Values_[j]);
400 if (scale<Epetra_MinDouble) {
401 if (scale==0.0) ierr = 1;
402 else if (ierr!=1) ierr = 2;
435 for (i=0; i <
NumMyCols_; i++) (*xp)[i] = 0.0;
438 int NumEntries =
GetRow(i);
449 double scale = (*xp)[i];
450 if (scale<Epetra_MinDouble) {
451 if (scale==0.0) ierr = 1;
452 else if (ierr!=1) ierr = 2;
456 (*xp)[i] = 1.0/scale;
477 MatDiagonalScale(
Amat_, petscX, PETSC_NULL);
479 ierr=VecDestroy(&petscX); CHKERRQ(ierr);
497 MatDiagonalScale(
Amat_, PETSC_NULL, petscX);
499 ierr=VecDestroy(&petscX); CHKERRQ(ierr);
#define EPETRA_CHK_ERR(a)
const double Epetra_MaxDouble
void UpdateFlops(int Flops_in) const
const Epetra_BlockMap & Map() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int ExtractView(double **A, int *MyLDA) const
int RightScale(const Epetra_Vector &x)
Scales the Epetra_PETScAIJMatrix on the right with a Epetra_Vector x.
double NormInf() const
Returns the infinity norm of the global matrix.
virtual ~Epetra_PETScAIJMatrix()
Epetra_PETScAIJMatrix Destructor.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
virtual const Epetra_Import * RowMatrixImporter() const
Returns the Epetra_Import object that contains the import operations for distributed operations.
int InvColSums(Epetra_Vector &x) const
Computes the sum of absolute values of the columns of the Epetra_PETScAIJMatrix, results returned in ...
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
int InvRowSums(Epetra_Vector &x) const
Computes the sum of absolute values of the rows of the Epetra_PETScAIJMatrix, results returned in x.
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
virtual void Print(std::ostream &os) const
Print method.
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_PETScAIJMatrix multiplied by a Epetra_MultiVector X in Y.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
const Epetra_Map & RowMatrixColMap() const
Returns the Column Map object needed for implementing Epetra_RowMatrix.
Epetra_PETScAIJMatrix(Mat Amat)
Epetra_PETScAIJMatrix constructor.
int MaxNumEntries() const
Returns the maximum of NumMyRowEntries() over all rows.
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_PETScAIJMatrix on the left with a Epetra_Vector x.
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_PETScAIJMatrix multiplied by a Epetra_MultiVector X in Y.
Epetra_Import * Importer_
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator (same as domain).
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int GetRow(int Row) const
Epetra_SerialComm * Comm_
Epetra_MultiVector * ImportVector_
double NormOne() const
Returns the one norm of the global matrix.
int ExtractView(double **V) const