Ifpack_OverlappingRowMatrix.cpp

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

Generated on Tue Oct 20 12:48:54 2009 for IFPACK by doxygen 1.4.7