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 
00166   // PETSc assumes the row number is global, whereas Trilinos assumes it's local.
00167   int globalRow = PetscRowStart_ + Row;
00168   assert(globalRow < PetscRowEnd_);
00169   ierr=MatGetRow(Amat_, globalRow, &nz, (const PetscInt**) &gcols, (const PetscScalar **) &vals);CHKERRQ(ierr);
00170 
00171   // I ripped this bit of code from PETSc's MatSetValues_MPIAIJ() in mpiaij.c.  The PETSc getrow returns everything in
00172   // global numbering, so we must convert to local numbering.
00173   if (strcmp(MatType_,MATMPIAIJ) == 0) {
00174     Mat_MPIAIJ  *aij = (Mat_MPIAIJ*)Amat_->data;
00175     lcols = (PetscInt *) malloc(nz * sizeof(int));
00176     if (!aij->colmap) {
00177       ierr = CreateColmap_MPIAIJ_Private(Amat_);CHKERRQ(ierr);
00178     }
00179     /*
00180       A PETSc parallel aij matrix uses two matrices to represent the local rows.
00181       The first matrix, A, is square and contains all local columns.
00182       The second matrix, B, is rectangular and contains all non-local columns.
00183 
00184       Matrix A:
00185       Local column ID's are mapped to global column id's by adding cmap.rstart.
00186       Given the global ID of a local column, the local ID is found by
00187       subtracting cmap.rstart.
00188 
00189       Matrix B:
00190       Non-local column ID's are mapped to global column id's by the local-to-
00191       global map garray.  Given the global ID of a local column, the local ID is
00192       found by the global-to-local map colmap.  colmap is either an array or
00193       hash table, the latter being the case when PETSC_USE_CTABLE is defined.
00194     */
00195     int offset = Amat_->cmap->n-1; //offset for non-local column indices
00196 
00197     for (int i=0; i<nz; i++) {
00198       if (gcols[i] >= Amat_->cmap->rstart && gcols[i] < Amat_->cmap->rend) {
00199         lcols[i] = gcols[i] - Amat_->cmap->rstart;
00200       } else {
00201 #       ifdef PETSC_USE_CTABLE
00202         ierr = PetscTableFind(aij->colmap,gcols[i]+1,lcols+i);CHKERRQ(ierr);
00203         lcols[i] = lcols[i] + offset;
00204 #       else
00205         lcols[i] = aij->colmap[gcols[i]] + offset;
00206 #       endif
00207       }
00208 
00209     } //for i=0; i<nz; i++)
00210   }
00211   else lcols = gcols;
00212 
00213   NumEntries = nz;
00214   if (NumEntries > Length) return(-1);
00215   for (int i=0; i<NumEntries; i++) {
00216     Indices[i] = lcols[i];
00217     Values[i] = vals[i];
00218   }
00219   MatRestoreRow(Amat_,globalRow,&nz,(const PetscInt**) &gcols, (const PetscScalar **) &vals);
00220   return(0);
00221 } //ExtractMyRowCopy()
00222 
00223 //==========================================================================
00224 
00225 int Epetra_PETScAIJMatrix::NumMyRowEntries(int Row, int & NumEntries) const 
00226 {
00227   int globalRow = PetscRowStart_ + Row;
00228   MatGetRow(Amat_, globalRow, &NumEntries, PETSC_NULL, PETSC_NULL);
00229   MatRestoreRow(Amat_,globalRow,&NumEntries, PETSC_NULL, PETSC_NULL);
00230   return(0);
00231 }
00232 
00233 //==============================================================================
00234 
00235 int Epetra_PETScAIJMatrix::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00236 {
00237 
00238   //TODO optimization: only get this diagonal once
00239   Vec petscDiag;
00240   double *vals=0;
00241   int length;
00242 
00243   int ierr=VecCreate(Comm_->Comm(),&petscDiag);CHKERRQ(ierr);
00244   VecSetSizes(petscDiag,NumMyRows_,NumGlobalRows_);
00245 # ifdef HAVE_MPI
00246   ierr = VecSetType(petscDiag,VECMPI);CHKERRQ(ierr);
00247 # else //TODO untested!!
00248   VecSetType(petscDiag,VECSEQ);
00249 # endif
00250 
00251   MatGetDiagonal(Amat_, petscDiag);
00252   VecGetArray(petscDiag,&vals);
00253   VecGetLocalSize(petscDiag,&length);
00254   for (int i=0; i<length; i++) Diagonal[i] = vals[i];
00255   VecRestoreArray(petscDiag,&vals);
00256   VecDestroy(petscDiag);
00257   return(0);
00258 }
00259 
00260 //=============================================================================
00261 
00262 int Epetra_PETScAIJMatrix::Multiply(bool TransA,
00263                                const Epetra_MultiVector& X,
00264                                Epetra_MultiVector& Y) const
00265 {
00266   (void)TransA;
00267   int NumVectors = X.NumVectors();
00268   if (NumVectors!=Y.NumVectors()) EPETRA_CHK_ERR(-1);  // X and Y must have same number of vectors
00269 
00270   double ** xptrs;
00271   double ** yptrs;
00272   X.ExtractView(&xptrs);
00273   Y.ExtractView(&yptrs);
00274   if (RowMatrixImporter()!=0) {
00275     if (ImportVector_!=0) {
00276       if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
00277     }
00278     if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(RowMatrixColMap(),NumVectors);
00279     ImportVector_->Import(X, *RowMatrixImporter(), Insert);
00280     ImportVector_->ExtractView(&xptrs);
00281   }
00282 
00283   double *vals=0;
00284   int length;
00285   Vec petscX, petscY;
00286   int ierr;
00287   for (int i=0; i<NumVectors; i++) {
00288 #   ifdef HAVE_MPI
00289     ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
00290     ierr=VecCreateMPIWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
00291 #   else //FIXME  untested
00292     ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
00293     ierr=VecCreateSeqWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
00294 #   endif
00295 
00296     ierr = MatMult(Amat_,petscX,petscY);CHKERRQ(ierr);
00297 
00298     ierr = VecGetArray(petscY,&vals);CHKERRQ(ierr);
00299     ierr = VecGetLocalSize(petscY,&length);CHKERRQ(ierr);
00300     for (int j=0; j<length; j++) yptrs[i][j] = vals[j];
00301     ierr = VecRestoreArray(petscY,&vals);CHKERRQ(ierr);
00302   }
00303 
00304   VecDestroy(petscX); VecDestroy(petscY);
00305   
00306   double flops = NumGlobalNonzeros();
00307   flops *= 2.0;
00308   flops *= (double) NumVectors;
00309   UpdateFlops(flops);
00310   return(0);
00311 } //Multiply()
00312 
00313 //=============================================================================
00314 
00315 int Epetra_PETScAIJMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, 
00316                             const Epetra_MultiVector& X,
00317                             Epetra_MultiVector& Y) const
00318 {
00319   (void)Upper;
00320   (void)Trans;
00321   (void)UnitDiagonal;
00322   (void)X;
00323   (void)Y;
00324   return(-1); // Not implemented
00325 }
00326 
00327 //=============================================================================
00328 // Utility routine to get the specified row and put it into Values_ and Indices_ arrays
00329 int Epetra_PETScAIJMatrix::MaxNumEntries() const {
00330   int NumEntries;
00331 
00332   if (MaxNumEntries_==-1) {
00333     for (int i=0; i < NumMyRows_; i++) {
00334       NumMyRowEntries(i, NumEntries);
00335       if (NumEntries>MaxNumEntries_) MaxNumEntries_ = NumEntries;
00336     }
00337   }
00338   return(MaxNumEntries_);
00339 }
00340 
00341 //=============================================================================
00342 // Utility routine to get the specified row and put it into Values_ and Indices_ arrays
00343 int Epetra_PETScAIJMatrix::GetRow(int Row) const {
00344 
00345   int NumEntries;
00346   int MaxNumEntries = Epetra_PETScAIJMatrix::MaxNumEntries();
00347 
00348   if (MaxNumEntries>0) {
00349     if (Values_==0) Values_ = new double[MaxNumEntries];
00350     if (Indices_==0) Indices_ = new int[MaxNumEntries];
00351   }
00352   Epetra_PETScAIJMatrix::ExtractMyRowCopy(Row, MaxNumEntries, NumEntries, Values_, Indices_);
00353   
00354   return(NumEntries);
00355 }
00356 
00357 //=============================================================================
00358 int Epetra_PETScAIJMatrix::InvRowSums(Epetra_Vector& x) const {
00359 //
00360 // Put inverse of the sum of absolute values of the ith row of A in x[i].
00361 //
00362 
00363   if (!OperatorRangeMap().SameAs(x.Map())) EPETRA_CHK_ERR(-2); // x must have the same distribution as the range of A
00364 
00365   int ierr = 0;
00366   int i, j;
00367   for (i=0; i < NumMyRows_; i++) {
00368     int NumEntries = GetRow(i); // Copies ith row of matrix into Values_ and Indices_
00369     double scale = 0.0;
00370     for (j=0; j < NumEntries; j++) scale += fabs(Values_[j]);
00371     if (scale<Epetra_MinDouble) {
00372       if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
00373       else if (ierr!=1) ierr = 2;
00374       x[i] = Epetra_MaxDouble;
00375     }
00376     else
00377       x[i] = 1.0/scale;
00378   }
00379   UpdateFlops(NumGlobalNonzeros());
00380   EPETRA_CHK_ERR(ierr);
00381   return(0);
00382 }
00383 
00384 //=============================================================================
00385 int Epetra_PETScAIJMatrix::InvColSums(Epetra_Vector& x) const {
00386 //
00387 // Put inverse of the sum of absolute values of the jth column of A in x[j].
00388 //
00389 
00390   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00391   if (!OperatorDomainMap().SameAs(x.Map())) EPETRA_CHK_ERR(-2); // x must have the same distribution as the domain of A
00392   
00393 
00394   Epetra_Vector * xp = 0;
00395   Epetra_Vector * x_tmp = 0;
00396   
00397 
00398   // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00399   if (RowMatrixImporter()!=0) {
00400     x_tmp = new Epetra_Vector(RowMatrixColMap()); // Create import vector if needed
00401     xp = x_tmp;
00402   }
00403   int ierr = 0;
00404   int i, j;
00405 
00406   for (i=0; i < NumMyCols_; i++) (*xp)[i] = 0.0;
00407 
00408   for (i=0; i < NumMyRows_; i++) {
00409     int NumEntries = GetRow(i);// Copies ith row of matrix into Values_ and Indices_
00410     for (j=0; j < NumEntries; j++) (*xp)[Indices_[j]] += fabs(Values_[j]);
00411   }
00412 
00413   if (RowMatrixImporter()!=0){
00414     x.Export(*x_tmp, *RowMatrixImporter(), Add); // Fill x with Values from import vector
00415     delete x_tmp;
00416     xp = &x;
00417   }
00418   // Invert values, don't allow them to get too large
00419   for (i=0; i < NumMyRows_; i++) {
00420     double scale = (*xp)[i];
00421     if (scale<Epetra_MinDouble) {
00422       if (scale==0.0) ierr = 1; // Set error to 1 to signal that zero rowsum found (supercedes ierr = 2)
00423       else if (ierr!=1) ierr = 2;
00424       (*xp)[i] = Epetra_MaxDouble;
00425     }
00426     else
00427       (*xp)[i] = 1.0/scale;
00428   }
00429   UpdateFlops(NumGlobalNonzeros());
00430   EPETRA_CHK_ERR(ierr);
00431   return(0);
00432 }
00433 
00434 //=============================================================================
00435 int Epetra_PETScAIJMatrix::LeftScale(const Epetra_Vector& X) {
00436 //
00437 // This function scales the ith row of A by x[i].
00438 //
00439   double *xptr;
00440   X.ExtractView(&xptr);
00441   Vec petscX;
00442 # ifdef HAVE_MPI
00443   int ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00444 # else //FIXME  untested
00445   int ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00446 # endif
00447 
00448   MatDiagonalScale(Amat_, petscX, PETSC_NULL);
00449 
00450   ierr=VecDestroy(petscX); CHKERRQ(ierr);
00451   return(0);
00452 } //LeftScale()
00453 
00454 //=============================================================================
00455 int Epetra_PETScAIJMatrix::RightScale(const Epetra_Vector& X) {
00456 //
00457 // This function scales the jth row of A by x[j].
00458 //
00459   double *xptr;
00460   X.ExtractView(&xptr);
00461   Vec petscX;
00462 # ifdef HAVE_MPI
00463   int ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00464 # else //FIXME  untested
00465   int ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
00466 # endif
00467 
00468   MatDiagonalScale(Amat_, PETSC_NULL, petscX);
00469 
00470   ierr=VecDestroy(petscX); CHKERRQ(ierr);
00471   return(0);
00472 } //RightScale()
00473 
00474 //=============================================================================
00475 
00476 double Epetra_PETScAIJMatrix::NormInf() const {
00477 
00478   if (NormInf_>-1.0) return(NormInf_);
00479 
00480   MatNorm(Amat_, NORM_INFINITY,&NormInf_);
00481   UpdateFlops(NumGlobalNonzeros());
00482 
00483   return(NormInf_);
00484 }
00485 
00486 //=============================================================================
00487 
00488 double Epetra_PETScAIJMatrix::NormOne() const {
00489 
00490   if (NormOne_>-1.0) return(NormOne_);
00491   if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.
00492 
00493   MatNorm(Amat_, NORM_1,&NormOne_);
00494   UpdateFlops(NumGlobalNonzeros());
00495 
00496   return(NormOne_);
00497 }
00498 
00499 //=============================================================================
00500 
00501 void Epetra_PETScAIJMatrix::Print(ostream& os) const {
00502   return;
00503 }
00504 
00505 #endif /*HAVE_PETSC*/

Generated on Wed May 12 21:24:46 2010 for EpetraExt by  doxygen 1.4.7