Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_SubdomainFilter.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //              Ifpack_SubdomainFilter
00005 //              Author: Radu Popescu <radu.popescu@epfl.ch>
00006 //              Copyright 2011 EPFL - CMCS
00007 // 
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are
00010 // met:
00011 //
00012 // 1. Redistributions of source code must retain the above copyright
00013 // notice, this list of conditions and the following disclaimer.
00014 //
00015 // 2. Redistributions in binary form must reproduce the above copyright
00016 // notice, this list of conditions and the following disclaimer in the
00017 // documentation and/or other materials provided with the distribution.
00018 //
00019 // 3. Neither the name of the copyright holder nor the names of the
00020 // contributors may be used to endorse or promote products derived from
00021 // this software without specific prior written permission.
00022 //
00023 // THIS SOFTWARE IS PROVIDED BY EPFL - CMCS "AS IS" AND ANY EXPRESS OR
00024 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
00025 // WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00026 // ARE DISCLAIMED. IN NO EVENT SHALL EPFL - CMCS OR THE CONTRIBUTORS
00027 // BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
00028 // OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
00029 // OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
00030 // BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
00031 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
00032 // OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
00033 // EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00034 //
00035 // ************************************************************************
00036 //@HEADER
00037 
00038 #ifdef IFPACK_SUBCOMM_CODE
00039 
00040 #include <vector>
00041 #include <algorithm>
00042 
00043 #include "Ifpack_ConfigDefs.h"
00044 
00045 #include "Epetra_MultiVector.h"
00046 #include "Epetra_Vector.h"
00047 #include "Epetra_IntVector.h"
00048 #include "Epetra_RowMatrix.h"
00049 #include "Epetra_Map.h"
00050 #include "Epetra_BlockMap.h"
00051 #include "Ifpack_SubdomainFilter.h"
00052 #include "Ifpack_OverlappingRowMatrix.h"
00053 #include "Epetra_Import.h"
00054 #include "Epetra_Export.h"
00055 #include "Epetra_CrsMatrix.h"
00056 #include "Epetra_BLAS_wrappers.h"
00057 
00058 using namespace Teuchos;
00059 
00060 //==============================================================================
00061 Ifpack_SubdomainFilter::Ifpack_SubdomainFilter(const RefCountPtr<const Epetra_RowMatrix>& Matrix,int subdomainId) :
00062   Matrix_(Matrix),
00063   NumMyRows_(0),
00064   NumMyNonzeros_(0),
00065   NumGlobalRows_(0),
00066   NumGlobalCols_(0),
00067   MaxNumEntries_(0),
00068   MaxNumEntriesA_(0)
00069 {
00070   sprintf(Label_,"%s","Ifpack_SubdomainFilter");
00071 
00072   ImportVector_ = 0;
00073   ExportVector_ = 0;
00074 
00075   ovA_ = dynamic_cast<const Ifpack_OverlappingRowMatrix*>(&*Matrix_);
00076 
00077   if (ovA_) {
00078     Acrs_=dynamic_cast<const Epetra_CrsMatrix*>(&ovA_->A());
00079     NumMyRowsA_ = ovA_->A().NumMyRows();
00080     NumMyRowsB_ = ovA_->B().NumMyRows();
00081   } else {
00082     NumMyRowsA_ = Matrix->NumMyRows();
00083     NumMyRowsB_ = 0;
00084   }
00085   
00086 #ifdef HAVE_MPI
00087   const Epetra_MpiComm *pComm = dynamic_cast<const Epetra_MpiComm*>( &(Matrix->Comm()) );
00088   assert(pComm != NULL);
00089   MPI_Comm_split(pComm->Comm(), subdomainId, pComm->MyPID(), &subdomainMPIComm_);
00090   SubComm_ = rcp( new Epetra_MpiComm(subdomainMPIComm_) );
00091 #else
00092   SubComm_ = rcp( new Epetra_SerialComm );
00093 #endif
00094 
00095   NumMyRows_ = Matrix->NumMyRows();
00096   SubComm_->SumAll(&NumMyRows_,&NumGlobalRows_,1);
00097 
00098   // Get the row GID corresponding to the process
00099   const Epetra_Map &globRowMap = Matrix->RowMatrixRowMap();
00100   std::vector<int> myRowsGID(NumMyRows_);
00101   globRowMap.MyGlobalElements(&myRowsGID[0]);
00102 
00103   // Gather the GID of all rows in the subdomain
00104   // Get the maximum number of local rows on each processor in the subdomain
00105   // Allocate a vector large enough (Proc cout * max number local rows
00106   // Gather all the local row indices. If a process has less than max num
00107   // local rows, pad the local vector with -1;
00108   // After the gather, eliminate the -1 elements from the target vector
00109   //
00110   // TODO: this section should be rewritten to avoid the GatherAll call.
00111   //       Right now it's not memory efficient. It could be moved to
00112   //       Epetra_Map since it a map operation
00113 
00114   int MaxLocalRows = 0;
00115   SubComm_->MaxAll(&NumMyRows_, &MaxLocalRows, 1);
00116   
00117   std::vector<int> SubdomainRowGID(SubComm_->NumProc() * MaxLocalRows);
00118   myRowsGID.resize(MaxLocalRows, -1);
00119 
00120   SubComm_->GatherAll(&myRowsGID[0],
00121                       &SubdomainRowGID[0],
00122                       MaxLocalRows);
00123 
00124   std::sort(SubdomainRowGID.begin(), SubdomainRowGID.end());
00125 
00126   int PaddingSize = SubdomainRowGID.size() - NumGlobalRows_;
00127   SubdomainRowGID.erase(SubdomainRowGID.begin(),
00128                          SubdomainRowGID.begin() + PaddingSize);
00129 
00130   // Get the col GID corresponding to the process
00131   const Epetra_Map &globColMap = Matrix->RowMatrixColMap();
00132   NumMyCols_ = globColMap.NumMyElements();
00133   std::vector<int> myGlobalCols(NumMyCols_);
00134   globColMap.MyGlobalElements(&myGlobalCols[0]);
00135 
00136   // Eliminate column GID that are outside the subdomain
00137   std::vector<int> filteredColGID;
00138   for (int i = 0; i < NumMyCols_; ++i) {
00139     if (std::find(SubdomainRowGID.begin(), SubdomainRowGID.end(),
00140       myGlobalCols[i]) != SubdomainRowGID.end()) {
00141       filteredColGID.push_back(myGlobalCols[i]);
00142     }
00143   }
00144   NumMyCols_ = filteredColGID.size();
00145 
00146   // Create the maps with the reindexed GIDs
00147   Map_ = rcp( new Epetra_Map(-1, NumMyRows_, &myRowsGID[0], globRowMap.IndexBase(), *SubComm_) );
00148   colMap_ = rcp( new Epetra_Map(-1, NumMyCols_, &filteredColGID[0], globColMap.IndexBase(), *SubComm_) );
00149   NumGlobalCols_ = NumGlobalRows_;
00150 
00151   NumEntries_.resize(NumMyRows_);
00152   MaxNumEntriesA_ = Matrix->MaxNumEntries();
00153   MaxNumEntries_ = Matrix->MaxNumEntries();
00154 
00155   Diagonal_ = rcp( new Epetra_Vector(*Map_) );
00156   if (Diagonal_ == Teuchos::null) IFPACK_CHK_ERRV(-5);
00157 
00158   Indices_.resize(MaxNumEntries_);
00159   Values_.resize(MaxNumEntries_);
00160 
00161   Ac_LIDMap_ = 0;
00162   Bc_LIDMap_ = 0;
00163   Ar_LIDMap_ = 0;
00164   Br_LIDMap_ = 0;
00165 
00166   if(ovA_){
00167     Ac_LIDMap_ = new int[ovA_->A().NumMyCols() + 1];
00168     for(int i = 0; i < ovA_->A().NumMyCols(); i++) {
00169       Ac_LIDMap_[i] = colMap_->LID(ovA_->A().RowMatrixColMap().GID(i));    
00170     }
00171     Bc_LIDMap_ = new int[ovA_->B().NumMyCols() + 1];
00172     for(int i = 0; i < ovA_->B().NumMyCols(); i++) {
00173       Bc_LIDMap_[i] = colMap_->LID(ovA_->B().RowMatrixColMap().GID(i));
00174     }
00175     Ar_LIDMap_ = new int[ovA_->A().NumMyRows() + 1];
00176     for(int i = 0; i < ovA_->A().NumMyRows(); i++) {
00177       Ar_LIDMap_[i] = Map_->LID(ovA_->A().RowMatrixRowMap().GID(i));    
00178     }
00179     Br_LIDMap_ = new int[ovA_->B().NumMyRows() + 1];
00180     for(int i = 0; i < ovA_->B().NumMyRows(); i++) {
00181       Br_LIDMap_[i] = Map_->LID(ovA_->B().RowMatrixRowMap().GID(i));
00182     }
00183   }
00184 
00185   int ActualMaxNumEntries = 0;
00186 
00187   for (int i = 0 ; i < NumMyRows_; ++i) {
00188     NumEntries_[i] = 0;
00189     int Nnz;
00190     IFPACK_CHK_ERRV(ExtractMyRowCopy(i, MaxNumEntries_, Nnz, &Values_[0], &Indices_[0]));
00191 
00192     for (int j = 0 ; j < Nnz; ++j) {
00193       if (Indices_[j] == i) {
00194         (*Diagonal_)[i] = Values_[j];
00195       }
00196     }
00197 
00198     if (Nnz > ActualMaxNumEntries) {
00199       ActualMaxNumEntries = Nnz;
00200     }
00201 
00202     NumMyNonzeros_ += Nnz;
00203     NumEntries_[i] = Nnz;
00204   }
00205 
00206   SubComm_->SumAll(&NumMyNonzeros_, &NumGlobalNonzeros_, 1);
00207   MaxNumEntries_ = ActualMaxNumEntries;
00208 
00209   int gpid = Matrix->Comm().MyPID();
00210   Exporter_ = null;
00211   Importer_ = null;
00212 
00213   if (!(RowMatrixRowMap().SameAs(OperatorRangeMap()))) {
00214     try{
00215       Exporter_ = rcp(new Epetra_Export(RowMatrixRowMap(), OperatorRangeMap()));
00216     }
00217     catch(...) {
00218       printf("** * gpid %d: Ifpack_SubdomainFilter ctor: problem creating Exporter_ * **\n\n",gpid);
00219     }
00220   }
00221 
00222   if (!(*colMap_).SameAs(*Map_)) {
00223     try{
00224       Importer_ = rcp(new Epetra_Import(*colMap_, *Map_));
00225     }
00226     catch(...) {
00227       printf("** * gpid %d: Ifpack_SubdomainFilter ctor: problem creating Importer_ * **\n\n",gpid);
00228     }
00229   }
00230 
00231 } //Ifpack_SubdomainFilter() ctor
00232 
00233 //==============================================================================
00234 Ifpack_SubdomainFilter::~Ifpack_SubdomainFilter()
00235 {
00236     if(Ac_LIDMap_) delete [] Ac_LIDMap_;
00237     if(Bc_LIDMap_) delete [] Bc_LIDMap_;
00238     if(Ar_LIDMap_) delete [] Ar_LIDMap_;
00239     if(Br_LIDMap_) delete [] Br_LIDMap_;
00240 
00241     if(ImportVector_) delete ImportVector_; 
00242 } //Ifpack_SubdomainFilter destructor
00243 
00244 //==============================================================================
00245 int Ifpack_SubdomainFilter:: ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, 
00246     double *Values, int * Indices) const
00247 {
00248   if ((MyRow < 0) || (MyRow >= NumMyRows_)) {
00249     IFPACK_CHK_ERR(-1);
00250   }
00251 
00252   assert(Length >= NumEntries_[MyRow]);
00253 
00254   int ierr;
00255   if (ovA_) {
00256     int LocRow = ovA_->A().RowMatrixRowMap().LID(Map_->GID(MyRow));      
00257     if (LocRow != -1) { 
00258       ierr = ovA_->A().ExtractMyRowCopy(LocRow, Length, NumEntries, Values, Indices);
00259       for(int j = 0;j < NumEntries; j++) {
00260         Indices[j] = Ac_LIDMap_[Indices[j]];
00261       }
00262     }
00263     else {
00264       LocRow = ovA_->B().RowMatrixRowMap().LID(Map_->GID(MyRow));      
00265       ierr = ovA_->B().ExtractMyRowCopy(LocRow, Length, NumEntries, Values, Indices);
00266       for(int j = 0; j < NumEntries; j++) {
00267         Indices[j] = Bc_LIDMap_[Indices[j]];
00268       }
00269     }
00270   } else {
00271     int Nnz = 0;
00272     ierr = Matrix_->ExtractMyRowCopy(MyRow, MaxNumEntriesA_, Nnz, &Values_[0], &Indices_[0]);
00273     IFPACK_CHK_ERR(ierr);
00274 
00275     NumEntries = 0;
00276     for (int j = 0 ; j < Nnz ; ++j) {
00277       int subdomainLocalIndex = colMap_->LID(Matrix_->RowMatrixColMap().GID(Indices_[j]));
00278       if (subdomainLocalIndex != -1) {
00279         Indices[NumEntries] = subdomainLocalIndex;
00280         Values[NumEntries] = Values_[j];
00281         NumEntries++;
00282       }
00283     }
00284   }
00285 
00286   return(ierr);
00287 }
00288 
00289 //==============================================================================
00290 int Ifpack_SubdomainFilter::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00291 {
00292   if (!Diagonal.Map().SameAs(*Map_)) {
00293     IFPACK_CHK_ERR(-1);
00294   }
00295 
00296   Diagonal = *Diagonal_;
00297   return(0);
00298 }
00299 
00300 //==============================================================================
00301 int Ifpack_SubdomainFilter::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00302 {
00303     std::cout << "Ifpack_SubdomainFilter::Apply not implemented.\n";
00304     IFPACK_CHK_ERR(-1); // not implemented
00305 }
00306 
00307 //==============================================================================
00308 void Ifpack_SubdomainFilter::UpdateImportVector(int NumVectors) const
00309 {
00310 }
00311 
00312 //==============================================================================
00313 void Ifpack_SubdomainFilter::UpdateExportVector(int NumVectors) const
00314 {
00315 }
00316 
00317 //==============================================================================
00318 int Ifpack_SubdomainFilter::ApplyInverse(const Epetra_MultiVector& X,
00319      Epetra_MultiVector& Y) const
00320 {
00321   std::cout << "Ifpack_SubdomainFilter::ApplyInverse not implemented.\n";
00322   IFPACK_CHK_ERR(-1); // not implemented
00323 }
00324 
00325 //==============================================================================
00326 const Epetra_BlockMap& Ifpack_SubdomainFilter::Map() const
00327 {
00328   return(*Map_);
00329 }
00330 
00331 #endif //ifdef IFPACK_SUBCOMM_CODE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines