52 #include "Epetra_Import.h" 53 #include "Epetra_Export.h" 54 #include "Epetra_Vector.h" 55 #include "Epetra_MultiVector.h" 61 Epetra_Object(
"Epetra::PETScAIJMatrix"),
72 NumGlobalNonzeros_(0),
83 PetscObjectGetComm( (PetscObject)Amat, &comm);
84 Comm_ =
new Epetra_MpiComm(comm);
86 Comm_ =
new Epetra_SerialComm();
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_);
160 Importer_ =
new Epetra_Import(*ColMap_, *DomainMap_);
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);
277 # else //TODO untested!! 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);
293 const Epetra_MultiVector& X,
294 Epetra_MultiVector& Y)
const 297 int NumVectors = X.NumVectors();
298 if (NumVectors!=Y.NumVectors()) EPETRA_CHK_ERR(-1);
302 X.ExtractView(&xptrs);
303 Y.ExtractView(&yptrs);
317 for (
int i=0; i<NumVectors; i++) {
319 ierr=VecCreateMPIWithArray(
Comm_->Comm(),1,X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
320 ierr=VecCreateMPIWithArray(
Comm_->Comm(),1, Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
321 # else //FIXME untested 322 ierr=VecCreateSeqWithArray(
Comm_->Comm(),1, X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
323 ierr=VecCreateSeqWithArray(
Comm_->Comm(),1, Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
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;
346 const Epetra_MultiVector& X,
347 Epetra_MultiVector& Y)
const 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;
403 x[i] = Epetra_MaxDouble;
409 EPETRA_CHK_ERR(ierr);
419 if (!
Filled()) EPETRA_CHK_ERR(-1);
423 Epetra_Vector * xp = 0;
424 Epetra_Vector * x_tmp = 0;
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;
453 (*xp)[i] = Epetra_MaxDouble;
456 (*xp)[i] = 1.0/scale;
459 EPETRA_CHK_ERR(ierr);
469 X.ExtractView(&xptr);
472 int ierr=VecCreateMPIWithArray(
Comm_->Comm(),1, X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
473 # else //FIXME untested 474 int ierr=VecCreateSeqWithArray(
Comm_->Comm(),1, X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
477 MatDiagonalScale(
Amat_, petscX, PETSC_NULL);
479 ierr=VecDestroy(&petscX); CHKERRQ(ierr);
489 X.ExtractView(&xptr);
492 int ierr=VecCreateMPIWithArray(
Comm_->Comm(),1,X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
493 # else //FIXME untested 494 int ierr=VecCreateSeqWithArray(
Comm_->Comm(),1,X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
497 MatDiagonalScale(
Amat_, PETSC_NULL, petscX);
499 ierr=VecDestroy(&petscX); CHKERRQ(ierr);
518 if (!
Filled()) EPETRA_CHK_ERR(-1);
Epetra_SerialComm * Comm_
Epetra_PETScAIJMatrix(Mat Amat)
Epetra_PETScAIJMatrix constructor.
virtual const Epetra_Import * RowMatrixImporter() const
Returns the Epetra_Import object that contains the import operations for distributed operations...
int MaxNumEntries() const
Returns the maximum of NumMyRowEntries() over all rows.
int GetRow(int Row) const
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.
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
int RightScale(const Epetra_Vector &x)
Scales the Epetra_PETScAIJMatrix on the right with a Epetra_Vector x.
double NormOne() const
Returns the one norm of the global matrix.
virtual ~Epetra_PETScAIJMatrix()
Epetra_PETScAIJMatrix Destructor.
Epetra_MultiVector * ImportVector_
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_PETScAIJMatrix on the left with a Epetra_Vector x.
Epetra_Import * Importer_
int InvColSums(Epetra_Vector &x) const
Computes the sum of absolute values of the columns of the Epetra_PETScAIJMatrix, results returned in ...
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
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 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 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...
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
double NormInf() const
Returns the infinity norm of the global matrix.
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator (same as domain)...
const Epetra_Map & RowMatrixColMap() const
Returns the Column Map object needed for implementing Epetra_RowMatrix.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
virtual void Print(std::ostream &os) const
Print method.