Ifpack_NodeFilter.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 #ifdef IFPACK_NODE_AWARE_CODE
00031 
00032 #include "Ifpack_ConfigDefs.h"
00033 
00034 #include "Epetra_MultiVector.h"
00035 #include "Epetra_Vector.h"
00036 #include "Epetra_IntVector.h"
00037 #include "Epetra_RowMatrix.h"
00038 #include "Epetra_Map.h"
00039 #include "Epetra_BlockMap.h"
00040 #include "Ifpack_NodeFilter.h"
00041 #include "Ifpack_OverlappingRowMatrix.h"
00042 #include "Epetra_Import.h"
00043 #include "Epetra_Export.h"
00044 #include "Epetra_CrsMatrix.h"
00045 #include "Epetra_BLAS_wrappers.h"
00046 
00047 using namespace Teuchos;
00048 
00049 #define OLD_AND_BUSTED
00050 
00051 //==============================================================================
00052 Ifpack_NodeFilter::Ifpack_NodeFilter(const RefCountPtr<const Epetra_RowMatrix>& Matrix,int nodeID) :
00053   Matrix_(Matrix),
00054   NumMyRows_(0),
00055   NumMyNonzeros_(0),
00056   NumGlobalRows_(0),
00057   MaxNumEntries_(0),
00058   MaxNumEntriesA_(0)
00059 {
00060   sprintf(Label_,"%s","Ifpack_NodeFilter");
00061 
00062   //ImportVector_=null;
00063   //ExportVector_=null;
00064   ImportVector_=0;
00065   ExportVector_=0;
00066 
00067   ovA_ = dynamic_cast<const Ifpack_OverlappingRowMatrix*>(&*Matrix_);
00068   //assert(ovA_ != 0);
00069 
00070   if (ovA_) {
00071     Acrs_=dynamic_cast<const Epetra_CrsMatrix*>(&ovA_->A());
00072     NumMyRowsA_ = ovA_->A().NumMyRows();
00073     NumMyRowsB_ = ovA_->B().NumMyRows();
00074   } else {
00075     NumMyRowsA_ = Matrix->NumMyRows();
00076     NumMyRowsB_ = 0;
00077   }
00078   
00079 #ifdef HAVE_MPI
00080   const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &(Matrix->Comm()) );
00081   assert(pComm != NULL);
00082   MPI_Comm_split(pComm->Comm(),nodeID,pComm->MyPID(),&nodeMPIComm_);
00083   SubComm_ = rcp( new Epetra_MpiComm(nodeMPIComm_) );
00084 #else
00085   SubComm_ = rcp( new Epetra_SerialComm );
00086 #endif
00087 
00088   NumMyRows_ = Matrix->NumMyRows();
00089   SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
00090   NumMyCols_ = Matrix->NumMyCols();
00091 
00092   // build row map, based on the local communicator
00093   try {
00094     const Epetra_Map &globRowMap = Matrix->RowMatrixRowMap();
00095     int *myGlobalElts =  globRowMap.MyGlobalElements();
00096     int numMyElts = globRowMap.NumMyElements();
00097     Map_ = rcp( new Epetra_Map(-1,numMyElts,myGlobalElts,globRowMap.IndexBase(),*SubComm_) );
00098   }
00099   catch(...) {
00100     printf("** * Ifpack_NodeFilter ctor: problem creating row map * **\n\n");
00101   }
00102 
00103 
00104   // build the column map, but don't use a copy constructor, b/c local communicator SubComm_ is
00105   // different from that of Matrix.
00106   try {
00107     const Epetra_Map &globColMap = Matrix->RowMatrixColMap();
00108     int *myGlobalElts =  globColMap.MyGlobalElements();
00109     int numMyElts = globColMap.NumMyElements();
00110     colMap_ = rcp( new Epetra_Map(-1,numMyElts,myGlobalElts,globColMap.IndexBase(),*SubComm_) );
00111   }
00112   catch(...) {
00113     printf("** * Ifpack_NodeFilter ctor: problem creating col map * **\n\n");
00114   }
00115 
00116   // NumEntries_ will contain the actual number of nonzeros
00117   // for each localized row (that is, without external nodes,
00118   // and always with the diagonal entry)
00119   NumEntries_.resize(NumMyRows_);
00120 
00121   // want to store the diagonal vector. FIXME: am I really useful?
00122   Diagonal_ = rcp( new Epetra_Vector(*Map_) );
00123   if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
00124 
00125   // store this for future access to ExtractMyRowCopy().
00126   // This is the # of nonzeros in the non-local matrix
00127   MaxNumEntriesA_ = Matrix->MaxNumEntries();
00128   // tentative value for MaxNumEntries. This is the number of
00129   // nonzeros in the local matrix
00130   MaxNumEntries_ = Matrix->MaxNumEntries();
00131 
00132   // ExtractMyRowCopy() will use these vectors
00133   Indices_.resize(MaxNumEntries_);
00134   Values_.resize(MaxNumEntries_);
00135 
00136   // now compute:
00137   // - the number of nonzero per row
00138   // - the total number of nonzeros
00139   // - the diagonal entries
00140 
00141 
00142   Ac_LIDMap_ = 0;
00143   Bc_LIDMap_ = 0;
00144   Ar_LIDMap_ = 0;
00145   Br_LIDMap_ = 0;
00146   tempX_ = 0;
00147   tempY_ = 0;
00148   // CMS: [A|B]-Local to Overlap-Local Column Indices
00149   if(ovA_){
00150     Ac_LIDMap_=new int[ovA_->A().NumMyCols()+1];
00151     for(int i=0;i<ovA_->A().NumMyCols();i++) Ac_LIDMap_[i]=colMap_->LID(ovA_->A().RowMatrixColMap().GID(i));    
00152     Bc_LIDMap_=new int[ovA_->B().NumMyCols()+1];
00153     for(int i=0;i<ovA_->B().NumMyCols();i++) Bc_LIDMap_[i]=colMap_->LID(ovA_->B().RowMatrixColMap().GID(i));
00154 
00155     Ar_LIDMap_=new int[ovA_->A().NumMyRows()+1];
00156     for(int i=0;i<ovA_->A().NumMyRows();i++) Ar_LIDMap_[i]=Map_->LID(ovA_->A().RowMatrixRowMap().GID(i));    
00157     Br_LIDMap_=new int[ovA_->B().NumMyRows()+1];
00158     for(int i=0;i<ovA_->B().NumMyRows();i++) Br_LIDMap_[i]=Map_->LID(ovA_->B().RowMatrixRowMap().GID(i));
00159 
00160 
00161 #ifndef OLD_AND_BUSTED
00162     NumMyColsA_=ovA_->A().NumMyCols();
00163     tempX_=new double[NumMyColsA_];
00164     tempY_=new double[NumMyRowsA_];     
00165 #else
00166     tempX_=0;
00167     tempY_=0;
00168 #endif
00169   }
00170   // end CMS
00171 
00172 
00173   // compute nonzeros (total and per-row), and store the
00174   // diagonal entries (already modified)
00175   int ActualMaxNumEntries = 0;
00176 
00177   for (int i = 0 ; i < NumMyRows_ ; ++i) {
00178     
00179     NumEntries_[i] = 0;
00180     int Nnz, NewNnz = 0;
00181     IFPACK_CHK_ERRV(ExtractMyRowCopy(i,MaxNumEntries_,Nnz,&Values_[0],&Indices_[0]));
00182 
00183     for (int j = 0 ; j < Nnz ; ++j) {
00184       NewNnz++;
00185       if (Indices_[j] == i) (*Diagonal_)[i] = Values_[j];
00186     }
00187 
00188     if (NewNnz > ActualMaxNumEntries)
00189       ActualMaxNumEntries = NewNnz;
00190 
00191     NumMyNonzeros_ += NewNnz;
00192     NumEntries_[i] = NewNnz;
00193 
00194   }
00195 
00196   SubComm_->SumAll(&NumMyNonzeros_,&NumGlobalNonzeros_,1);
00197   MaxNumEntries_ = ActualMaxNumEntries;
00198 
00199   int gpid = Matrix->Comm().MyPID();
00200   Exporter_ = null;
00201   Importer_ = null;
00202   // Check if non-trivial import/export operators
00203   if (!(RowMatrixRowMap().SameAs(OperatorRangeMap()))) {
00204     try{Exporter_ = rcp(new Epetra_Export(RowMatrixRowMap(), OperatorRangeMap()));}
00205     catch(...) {
00206       printf("** * gpid %d: Ifpack_NodeFilter ctor: problem creating Exporter_ * **\n\n",gpid);
00207     }
00208   }
00209   if (!(*colMap_).SameAs(*Map_)) {
00210     //TODO change this to RCP
00211     try{Importer_ = rcp(new Epetra_Import(*colMap_, *Map_));}
00212     catch(...) {
00213       printf("** * gpid %d: Ifpack_NodeFilter ctor: problem creating Importer_ * **\n\n",gpid);
00214     }
00215   }
00216 
00217 } //Ifpack_NodeFilter() ctor
00218 
00219 //==============================================================================
00220 int Ifpack_NodeFilter::
00221 ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00222          double *Values, int * Indices) const
00223 {
00224   if ((MyRow < 0) || (MyRow >= NumMyRows_)) {
00225     IFPACK_CHK_ERR(-1); // range not valid
00226   }
00227 
00228   if (Length < NumEntries_[MyRow])
00229     assert(1==0);
00230     //return(-1);
00231 
00232   //TODO  will we have problems in the original use case?
00233   // always extract using the object Values_ and Indices_.
00234   // This is because I need more space than that given by
00235   // the user (for the external nodes)
00236 
00237   int ierr;
00238   if (ovA_) {
00239     int LocRow=ovA_->A().RowMatrixRowMap().LID(Map_->GID(MyRow));      
00240     if (LocRow != -1) { 
00241       ierr=ovA_->A().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
00242       for(int j=0;j<NumEntries;j++)
00243     Indices[j]=Ac_LIDMap_[Indices[j]];
00244       
00245     }
00246     else {
00247       LocRow=ovA_->B().RowMatrixRowMap().LID(Map_->GID(MyRow));      
00248       ierr=ovA_->B().ExtractMyRowCopy(LocRow,Length,NumEntries,Values,Indices);
00249       for(int j=0;j<NumEntries;j++)
00250     Indices[j]=Bc_LIDMap_[Indices[j]];
00251     }
00252   } else {
00253     int Nnz;
00254     ierr = Matrix_->ExtractMyRowCopy(MyRow,MaxNumEntriesA_,Nnz, &Values_[0],&Indices_[0]);
00255     IFPACK_CHK_ERR(ierr);
00256 
00257     NumEntries = 0;
00258     for (int j = 0 ; j < Nnz ; ++j) {
00259       // only local indices
00260       if (Indices_[j] < NumMyRows_ ) {
00261         Indices[NumEntries] = Indices_[j];
00262         Values[NumEntries] = Values_[j];
00263         ++NumEntries;
00264       }
00265     }
00266   }
00267     
00268   return(ierr);
00269 
00270 }
00271 
00272 //==============================================================================
00273 int Ifpack_NodeFilter::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00274 {
00275   if (!Diagonal.Map().SameAs(*Map_))
00276     IFPACK_CHK_ERR(-1);
00277   Diagonal = *Diagonal_;
00278   return(0);
00279 }
00280 
00281 //==============================================================================
00282 /*
00283 //old Apply (no communication)
00284 int Ifpack_NodeFilter::Apply(const Epetra_MultiVector& X,
00285       Epetra_MultiVector& Y) const 
00286 {
00287 
00288   // skip expensive checks, I suppose input data are ok
00289 
00290   Y.PutScalar(0.0);
00291   int NumVectors = Y.NumVectors();
00292 
00293   double** X_ptr;
00294   double** Y_ptr;
00295   X.ExtractView(&X_ptr);
00296   Y.ExtractView(&Y_ptr);
00297 
00298   for (int i = 0 ; i < NumRows_ ; ++i) {
00299     
00300     int Nnz;
00301     int ierr = Matrix_->ExtractMyRowCopy(i,MaxNumEntriesA_,Nnz,&Values_[0],
00302                                          &Indices_[0]);
00303     IFPACK_CHK_ERR(ierr);
00304 
00305     for (int j = 0 ; j < Nnz ; ++j) {
00306       if (Indices_[j] < NumRows_ ) {
00307     for (int k = 0 ; k < NumVectors ; ++k)
00308       Y_ptr[k][i] += Values_[j] * X_ptr[k][Indices_[j]];
00309       }
00310     }
00311   }
00312 
00313   return(0);
00314 }
00315 */
00316 int Ifpack_NodeFilter::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00317   //
00318   // This function forms the product Y = A * X.
00319   //
00320 
00321   int NumEntries;
00322 
00323   int NumVectors = X.NumVectors();
00324   if (NumVectors!=Y.NumVectors()) {
00325     EPETRA_CHK_ERR(-1); // Need same number of vectors in each MV
00326   }
00327 
00328   // Make sure Import and Export Vectors are compatible
00329   UpdateImportVector(NumVectors);
00330   UpdateExportVector(NumVectors);
00331 
00332   double ** Xp = (double**) X.Pointers();
00333   double ** Yp = (double**) Y.Pointers();
00334 
00335 
00336   // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00337   if (Importer()!=0) {
00338     EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer(), Insert));
00339     Xp = (double**)ImportVector_->Pointers();
00340   }
00341 
00342   // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00343   if (Exporter()!=0) {
00344     Yp = (double**)ExportVector_->Pointers();
00345   }
00346 
00347   // Do actual computation
00348   assert(ovA_!=0);
00349   int *MyRows;
00350   double *MyValues;
00351   int *MyIndices;
00352 
00353   if(Acrs_){
00354     // A rows - CrsMatrix Case
00355     IFPACK_CHK_ERR(Acrs_->ExtractCrsDataPointers(MyRows,MyIndices,MyValues));
00356     //special case NumVectors==1
00357     if (NumVectors==1) {
00358 #ifdef OLD_AND_BUSTED
00359 
00360       for(int i=0;i<NumMyRowsA_;i++) {
00361         int LocRow=Ar_LIDMap_[i];
00362         double sum = 0.0;
00363         for(int j = MyRows[i]; j < MyRows[i+1]; j++)
00364           sum += MyValues[j]*Xp[0][Ac_LIDMap_[MyIndices[j]]];          
00365         Yp[0][LocRow] = sum;
00366       }
00367 #else      
00368       int izero=0;
00369       for(int i=0;i<NumMyColsA_;i++) tempX_[i]=Xp[0][Ac_LIDMap_[i]];
00370       EPETRA_DCRSMV_F77(&izero,&NumMyRowsA_,&NumMyRowsA_,MyValues,MyIndices,MyRows,tempX_,tempY_);
00371 
00372       /*      for(int i=0;i<NumMyRowsA_;i++) {
00373         double sum = 0.0;
00374         for(int j = MyRows[i]; j < MyRows[i+1]; j++)
00375           sum += MyValues[j]*tempX_[MyIndices[j]];
00376         tempY_[i] = sum;
00377         }*/
00378       
00379       for(int i=0;i<NumMyRowsA_;i++) Yp[0][Ar_LIDMap_[i]]=tempY_[i];    
00380 #endif
00381 
00382     }
00383     else {
00384       for(int i=0;i<NumMyRowsA_;i++) {
00385         int LocRow=Ar_LIDMap_[i];
00386         for (int k=0; k<NumVectors; k++) {
00387           double sum = 0.0;
00388           for(int j = MyRows[i]; j < MyRows[i+1]; j++)
00389             sum += MyValues[j]*Xp[k][Ac_LIDMap_[MyIndices[j]]];          
00390           Yp[k][LocRow] = sum;
00391         }
00392       }
00393     }
00394   }
00395   else{
00396     // A rows - RowMatrix Case
00397     MyValues=&Values_[0];
00398     MyIndices=&Indices_[0];
00399     for(int i=0;i<NumMyRowsA_;i++) {
00400       ovA_->A().ExtractMyRowCopy(i,MaxNumEntries_,NumEntries,MyValues,MyIndices);
00401       int LocRow=Ar_LIDMap_[i];
00402       for (int k=0; k<NumVectors; k++) {
00403         double sum = 0.0;
00404         for(int j = 0; j < NumEntries; j++)
00405         sum += MyValues[j]*Xp[k][Ac_LIDMap_[MyIndices[j]]];          
00406         Yp[k][LocRow] = sum;
00407       }
00408     }
00409   }
00410 
00411   // B rows, always CrsMatrix
00412   IFPACK_CHK_ERR(ovA_->B().ExtractCrsDataPointers(MyRows,MyIndices,MyValues));
00413   //special case NumVectors==1
00414   if (NumVectors==1) {
00415     for(int i=0;i<NumMyRowsB_;i++) {
00416       int LocRow=Br_LIDMap_[i];
00417       double sum = 0.0;
00418       for(int j = MyRows[i]; j < MyRows[i+1]; j++)
00419         sum += MyValues[j]*Xp[0][Bc_LIDMap_[MyIndices[j]]];          
00420       Yp[0][LocRow] = sum;
00421     }
00422   } else {
00423     for(int i=0;i<NumMyRowsB_;i++) {
00424       int LocRow=Br_LIDMap_[i];
00425       for (int k=0; k<NumVectors; k++) { //FIXME optimization, check for NumVectors=1
00426         double sum = 0.0;
00427         for(int j = MyRows[i]; j < MyRows[i+1]; j++)
00428           sum += MyValues[j]*Xp[k][Bc_LIDMap_[MyIndices[j]]];          
00429         Yp[k][LocRow] = sum;
00430       }
00431     }
00432   }
00433 
00434   if (Exporter()!=0) {
00435     Y.PutScalar(0.0);  // Make sure target is zero
00436     Y.Export(*ExportVector_, *Exporter(), Add); // Fill Y with Values from export vector
00437   }
00438   // Handle case of rangemap being a local replicated map
00439   if (!OperatorRangeMap().DistributedGlobal() && Comm().NumProc()>1) EPETRA_CHK_ERR(Y.Reduce());
00440 
00441   return(0);
00442 } //Apply
00443 
00444 //==============================================================================
00445 
00446 void Ifpack_NodeFilter::UpdateImportVector(int NumVectors) const {
00447   if(Importer() != 0) {
00448     if(ImportVector_ != 0) {
00449       if(ImportVector_->NumVectors() != NumVectors) {
00450      delete ImportVector_;
00451      ImportVector_= 0;
00452       }
00453     }
00454     if(ImportVector_ == 0)
00455       ImportVector_ = new Epetra_MultiVector(Importer_->TargetMap(),NumVectors); // Create Import vector if needed
00456   }
00457   return;
00458 /*
00459   if(ImportVector_ == null || ImportVector_->NumVectors() != NumVectors)
00460     ImportVector_ = rcp(new Epetra_MultiVector(Importer_->TargetMap(),NumVectors));
00461 */
00462 }
00463 
00464 //=======================================================================================================
00465 
00466 void Ifpack_NodeFilter::UpdateExportVector(int NumVectors) const {
00467   if(Exporter() != 0) {
00468     if(ExportVector_ != 0) {
00469       if(ExportVector_->NumVectors() != NumVectors) {
00470      delete ExportVector_;
00471      ExportVector_= 0;
00472       }
00473     }
00474     if(ExportVector_ == 0)
00475       ExportVector_ = new Epetra_MultiVector(Exporter_->SourceMap(),NumVectors); // Create Export vector if needed
00476   }
00477   return;
00478 /*
00479   if(ExportVector_ == null || ExportVector_->NumVectors() != NumVectors)
00480     ExportVector_ = rcp(new Epetra_MultiVector(Exporter_->SourceMap(),NumVectors));
00481 */
00482 }
00483 
00484 //=======================================================================================================
00485 int Ifpack_NodeFilter::ApplyInverse(const Epetra_MultiVector& X,
00486          Epetra_MultiVector& Y) const
00487 {
00488   IFPACK_CHK_ERR(-1); // not implemented
00489 }
00490 
00491 //==============================================================================
00492 const Epetra_BlockMap& Ifpack_NodeFilter::Map() const
00493 {
00494   return(*Map_);
00495 
00496 }
00497 
00498 #endif //ifdef IFPACK_NODE_AWARE_CODE

Generated on Wed May 12 21:30:18 2010 for IFPACK by  doxygen 1.4.7