EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_PETScAIJMatrix.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 /*#############################################################################
00030 # CVS File Information
00031 #    Current revision: $Revision$
00032 #    Last modified:    $Date$
00033 #    Modified by:      $Author$
00034 #############################################################################*/
00035 
00036 #include "EpetraExt_ConfigDefs.h"
00037 #ifdef HAVE_PETSC
00038 
00039 #include "Epetra_Import.h"
00040 #include "Epetra_Export.h"
00041 #include "Epetra_Vector.h"
00042 #include "Epetra_MultiVector.h"
00043 #include "EpetraExt_PETScAIJMatrix.h"
00044 
00045 //==============================================================================
00046 Epetra_PETScAIJMatrix::Epetra_PETScAIJMatrix(Mat Amat)
00047   : Epetra_Object("Epetra::PETScAIJMatrix"),
00048     Amat_(Amat),
00049     Values_(0),
00050     Indices_(0),
00051     MaxNumEntries_(-1),
00052     ImportVector_(0),
00053     NormInf_(-1.0),
00054     NormOne_(-1.0)
00055 {
00056 #ifdef HAVE_MPI
00057   MPI_Comm comm;
00058   PetscObjectGetComm( (PetscObject)Amat, &comm);
00059   Comm_ = new Epetra_MpiComm(comm);
00060 #else
00061   Comm_ = new Epetra_SerialComm();
00062 #endif  
00063   int ierr;
00064   char errMsg[80];
00065   MatGetType(Amat, &MatType_);
00066   if ( strcmp(MatType_,MATSEQAIJ) != 0 && strcmp(MatType_,MATMPIAIJ) != 0 ) {
00067     sprintf(errMsg,"PETSc matrix must be either seqaij or mpiaij (but it is %s)",MatType_);
00068     throw Comm_->ReportError(errMsg,-1);
00069   }
00070   petscMatrixType mt;
00071   Mat_MPIAIJ* aij=0;
00072   if (strcmp(MatType_,MATMPIAIJ) == 0) {
00073     mt = PETSC_MPI_AIJ;
00074     aij = (Mat_MPIAIJ*)Amat->data;
00075   }
00076   else if (strcmp(MatType_,MATSEQAIJ) == 0) {
00077     mt = PETSC_SEQ_AIJ;
00078   }
00079   int numLocalRows, numLocalCols;
00080   ierr = MatGetLocalSize(Amat,&numLocalRows,&numLocalCols);
00081   if (ierr) {
00082     sprintf(errMsg,"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetLocalSize() returned error code %d",__LINE__,ierr);
00083     throw Comm_->ReportError(errMsg,-1);
00084   }
00085   NumMyRows_ = numLocalRows;
00086   NumMyCols_ = numLocalCols; //numLocalCols is the total # of unique columns in the local matrix (the diagonal block)
00087   //TODO what happens if some columns are empty?
00088   if (mt == PETSC_MPI_AIJ)
00089     NumMyCols_ += aij->B->cmap->n;
00090   MatInfo info;
00091   ierr = MatGetInfo(Amat,MAT_LOCAL,&info);
00092   if (ierr) {
00093     sprintf(errMsg,"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetInfo() returned error code %d",__LINE__,ierr);
00094     throw Comm_->ReportError(errMsg,-1);
00095   }
00096   NumMyNonzeros_ = (int) info.nz_used; //PETSc stores nnz as double
00097   Comm_->SumAll(&(info.nz_used), &NumGlobalNonzeros_, 1);
00098 
00099   //The PETSc documentation warns that this may not be robust.
00100   //In particular, this will break if the ordering is not contiguous!
00101   int rowStart, rowEnd;
00102   ierr = MatGetOwnershipRange(Amat,&rowStart,&rowEnd);
00103   if (ierr) {
00104     sprintf(errMsg,"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetOwnershipRange() returned error code %d",__LINE__,ierr);
00105     throw Comm_->ReportError(errMsg,-1);
00106   }
00107 
00108   PetscRowStart_ = rowStart;
00109   PetscRowEnd_   = rowEnd;
00110 
00111   int* MyGlobalElements = new int[rowEnd-rowStart];
00112   for (int i=0; i<rowEnd-rowStart; i++)
00113     MyGlobalElements[i] = rowStart+i;
00114 
00115   ierr = MatGetInfo(Amat,MAT_GLOBAL_SUM,&info);
00116   if (ierr) {
00117     sprintf(errMsg,"EpetraExt_PETScAIJMatrix.cpp, line %d, MatGetInfo() returned error code %d",__LINE__,ierr);
00118     throw Comm_->ReportError(errMsg,-1);
00119   }
00120   int tmp;
00121   ierr = MatGetSize(Amat,&NumGlobalRows_,&tmp);
00122 
00123   DomainMap_ = new Epetra_Map(NumGlobalRows_, NumMyRows_, MyGlobalElements, 0, *Comm_);
00124 
00125   // get the GIDs of the non-local columns
00126   //FIXME what if the matrix is sequential?
00127 
00128   int * ColGIDs = new int[NumMyCols_];
00129   for (int i=0; i<numLocalCols; i++) ColGIDs[i] = MyGlobalElements[i];
00130   for (int i=numLocalCols; i<NumMyCols_; i++) ColGIDs[i] = aij->garray[i-numLocalCols];
00131 
00132   ColMap_ = new Epetra_Map(-1, NumMyCols_, ColGIDs, 0, *Comm_);
00133 
00134   Importer_ = new Epetra_Import(*ColMap_, *DomainMap_);
00135 
00136   delete [] MyGlobalElements;
00137   delete [] ColGIDs;
00138 } //Epetra_PETScAIJMatrix(Mat Amat)
00139 
00140 //==============================================================================
00141 
00142 Epetra_PETScAIJMatrix::~Epetra_PETScAIJMatrix(){
00143   if (ImportVector_!=0) delete ImportVector_;
00144   delete Importer_;
00145   delete ColMap_;
00146   delete DomainMap_;
00147   delete Comm_;
00148 
00149   if (Values_!=0) {delete [] Values_; Values_=0;}
00150   if (Indices_!=0) {delete [] Indices_; Indices_=0;}
00151 } //Epetra_PETScAIJMatrix dtor
00152 
00153 //==========================================================================
00154 
00155 extern "C" {
00156   PetscErrorCode CreateColmap_MPIAIJ_Private(Mat);
00157 }
00158 
00159 int Epetra_PETScAIJMatrix::ExtractMyRowCopy(int Row, int Length, int & NumEntries, double * Values,
00160                      int * Indices) const 
00161 {
00162   int nz;
00163   PetscInt *gcols, *lcols, ierr;
00164   PetscScalar *vals;
00165   bool alloc=false;
00166 
00167   // PETSc assumes the row number is global, whereas Trilinos assumes it's local.
00168   int globalRow = PetscRowStart_ + Row;
00169   assert(globalRow < PetscRowEnd_);
00170   ierr=MatGetRow(Amat_, globalRow, &nz, (const PetscInt**) &gcols, (const PetscScalar **) &vals);CHKERRQ(ierr);
00171 
00172   // I ripped this bit of code from PETSc's MatSetValues_MPIAIJ() in mpiaij.c.  The PETSc getrow returns everything in
00173   // global numbering, so we must convert to local numbering.
00174   if (strcmp(MatType_,MATMPIAIJ) == 0) {
00175     Mat_MPIAIJ  *aij = (Mat_MPIAIJ*)Amat_->data;
00176     lcols = (PetscInt *) malloc(nz * sizeof(int));
00177     alloc=true;
00178     if (!aij->colmap) {
00179       ierr = CreateColmap_MPIAIJ_Private(Amat_);CHKERRQ(ierr);
00180     }
00181     /*
00182       A PETSc parallel aij matrix uses two matrices to represent the local rows.
00183       The first matrix, A, is square and contains all local columns.
00184       The second matrix, B, is rectangular and contains all non-local columns.
00185 
00186       Matrix A:
00187       Local column ID's are mapped to global column id's by adding cmap.rstart.
00188       Given the global ID of a local column, the local ID is found by
00189       subtracting cmap.rstart.
00190 
00191       Matrix B:
00192       Non-local column ID's are mapped to global column id's by the local-to-
00193       global map garray.  Given the global ID of a local column, the local ID is
00194       found by the global-to-local map colmap.  colmap is either an array or
00195       hash table, the latter being the case when PETSC_USE_CTABLE is defined.
00196     */
00197     int offset = Amat_->cmap->n-1; //offset for non-local column indices
00198 
00199     for (int i=0; i<nz; i++) {
00200       if (gcols[i] >= Amat_->cmap->rstart && gcols[i] < Amat_->cmap->rend) {
00201         lcols[i] = gcols[i] - Amat_->cmap->rstart;
00202       } else {
00203 #       ifdef PETSC_USE_CTABLE
00204         ierr = PetscTableFind(aij->colmap,gcols[i]+1,lcols+i);CHKERRQ(ierr);
00205         lcols[i] = lcols[i] + offset;
00206 #       else
00207         lcols[i] = aij->colmap[gcols[i]] + offset;
00208 #       endif
00209       }
00210 
00211     } //for i=0; i<nz; i++)
00212   }
00213   else lcols = gcols;
00214 
00215   NumEntries = nz;
00216   if (NumEntries > Length) return(-1);
00217   for (int i=0; i<NumEntries; i++) {
00218     Indices[i] = lcols[i];
00219     Values[i] = vals[i];
00220   }
00221   if (alloc) free(lcols);
00222   MatRestoreRow(Amat_,globalRow,&nz,(const PetscInt**) &gcols, (const PetscScalar **) &vals);
00223   return(0);
00224 } //ExtractMyRowCopy()
00225 
00226 //==========================================================================
00227 
00228 int Epetra_PETScAIJMatrix::NumMyRowEntries(int Row, int & NumEntries) const 
00229 {
00230   int globalRow = PetscRowStart_ + Row;
00231   MatGetRow(Amat_, globalRow, &NumEntries, PETSC_NULL, PETSC_NULL);
00232   MatRestoreRow(Amat_,globalRow,&NumEntries, PETSC_NULL, PETSC_NULL);
00233   return(0);
00234 }
00235 
00236 //==============================================================================
00237 
00238 int Epetra_PETScAIJMatrix::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00239 {
00240 
00241   //TODO optimization: only get this diagonal once
00242   Vec petscDiag;
00243   double *vals=0;
00244   int length;
00245 
00246   int ierr=VecCreate(Comm_->Comm(),&petscDiag);CHKERRQ(ierr);
00247   VecSetSizes(petscDiag,NumMyRows_,NumGlobalRows_);
00248 # ifdef HAVE_MPI
00249   ierr = VecSetType(petscDiag,VECMPI);CHKERRQ(ierr);
00250 # else //TODO untested!!
00251   VecSetType(petscDiag,VECSEQ);
00252 # endif
00253 
00254   MatGetDiagonal(Amat_, petscDiag);
00255   VecGetArray(petscDiag,&vals);
00256   VecGetLocalSize(petscDiag,&length);
00257   for (int i=0; i<length; i++) Diagonal[i] = vals[i];
00258   VecRestoreArray(petscDiag,&vals);
00259   VecDestroy(petscDiag);
00260   return(0);
00261 }
00262 
00263 //=============================================================================
00264 
00265 int Epetra_PETScAIJMatrix::Multiply(bool TransA,
00266                                const Epetra_MultiVector& X,
00267                                Epetra_MultiVector& Y) const
00268 {
00269   (void)TransA;
00270   int NumVectors = X.NumVectors();
00271   if (NumVectors!=Y.NumVectors()) EPETRA_CHK_ERR(-1);  // X and Y must have same number of vectors
00272 
00273   double ** xptrs;
00274   double ** yptrs;
00275   X.ExtractView(&xptrs);
00276   Y.ExtractView(&yptrs);
00277   if (RowMatrixImporter()!=0) {
00278     if (ImportVector_!=0) {
00279       if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
00280     }
00281     if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(RowMatrixColMap(),NumVectors);
00282     ImportVector_->Import(X, *RowMatrixImporter(), Insert);
00283     ImportVector_->ExtractView(&xptrs);
00284   }
00285 
00286   double *vals=0;
00287   int length;
00288   Vec petscX, petscY;
00289   int ierr;
00290   for (int i=0; i<NumVectors; i++) {
00291 #   ifdef HAVE_MPI
00292     ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
00293     ierr=VecCreateMPIWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
00294 #   else //FIXME  untested
00295     ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
00296     ierr=VecCreateSeqWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
00297 #   endif
00298 
00299     ierr = MatMult(Amat_,petscX,petscY);CHKERRQ(ierr);
00300 
00301     ierr = VecGetArray(petscY,&vals);CHKERRQ(ierr);
00302     ierr = VecGetLocalSize(petscY,&length);CHKERRQ(ierr);
00303     for (int j=0; j<length; j++) yptrs[i][j] = vals[j];
00304     ierr = VecRestoreArray(petscY,&vals);CHKERRQ(ierr);
00305   }
00306 
00307   VecDestroy(petscX); VecDestroy(petscY);
00308   
00309   double flops = NumGlobalNonzeros();
00310   flops *= 2.0;
00311   flops *= (double) NumVectors;
00312   UpdateFlops(flops);
00313   return(0);
00314 } //Multiply()
00315 
00316 //=============================================================================
00317 
00318 int Epetra_PETScAIJMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00319                             const Epetra_MultiVector& X,
00320                             Epetra_MultiVector& Y) const
00321 {
00322   (void)Upper;
00323   (void)Trans;
00324   (void)UnitDiagonal;
00325   (void)X;
00326   (void)Y;
00327   return(-1); // Not implemented
00328 }
00329 
00330 //=============================================================================
00331 // Utility routine to get the specified row and put it into Values_ and Indices_ arrays
00332 int Epetra_PETScAIJMatrix::MaxNumEntries() const {
00333   int NumEntries;
00334 
00335   if (MaxNumEntries_==-1) {
00336     for (int i=0; i < NumMyRows_; i++) {
00337       NumMyRowEntries(i, NumEntries);
00338       if (NumEntries>MaxNumEntries_) MaxNumEntries_ = NumEntries;
00339     }
00340   }
00341   return(MaxNumEntries_);
00342 }
00343 
00344 //=============================================================================
00345 // Utility routine to get the specified row and put it into Values_ and Indices_ arrays
00346 int Epetra_PETScAIJMatrix::GetRow(int Row) const {
00347 
00348   int NumEntries;
00349   int MaxNumEntries = Epetra_PETScAIJMatrix::MaxNumEntries();
00350 
00351   if (MaxNumEntries>0) {
00352     if (Values_==0) Values_ = new double[MaxNumEntries];
00353     if (Indices_==0) Indices_ = new int[MaxNumEntries];
00354   }
00355   Epetra_PETScAIJMatrix::ExtractMyRowCopy(Row, MaxNumEntries, NumEntries, Values_, Indices_);
00356   
00357   return(NumEntries);
00358 }
00359 
00360 //=============================================================================
00361 int Epetra_PETScAIJMatrix::InvRowSums(Epetra_Vector& x) const {
00362 //
00363 // Put inverse of the sum of absolute values of the ith row of A in x[i].
00364 //
00365 
00366   if (!OperatorRangeMap().SameAs(x.Map())) EPETRA_CHK_ERR(-2); // x must have the same distribution as the range of A
00367 
00368   int ierr = 0;
00369   int i, j;
00370   for (i=0; i < NumMyRows_; i++) {
00371     int NumEntries = GetRow(i); // Copies ith row of matrix into Values_ and Indices_
00372     double scale = 0.0;
00373     for (j=0; j < NumEntries; j++) scale += fabs(Values_[j]);
00374     if (scale<Epetra_MinDouble) {
00375       if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
00376       else if (ierr!=1) ierr = 2;
00377       x[i] = Epetra_MaxDouble;
00378     }
00379     else
00380       x[i] = 1.0/scale;
00381   }
00382   UpdateFlops(NumGlobalNonzeros());
00383   EPETRA_CHK_ERR(ierr);
00384   return(0);
00385 }
00386 
00387 //=============================================================================
00388 int Epetra_PETScAIJMatrix::InvColSums(Epetra_Vector& x) const {
00389 //
00390 // Put inverse of the sum of absolute values of the jth column of A in x[j].
00391 //
00392 
00393   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00394   if (!OperatorDomainMap().SameAs(x.Map())) EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
00395   
00396 
00397   Epetra_Vector * xp = 0;
00398   Epetra_Vector * x_tmp = 0;
00399   
00400 
00401   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00402   if (RowMatrixImporter()!=0) {
00403     x_tmp = new Epetra_Vector(RowMatrixColMap()); // Create import vector if needed
00404     xp = x_tmp;
00405   }
00406   int ierr = 0;
00407   int i, j;
00408 
00409   for (i=0; i < NumMyCols_; i++) (*xp)[i] = 0.0;
00410 
00411   for (i=0; i < NumMyRows_; i++) {
00412     int NumEntries = GetRow(i);// Copies ith row of matrix into Values_ and Indices_
00413     for (j=0; j < NumEntries; j++) (*xp)[Indices_[j]] += fabs(Values_[j]);
00414   }
00415 
00416   if (RowMatrixImporter()!=0){
00417     x.Export(*x_tmp, *RowMatrixImporter(), Add); // Fill x with Values from import vector
00418     delete x_tmp;
00419     xp = &x;
00420   }
00421   // Invert values, don't allow them to get too large
00422   for (i=0; i < NumMyRows_; i++) {
00423     double scale = (*xp)[i];
00424     if (scale<Epetra_MinDouble) {
00425       if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
00426       else if (ierr!=1) ierr = 2;
00427       (*xp)[i] = Epetra_MaxDouble;
00428     }
00429     else
00430       (*xp)[i] = 1.0/scale;
00431   }
00432   UpdateFlops(NumGlobalNonzeros());
00433   EPETRA_CHK_ERR(ierr);
00434   return(0);
00435 }
00436 
00437 //=============================================================================
00438 int Epetra_PETScAIJMatrix::LeftScale(const Epetra_Vector& X) {
00439 //
00440 // This function scales the ith row of A by x[i].
00441 //
00442   double *xptr;
00443   X.ExtractView(&xptr);
00444   Vec petscX;
00445 # ifdef HAVE_MPI
00446   int ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00447 # else //FIXME  untested
00448   int ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00449 # endif
00450 
00451   MatDiagonalScale(Amat_, petscX, PETSC_NULL);
00452 
00453   ierr=VecDestroy(petscX); CHKERRQ(ierr);
00454   return(0);
00455 } //LeftScale()
00456 
00457 //=============================================================================
00458 int Epetra_PETScAIJMatrix::RightScale(const Epetra_Vector& X) {
00459 //
00460 // This function scales the jth row of A by x[j].
00461 //
00462   double *xptr;
00463   X.ExtractView(&xptr);
00464   Vec petscX;
00465 # ifdef HAVE_MPI
00466   int ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00467 # else //FIXME  untested
00468   int ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00469 # endif
00470 
00471   MatDiagonalScale(Amat_, PETSC_NULL, petscX);
00472 
00473   ierr=VecDestroy(petscX); CHKERRQ(ierr);
00474   return(0);
00475 } //RightScale()
00476 
00477 //=============================================================================
00478 
00479 double Epetra_PETScAIJMatrix::NormInf() const {
00480 
00481   if (NormInf_>-1.0) return(NormInf_);
00482 
00483   MatNorm(Amat_, NORM_INFINITY,&NormInf_);
00484   UpdateFlops(NumGlobalNonzeros());
00485 
00486   return(NormInf_);
00487 }
00488 
00489 //=============================================================================
00490 
00491 double Epetra_PETScAIJMatrix::NormOne() const {
00492 
00493   if (NormOne_>-1.0) return(NormOne_);
00494   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00495 
00496   MatNorm(Amat_, NORM_1,&NormOne_);
00497   UpdateFlops(NumGlobalNonzeros());
00498 
00499   return(NormOne_);
00500 }
00501 
00502 //=============================================================================
00503 
00504 void Epetra_PETScAIJMatrix::Print(ostream& os) const {
00505   return;
00506 }
00507 
00508 #endif /*HAVE_PETSC*/
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines