Ifpack_OverlappingRowMatrix.cpp

Go to the documentation of this file.
00001 #include "Ifpack_ConfigDefs.h"
00002 #include "Ifpack_OverlappingRowMatrix.h"
00003 #include "Epetra_RowMatrix.h"
00004 #include "Epetra_Map.h"
00005 #include "Epetra_CrsMatrix.h"
00006 #include "Epetra_Import.h"
00007 #include "Epetra_Comm.h"
00008 #include "Epetra_MultiVector.h"
00009 
00010 // ====================================================================== 
00011 Ifpack_OverlappingRowMatrix::
00012 Ifpack_OverlappingRowMatrix(const Epetra_RowMatrix* Matrix,
00013                             int OverlapLevel)  :
00014   Map_(0),
00015   Importer_(0),
00016   Matrix_(Matrix),
00017   ExtMatrix_(0),
00018   ExtMap_(0),
00019   ExtImporter_(0),
00020   OverlapLevel_(OverlapLevel)
00021 {
00022   // should not be here if no overlap
00023   if (OverlapLevel == 0)
00024     IFPACK_CHK_ERRV(-1);
00025 
00026   // nothing to do as well with one process
00027   if (Comm().NumProc() == 1)
00028     IFPACK_CHK_ERRV(-1);
00029   
00030   NumMyRowsA_ = A().NumMyRows();
00031 
00032   // construct the external matrix
00033   vector<int> ExtElements; 
00034 
00035   Epetra_Map* TmpMap = 0;
00036   Epetra_CrsMatrix* TmpMatrix = 0; 
00037   Epetra_Import* TmpImporter = 0; 
00038 
00039   // importing rows corresponding to elements that are 
00040   // in ColMap, but not in RowMap 
00041   const Epetra_Map *RowMap; 
00042   const Epetra_Map *ColMap; 
00043 
00044   for (int overlap = 0 ; overlap < OverlapLevel ; ++overlap) {
00045     if (TmpMatrix) {
00046       RowMap = &(TmpMatrix->RowMatrixRowMap()); 
00047       ColMap = &(TmpMatrix->RowMatrixColMap()); 
00048     }
00049     else {
00050       RowMap = &(A().RowMatrixRowMap()); 
00051       ColMap = &(A().RowMatrixColMap()); 
00052     }
00053 
00054     int size = ColMap->NumMyElements() - RowMap->NumMyElements(); 
00055     vector<int> list(size); 
00056 
00057     int count = 0; 
00058 
00059     // define the set of rows that are in ColMap but not in RowMap 
00060     for (int i = 0 ; i < ColMap->NumMyElements() ; ++i) { 
00061       int GID = ColMap->GID(i); 
00062       if (A().RowMatrixRowMap().LID(GID) == -1) { 
00063         vector<int>::iterator pos 
00064           = find(ExtElements.begin(),ExtElements.end(),GID); 
00065         if (pos == ExtElements.end()) { 
00066           ExtElements.push_back(GID);
00067           list[count] = GID; 
00068           ++count; 
00069         } 
00070       } 
00071     } 
00072 
00073     if (TmpMap)
00074       delete TmpMap;
00075     TmpMap = new Epetra_Map(-1,count, &list[0],0,Comm()); 
00076 
00077     if (TmpMatrix)
00078       delete TmpMatrix;
00079     TmpMatrix = new Epetra_CrsMatrix(Copy,*TmpMap,0); 
00080 
00081     if (TmpImporter)
00082       delete TmpImporter;
00083     TmpImporter = new Epetra_Import(*TmpMap,A().RowMatrixRowMap()); 
00084 
00085     TmpMatrix->Import(A(),*TmpImporter,Insert); 
00086     TmpMatrix->FillComplete(A().OperatorDomainMap(),*TmpMap); 
00087 
00088   }
00089 
00090   // clean up memory
00091   delete TmpMap;
00092   delete TmpMatrix;
00093   delete TmpImporter;
00094 
00095   // build the map containing all the nodes (original
00096   // matrix + extended matrix)
00097   vector<int> list(NumMyRowsA_ + ExtElements.size());
00098   for (int i = 0 ; i < NumMyRowsA_ ; ++i)
00099     list[i] = A().RowMatrixRowMap().GID(i);
00100   for (int i = 0 ; i < (int)ExtElements.size() ; ++i)
00101     list[i + NumMyRowsA_] = ExtElements[i];
00102 
00103   Map_ = new Epetra_Map(-1, NumMyRowsA_ + ExtElements.size(),
00104                         &list[0], 0, Comm());
00105   // now build the map corresponding to all the external nodes
00106   // (with respect to A().RowMatrixRowMap().
00107   ExtMap_ = new Epetra_Map(-1,ExtElements.size(),
00108                            &ExtElements[0],0,A().Comm());
00109   ExtMatrix_ = new Epetra_CrsMatrix(Copy,*ExtMap_,*Map_,0); 
00110 
00111   ExtImporter_ = new Epetra_Import(*ExtMap_,A().RowMatrixRowMap()); 
00112   ExtMatrix_->Import(A(),*ExtImporter_,Insert); 
00113   ExtMatrix_->FillComplete(A().OperatorDomainMap(),*Map_);
00114 
00115   Importer_ = new Epetra_Import(*Map_,A().RowMatrixRowMap());
00116 
00117   // fix indices for overlapping matrix
00118   NumMyRowsB_ = B().NumMyRows();
00119   NumMyRows_ = NumMyRowsA_ + NumMyRowsB_;
00120   NumMyCols_ = NumMyRows_;
00121 
00122   NumMyDiagonals_ = A().NumMyDiagonals() + B().NumMyDiagonals();
00123 
00124   NumMyNonzeros_ = A().NumMyNonzeros() + B().NumMyNonzeros();
00125   Comm().SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
00126   MaxNumEntries_ = A().MaxNumEntries();
00127 
00128   if (MaxNumEntries_ < B().MaxNumEntries())
00129     MaxNumEntries_ = B().MaxNumEntries();
00130 
00131 }
00132 
00133 // ======================================================================
00134 Ifpack_OverlappingRowMatrix::~Ifpack_OverlappingRowMatrix()
00135 {
00136   if (Map_)
00137     delete Map_;
00138   if (Importer_)
00139     delete Importer_;
00140   if (ExtMatrix_)
00141     delete ExtMatrix_;
00142   if (ExtMap_)
00143     delete ExtMap_;
00144   if (ExtImporter_)
00145     delete ExtImporter_;
00146 }
00147 
00148 // ======================================================================
00149 int Ifpack_OverlappingRowMatrix::
00150 NumMyRowEntries(int MyRow, int & NumEntries) const
00151 {
00152   if (MyRow < NumMyRowsA_)
00153     return(A().NumMyRowEntries(MyRow,NumEntries));
00154   else
00155     return(B().NumMyRowEntries(MyRow - NumMyRowsA_, NumEntries));
00156 }
00157 
00158 // ======================================================================
00159 int Ifpack_OverlappingRowMatrix::
00160 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, 
00161                  int * Indices) const
00162 {
00163   int ierr;
00164   if (MyRow < NumMyRowsA_)
00165     ierr = A().ExtractMyRowCopy(MyRow,Length,NumEntries,Values,Indices);
00166   else
00167     ierr = B().ExtractMyRowCopy(MyRow - NumMyRowsA_,Length,NumEntries,
00168                                 Values,Indices);
00169 
00170   IFPACK_RETURN(ierr);
00171 }
00172 
00173 // ======================================================================
00174 int Ifpack_OverlappingRowMatrix::
00175 ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00176 {
00177   IFPACK_CHK_ERR(-1);
00178 }
00179 
00180 
00181 // ======================================================================
00182 int Ifpack_OverlappingRowMatrix::
00183 Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00184 {
00185   int NumVectors = X.NumVectors();
00186   vector<int> Ind(MaxNumEntries_);
00187   vector<double> Val(MaxNumEntries_);
00188 
00189   Y.PutScalar(0.0);
00190 
00191   // matvec with A (local rows)
00192   for (int i = 0 ; i < NumMyRowsA_ ; ++i) {
00193     for (int k = 0 ; k < NumVectors ; ++k) {
00194       int Nnz;
00195       IFPACK_CHK_ERR(A().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
00196                                           &Val[0], &Ind[0]));
00197       for (int j = 0 ; j < Nnz ; ++j) {
00198         Y[k][i] += Val[j] * X[k][Ind[j]];
00199       }
00200     }
00201   }
00202 
00203   // matvec with B (overlapping rows)
00204   for (int i = 0 ; i < NumMyRowsB_ ; ++i) {
00205     for (int k = 0 ; k < NumVectors ; ++k) {
00206       int Nnz;
00207       IFPACK_CHK_ERR(B().ExtractMyRowCopy(i,MaxNumEntries_,Nnz, 
00208                                           &Val[0], &Ind[0]));
00209       for (int j = 0 ; j < Nnz ; ++j) {
00210         Y[k][i + NumMyRowsA_] += Val[j] * X[k][Ind[j]];
00211       }
00212     }
00213   }
00214   return(0);
00215 }
00216 
00217 // ======================================================================
00218 int Ifpack_OverlappingRowMatrix::
00219 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00220 {
00221   IFPACK_CHK_ERR(Multiply(UseTranspose(),X,Y));
00222   return(0);
00223 }
00224 
00225 // ======================================================================
00226 int Ifpack_OverlappingRowMatrix::
00227 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00228 {
00229   IFPACK_CHK_ERR(-1);
00230 }
00231 
00232 // ======================================================================
00233 Epetra_RowMatrix& Ifpack_OverlappingRowMatrix::B() const
00234 {
00235   return(*ExtMatrix_);
00236 }
00237 
00238 // ======================================================================
00239 const Epetra_BlockMap& Ifpack_OverlappingRowMatrix::Map() const
00240 {
00241   return(*Map_);
00242 }
00243 
00244 // ======================================================================
00245 int Ifpack_OverlappingRowMatrix::
00246 ImportMultiVector(const Epetra_MultiVector& X, Epetra_MultiVector& OvX,
00247                   Epetra_CombineMode CM)
00248 {
00249   OvX.Import(X,*Importer_,CM);
00250   return(0);
00251 }
00252 
00253 // ======================================================================
00254 int Ifpack_OverlappingRowMatrix::
00255 ExportMultiVector(const Epetra_MultiVector& OvX, Epetra_MultiVector& X,
00256                   Epetra_CombineMode CM)
00257 {
00258   X.Export(OvX,*Importer_,CM);
00259   return(0);
00260 }
00261 

Generated on Thu Sep 18 12:37:22 2008 for Ifpack Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1