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