IFPACK Development
Ifpack_AdditiveSchwarz.h
00001 #ifndef IFPACK_ADDITIVESCHWARZ_H
00002 #define IFPACK_ADDITIVESCHWARZ_H
00003 
00004 #include "Ifpack_ConfigDefs.h"
00005 #include "Ifpack_Preconditioner.h"
00006 #include "Ifpack_ConfigDefs.h"
00007 #include "Ifpack_Preconditioner.h"
00008 #include "Ifpack_Reordering.h"
00009 #include "Ifpack_RCMReordering.h"
00010 #include "Ifpack_METISReordering.h"
00011 #include "Ifpack_LocalFilter.h"
00012 #include "Ifpack_NodeFilter.h"
00013 #ifdef IFPACK_SUBCOMM_CODE
00014 #include "Ifpack_SubdomainFilter.h"
00015 #endif
00016 #include "Ifpack_SingletonFilter.h"
00017 #include "Ifpack_ReorderFilter.h"
00018 #include "Ifpack_Utils.h"
00019 #include "Ifpack_OverlappingRowMatrix.h"
00020 #include "Epetra_CombineMode.h"
00021 #include "Epetra_MultiVector.h"
00022 #include "Epetra_Map.h"
00023 #include "Epetra_Comm.h"
00024 #include "Epetra_Time.h"
00025 #include "Epetra_LinearProblem.h"
00026 #include "Epetra_RowMatrix.h"
00027 #include "Epetra_CrsMatrix.h"
00028 #include "Teuchos_ParameterList.hpp"
00029 #include "Teuchos_RefCountPtr.hpp"
00030 
00031 #ifdef IFPACK_SUBCOMM_CODE
00032 #include "EpetraExt_Reindex_CrsMatrix.h"
00033 #include "EpetraExt_Reindex_MultiVector.h"
00034 #endif
00035 #ifdef IFPACK_NODE_AWARE_CODE
00036 #include "EpetraExt_OperatorOut.h"
00037 #include "EpetraExt_RowMatrixOut.h"
00038 #include "EpetraExt_BlockMapOut.h"
00039 #endif
00040 
00041 #ifdef HAVE_IFPACK_AMESOS
00042   #include "Ifpack_AMDReordering.h"
00043 #endif
00044 
00045 
00047 
00100 template<typename T>
00101 class Ifpack_AdditiveSchwarz : public virtual Ifpack_Preconditioner {
00102       
00103 public:
00104 
00106 
00107 
00118   Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
00119              int OverlapLevel_in = 0);
00120   
00122   virtual ~Ifpack_AdditiveSchwarz() {};
00124 
00126 
00128 
00137     virtual int SetUseTranspose(bool UseTranspose_in);
00139   
00141 
00143 
00153     virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00154 
00156 
00167     virtual int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00168 
00170     virtual double NormInf() const;
00172   
00174 
00176     virtual const char * Label() const;
00177 
00179     virtual bool UseTranspose() const;
00180 
00182     virtual bool HasNormInf() const;
00183 
00185     virtual const Epetra_Comm & Comm() const;
00186 
00188     virtual const Epetra_Map & OperatorDomainMap() const;
00189 
00191     virtual const Epetra_Map & OperatorRangeMap() const;
00193 
00195   virtual bool IsInitialized() const
00196   {
00197     return(IsInitialized_);
00198   }
00199 
00201   virtual bool IsComputed() const
00202   {
00203     return(IsComputed_);
00204   }
00205 
00207 
00226   virtual int SetParameters(Teuchos::ParameterList& List);
00227 
00228   // @}
00229 
00230   // @{ Query methods
00231   
00233   virtual int Initialize();
00234 
00236   virtual int Compute();
00237 
00239   virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
00240                          const int MaxIters = 1550,
00241                          const double Tol = 1e-9,
00242              Epetra_RowMatrix* Matrix_in = 0);
00243 
00245   virtual double Condest() const
00246   {
00247     return(Condest_);
00248   }
00249 
00251   virtual const Epetra_RowMatrix& Matrix() const
00252   {
00253     return(*Matrix_);
00254   }
00255 
00257   virtual bool IsOverlapping() const
00258   {
00259     return(IsOverlapping_);
00260   }
00261 
00263   virtual std::ostream& Print(std::ostream&) const;
00264   
00265   virtual const T* Inverse() const
00266   {
00267     return(&*Inverse_);
00268   }
00269 
00271   virtual int NumInitialize() const
00272   {
00273     return(NumInitialize_);
00274   }
00275 
00277   virtual int NumCompute() const
00278   {
00279     return(NumCompute_);
00280   }
00281 
00283   virtual int NumApplyInverse() const
00284   {
00285     return(NumApplyInverse_);
00286   }
00287 
00289   virtual double InitializeTime() const
00290   {
00291     return(InitializeTime_);
00292   }
00293 
00295   virtual double ComputeTime() const
00296   {
00297     return(ComputeTime_);
00298   }
00299 
00301   virtual double ApplyInverseTime() const
00302   {
00303     return(ApplyInverseTime_);
00304   }
00305 
00307   virtual double InitializeFlops() const
00308   {
00309     return(InitializeFlops_);
00310   }
00311 
00312   virtual double ComputeFlops() const
00313   {
00314     return(ComputeFlops_);
00315   }
00316 
00317   virtual double ApplyInverseFlops() const
00318   {
00319     return(ApplyInverseFlops_);
00320   }
00321 
00323   virtual int OverlapLevel() const
00324   {
00325     return(OverlapLevel_);
00326   }
00327 
00329   virtual const Teuchos::ParameterList& List() const
00330   {
00331     return(List_);
00332   }
00333 
00334 protected:
00335 
00336   // @}
00337 
00338   // @{ Internal merhods.
00339   
00341   Ifpack_AdditiveSchwarz(const Ifpack_AdditiveSchwarz& RHS)
00342   { }
00343 
00345   int Setup();
00346   
00347   // @}
00348 
00349   // @{ Internal data.
00350   
00352   Teuchos::RefCountPtr<const Epetra_RowMatrix> Matrix_;
00354   Teuchos::RefCountPtr<Ifpack_OverlappingRowMatrix> OverlappingMatrix_;
00356 /*
00357   //TODO if we choose to expose the node aware code, i.e., no ifdefs,
00358   //TODO then we should switch to this definition.
00359   Teuchos::RefCountPtr<Epetra_RowMatrix> LocalizedMatrix_;
00360 */
00361 #ifdef IFPACK_SUBCOMM_CODE
00362 
00363   Teuchos::RCP<Epetra_RowMatrix> LocalizedMatrix_;
00365   Teuchos::RCP<Epetra_CrsMatrix> SubdomainCrsMatrix_;
00367   Teuchos::RCP<Epetra_CrsMatrix> ReindexedCrsMatrix_;
00368 
00369   // The following data members are needed when doing ApplyInverse
00370   // with an Ifpack_SubdomainFilter local matrix
00371   Teuchos::RCP<Epetra_Map> tempMap_;
00372   Teuchos::RCP<Epetra_Map> tempDomainMap_;
00373   Teuchos::RCP<Epetra_Map> tempRangeMap_;
00374   Teuchos::RCP<EpetraExt::CrsMatrix_Reindex> SubdomainMatrixReindexer_;
00375   Teuchos::RCP<EpetraExt::MultiVector_Reindex> DomainVectorReindexer_;
00376   Teuchos::RCP<EpetraExt::MultiVector_Reindex> RangeVectorReindexer_;
00377   mutable Teuchos::RCP<Epetra_MultiVector> tempX_;
00378   mutable Teuchos::RCP<Epetra_MultiVector> tempY_;
00379 #else
00380 # ifdef IFPACK_NODE_AWARE_CODE
00381   Teuchos::RefCountPtr<Ifpack_NodeFilter> LocalizedMatrix_;
00382 # else
00383   Teuchos::RefCountPtr<Ifpack_LocalFilter> LocalizedMatrix_;
00384 # endif
00385 #endif
00386 #ifdef IFPACK_SUBCOMM_CODE
00387   // Some data members used for determining the subdomain id (color)
00389   int MpiRank_;
00391   int NumMpiProcs_;
00393   int NumMpiProcsPerSubdomain_;
00395   int NumSubdomains_;
00397   int SubdomainId_;
00398 #endif
00399 
00400   string Label_;
00402   bool IsInitialized_;
00404   bool IsComputed_;
00406   bool UseTranspose_;
00408   bool IsOverlapping_;
00410   int OverlapLevel_;
00412   Teuchos::ParameterList List_;
00414   Epetra_CombineMode CombineMode_;
00416   double Condest_;
00418   bool ComputeCondest_;
00420   bool UseReordering_;
00422   string ReorderingType_;
00424   Teuchos::RefCountPtr<Ifpack_Reordering> Reordering_;
00426   Teuchos::RefCountPtr<Ifpack_ReorderFilter> ReorderedLocalizedMatrix_;
00428   bool FilterSingletons_;
00430   Teuchos::RefCountPtr<Ifpack_SingletonFilter> SingletonFilter_;
00432   int NumInitialize_;
00434   int NumCompute_;
00436   mutable int NumApplyInverse_;
00438   double InitializeTime_;
00440   double ComputeTime_;
00442   mutable double ApplyInverseTime_;
00444   double InitializeFlops_;
00446   double ComputeFlops_;
00448   mutable double ApplyInverseFlops_;
00450   Teuchos::RefCountPtr<Epetra_Time> Time_;
00452   Teuchos::RefCountPtr<T> Inverse_;
00454 #ifdef IFPACK_SUBCOMM_CODE
00455   mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
00456   mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
00457 #endif
00458 # ifdef IFPACK_NODE_AWARE_CODE
00459   mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
00460   mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
00461 #endif
00462 
00463 }; // class Ifpack_AdditiveSchwarz<T>
00464 
00465 //==============================================================================
00466 template<typename T>
00467 Ifpack_AdditiveSchwarz<T>::
00468 Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
00469                        int OverlapLevel_in) :
00470 #ifdef IFPACK_SUBCOMM_CODE
00471   MpiRank_(0),
00472   NumMpiProcs_(1),
00473   NumMpiProcsPerSubdomain_(1),
00474   NumSubdomains_(1),
00475   SubdomainId_(0),
00476 #endif
00477   IsInitialized_(false),
00478   IsComputed_(false),
00479   UseTranspose_(false),
00480   IsOverlapping_(false),
00481   OverlapLevel_(OverlapLevel_in),
00482   CombineMode_(Zero),
00483   Condest_(-1.0),
00484   ComputeCondest_(true),
00485   UseReordering_(false),
00486   ReorderingType_("none"),
00487   FilterSingletons_(false),
00488   NumInitialize_(0),
00489   NumCompute_(0),
00490   NumApplyInverse_(0),
00491   InitializeTime_(0.0),
00492   ComputeTime_(0.0),
00493   ApplyInverseTime_(0.0),
00494   InitializeFlops_(0.0),
00495   ComputeFlops_(0.0),
00496   ApplyInverseFlops_(0.0)
00497 {
00498   // Construct a reference-counted pointer with the input matrix, don't manage the memory.
00499   Matrix_ = Teuchos::rcp( Matrix_in, false );
00500 
00501   if (Matrix_->Comm().NumProc() == 1)
00502     OverlapLevel_ = 0;
00503 
00504   if ((OverlapLevel_ != 0) && (Matrix_->Comm().NumProc() > 1))
00505     IsOverlapping_ = true;
00506   // Sets parameters to default values
00507   Teuchos::ParameterList List_in;
00508   SetParameters(List_in);
00509 }
00510 
00511 #ifdef IFPACK_NODE_AWARE_CODE
00512 extern int ML_NODE_ID;
00513 #endif
00514 
00515 //==============================================================================
00516 template<typename T>
00517 int Ifpack_AdditiveSchwarz<T>::Setup()
00518 {
00519 
00520   Epetra_RowMatrix* MatrixPtr;
00521 
00522 #ifndef IFPACK_SUBCOMM_CODE
00523 # ifdef IFPACK_NODE_AWARE_CODE
00524 /*
00525   sleep(3);
00526   if (Comm().MyPID() == 0) cout << "Printing out ovArowmap" << endl;
00527   Comm().Barrier();
00528 
00529   EpetraExt::BlockMapToMatrixMarketFile("ovArowmap",OverlappingMatrix_->RowMatrixRowMap());
00530   if (Comm().MyPID() == 0) cout << "Printing out ovAcolmap" << endl;
00531   Comm().Barrier();
00532   EpetraExt::BlockMapToMatrixMarketFile("ovAcolmap",OverlappingMatrix_->RowMatrixColMap());
00533   Comm().Barrier();
00534 */
00535 /*
00536   EpetraExt::RowMatrixToMatlabFile("ovA",*OverlappingMatrix_);
00537   fprintf(stderr,"p %d n %d matrix file done\n",Comm().MyPID(),ML_NODE_ID);
00538   Comm().Barrier();
00539 */
00540   int nodeID;
00541   try{ nodeID = List_.get("ML node id",0);}
00542   catch(...){fprintf(stderr,"%s","Ifpack_AdditiveSchwarz<T>::Setup(): no parameter \"ML node id\"\n\n");
00543              cout << List_ << endl;}
00544 # endif
00545 #endif
00546 
00547   try{
00548   if (OverlappingMatrix_ != Teuchos::null)
00549   {
00550 #ifdef IFPACK_SUBCOMM_CODE
00551     if (NumMpiProcsPerSubdomain_ > 1) {
00552       LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(OverlappingMatrix_, SubdomainId_));
00553     } else {
00554       LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(OverlappingMatrix_));
00555     }
00556 #else
00557 #   ifdef IFPACK_NODE_AWARE_CODE
00558     Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
00559     LocalizedMatrix_ = Teuchos::rcp(tt);
00560     //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
00561 #   else
00562     LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
00563 #   endif
00564 #endif
00565   }
00566   else
00567   {
00568 #ifdef IFPACK_SUBCOMM_CODE
00569     if (NumMpiProcsPerSubdomain_ > 1) {
00570       LocalizedMatrix_ = Teuchos::rcp(new Ifpack_SubdomainFilter(Matrix_, SubdomainId_));
00571     } else {
00572       LocalizedMatrix_ = Teuchos::rcp(new Ifpack_LocalFilter(Matrix_));
00573     }
00574 #else
00575 #   ifdef IFPACK_NODE_AWARE_CODE
00576     Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(Matrix_,nodeID); //FIXME
00577     LocalizedMatrix_ = Teuchos::rcp(tt);
00578     //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
00579 #   else
00580     LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
00581 #   endif
00582 #endif
00583   }
00584   }
00585   catch(...) {
00586      fprintf(stderr,"%s","AdditiveSchwarz Setup: problem creating local filter matrix.\n");
00587   }
00588 
00589   if (LocalizedMatrix_ == Teuchos::null)
00590     IFPACK_CHK_ERR(-5);
00591 
00592   // users may want to skip singleton check
00593   if (FilterSingletons_) {
00594     SingletonFilter_ = Teuchos::rcp( new Ifpack_SingletonFilter(LocalizedMatrix_) );
00595     MatrixPtr = &*SingletonFilter_;
00596   }
00597   else
00598     MatrixPtr = &*LocalizedMatrix_;
00599 
00600   if (UseReordering_) {
00601 
00602     // create reordering and compute it
00603     if (ReorderingType_ == "rcm")
00604       Reordering_ = Teuchos::rcp( new Ifpack_RCMReordering() );
00605     else if (ReorderingType_ == "metis")
00606       Reordering_ = Teuchos::rcp( new Ifpack_METISReordering() );
00607 #ifdef HAVE_IFPACK_AMESOS   
00608     else if (ReorderingType_ == "amd" )
00609       Reordering_ = Teuchos::rcp( new Ifpack_AMDReordering() );
00610 #endif
00611     else {
00612       cerr << "reordering type not correct (" << ReorderingType_ << ")" << endl;
00613       exit(EXIT_FAILURE);
00614     }
00615     if (Reordering_ == Teuchos::null) IFPACK_CHK_ERR(-5);
00616 
00617     IFPACK_CHK_ERR(Reordering_->SetParameters(List_));
00618     IFPACK_CHK_ERR(Reordering_->Compute(*MatrixPtr));
00619 
00620     // now create reordered localized matrix
00621     ReorderedLocalizedMatrix_ = 
00622       Teuchos::rcp( new Ifpack_ReorderFilter(Teuchos::rcp( MatrixPtr, false ), Reordering_) );
00623 
00624     if (ReorderedLocalizedMatrix_ == Teuchos::null) IFPACK_CHK_ERR(-5);
00625 
00626     MatrixPtr = &*ReorderedLocalizedMatrix_;
00627   }
00628 
00629 #ifdef IFPACK_SUBCOMM_CODE
00630   // The subdomain matrix needs to be reindexed by Amesos so we need to make a CrsMatrix
00631   // and then reindex it with EpetraExt.
00632   // The reindexing is done here because this feature is only implemented in Amesos_Klu,
00633   // and it is desired to have reindexing with other direct solvers in the Amesos package
00634   SubdomainCrsMatrix_.reset(new Epetra_CrsMatrix(Copy, MatrixPtr->RowMatrixRowMap(), -1));
00635   Teuchos::RCP<Epetra_Import> tempImporter = Teuchos::rcp(new Epetra_Import(SubdomainCrsMatrix_->Map(), MatrixPtr->Map()));
00636   SubdomainCrsMatrix_->Import(*MatrixPtr, *tempImporter, Insert);
00637   SubdomainCrsMatrix_->FillComplete();
00638 
00639   tempMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->RowMap().NumGlobalElements(),
00640                                SubdomainCrsMatrix_->RowMap().NumMyElements(),
00641                                0, SubdomainCrsMatrix_->Comm()));
00642   tempRangeMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorRangeMap().NumGlobalElements(),
00643                                     SubdomainCrsMatrix_->OperatorRangeMap().NumMyElements(),
00644                                     0, SubdomainCrsMatrix_->Comm()));
00645   tempDomainMap_.reset(new Epetra_Map(SubdomainCrsMatrix_->OperatorDomainMap().NumGlobalElements(),
00646                                      SubdomainCrsMatrix_->OperatorDomainMap().NumMyElements(),
00647                                      0, SubdomainCrsMatrix_->Comm()));
00648 
00649   SubdomainMatrixReindexer_.reset(new EpetraExt::CrsMatrix_Reindex(*tempMap_));
00650   DomainVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempDomainMap_));
00651   RangeVectorReindexer_.reset(new EpetraExt::MultiVector_Reindex(*tempRangeMap_));
00652 
00653   ReindexedCrsMatrix_.reset(&((*SubdomainMatrixReindexer_)(*SubdomainCrsMatrix_)), false);
00654   
00655   MatrixPtr = &*ReindexedCrsMatrix_;
00656 #endif
00657 
00658   Inverse_ = Teuchos::rcp( new T(MatrixPtr) );
00659 
00660   if (Inverse_ == Teuchos::null)
00661     IFPACK_CHK_ERR(-5);
00662 
00663   return(0);
00664 }
00665 
00666 //==============================================================================
00667 template<typename T>
00668 int Ifpack_AdditiveSchwarz<T>::SetParameters(Teuchos::ParameterList& List_in)
00669 {
00670 #ifdef IFPACK_SUBCOMM_CODE
00671   MpiRank_ = Matrix_->Comm().MyPID();
00672   NumMpiProcs_ = Matrix_->Comm().NumProc();
00673   NumMpiProcsPerSubdomain_ = List_in.get("subdomain: number-of-processors", 1);
00674   NumSubdomains_ = NumMpiProcs_ / NumMpiProcsPerSubdomain_;
00675   SubdomainId_ = MpiRank_ / NumMpiProcsPerSubdomain_;
00676 
00677   if (NumSubdomains_ == 1) {
00678       IsOverlapping_ = false;
00679   }
00680 
00681 #endif
00682   // compute the condition number each time Compute() is invoked.
00683   ComputeCondest_ = List_in.get("schwarz: compute condest", ComputeCondest_);
00684   // combine mode
00685   if( Teuchos::ParameterEntry *combineModeEntry = List_in.getEntryPtr("schwarz: combine mode") )
00686   {
00687     if( typeid(std::string) == combineModeEntry->getAny().type() )
00688     {
00689       std::string mode = List_in.get("schwarz: combine mode", "Add");
00690       if (mode == "Add")
00691         CombineMode_ = Add;
00692       else if (mode == "Zero")
00693         CombineMode_ = Zero;
00694       else if (mode == "Insert")
00695         CombineMode_ = Insert;
00696       else if (mode == "InsertAdd")
00697         CombineMode_ = InsertAdd;
00698       else if (mode == "Average")
00699         CombineMode_ = Average;
00700       else if (mode == "AbsMax")
00701         CombineMode_ = AbsMax;
00702       else
00703       {
00704         TEUCHOS_TEST_FOR_EXCEPTION(
00705           true,std::logic_error
00706           ,"Error, The (Epetra) combine mode of \""<<mode<<"\" is not valid!  Only the values"
00707           " \"Add\", \"Zero\", \"Insert\", \"InsertAdd\", \"Average\", and \"AbsMax\" are accepted!"
00708           );
00709       }
00710     }
00711     else if ( typeid(Epetra_CombineMode) == combineModeEntry->getAny().type() )
00712     {
00713       CombineMode_ = Teuchos::any_cast<Epetra_CombineMode>(combineModeEntry->getAny());
00714     }
00715     else
00716     {
00717       // Throw exception with good error message!
00718       Teuchos::getParameter<std::string>(List_in,"schwarz: combine mode");
00719     }
00720   }
00721   else
00722   {
00723     // Make the default be a string to be consistent with the valid parameters!
00724     List_in.get("schwarz: combine mode","Zero");
00725   }
00726   // type of reordering
00727   ReorderingType_ = List_in.get("schwarz: reordering type", ReorderingType_);
00728   if (ReorderingType_ == "none")
00729     UseReordering_ = false;
00730   else 
00731     UseReordering_ = true;
00732   // if true, filter singletons. NOTE: the filtered matrix can still have
00733   // singletons! A simple example: upper triangular matrix, if I remove
00734   // the lower node, I still get a matrix with a singleton! However, filter
00735   // singletons should help for PDE problems with Dirichlet BCs.
00736   FilterSingletons_ = List_in.get("schwarz: filter singletons", FilterSingletons_);
00737 
00738   // This copy may be needed by Amesos or other preconditioners.
00739   List_ = List_in;
00740 
00741   return(0);
00742 }
00743 
00744 //==============================================================================
00745 template<typename T>
00746 int Ifpack_AdditiveSchwarz<T>::Initialize()
00747 {
00748   IsInitialized_ = false;
00749   IsComputed_ = false; // values required
00750   Condest_ = -1.0; // zero-out condest
00751 
00752   if (Time_ == Teuchos::null)
00753     Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
00754 
00755   Time_->ResetStartTime();
00756 
00757   // compute the overlapping matrix if necessary
00758   if (IsOverlapping_) {
00759 #ifdef IFPACK_SUBCOMM_CODE
00760     if (NumMpiProcsPerSubdomain_ > 1) {
00761       OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, SubdomainId_) );
00762     } else {
00763       OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_));
00764     }
00765 #else
00766 #     ifdef IFPACK_NODE_AWARE_CODE
00767     int myNodeID;
00768     try{ myNodeID = List_.get("ML node id",-1);}
00769     catch(...){fprintf(stderr,"pid %d: no such entry (returned %d)\n",Comm().MyPID(),myNodeID);}
00770 /*
00771       cout << "pid " << Comm().MyPID()
00772            << ": calling Ifpack_OverlappingRowMatrix with myNodeID = "
00773            << myNodeID << ", OverlapLevel_ = " << OverlapLevel_ << endl;
00774 */
00775     OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, myNodeID) );
00776 #   else
00777     OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_) );
00778 #   endif
00779 #endif
00780 
00781     if (OverlappingMatrix_ == Teuchos::null) {
00782       IFPACK_CHK_ERR(-5);
00783     } 
00784   }
00785 
00786 # ifdef IFPACK_NODE_AWARE_CODE
00787 /*
00788   sleep(1);
00789   Comm().Barrier();
00790 */
00791 # endif
00792 
00793   IFPACK_CHK_ERR(Setup());
00794 
00795 # ifdef IFPACK_NODE_AWARE_CODE
00796 /*
00797   sleep(1);
00798   Comm().Barrier();
00799 */
00800 #endif
00801 
00802   if (Inverse_ == Teuchos::null)
00803     IFPACK_CHK_ERR(-5);
00804 
00805   if (LocalizedMatrix_ == Teuchos::null)
00806     IFPACK_CHK_ERR(-5);
00807 
00808   IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose()));
00809   IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
00810   IFPACK_CHK_ERR(Inverse_->Initialize());
00811 
00812   // Label is for Aztec-like solvers
00813   Label_ = "Ifpack_AdditiveSchwarz, ";
00814   if (UseTranspose())
00815     Label_ += ", transp";
00816   Label_ += ", ov = " + Ifpack_toString(OverlapLevel_)
00817     + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'";
00818 
00819   IsInitialized_ = true;
00820   ++NumInitialize_;
00821   InitializeTime_ += Time_->ElapsedTime();
00822 
00823 #ifdef IFPACK_FLOPCOUNTERS
00824   // count flops by summing up all the InitializeFlops() in each
00825   // Inverse. Each Inverse() can only give its flops -- it acts on one
00826   // process only
00827   double partial = Inverse_->InitializeFlops();
00828   double total;
00829   Comm().SumAll(&partial, &total, 1);
00830   InitializeFlops_ += total;
00831 #endif
00832 
00833   return(0);
00834 }
00835 
00836 //==============================================================================
00837 template<typename T>
00838 int Ifpack_AdditiveSchwarz<T>::Compute()
00839 {
00840   if (IsInitialized() == false)
00841     IFPACK_CHK_ERR(Initialize());
00842 
00843   Time_->ResetStartTime();
00844   IsComputed_ = false;
00845   Condest_ = -1.0;
00846   
00847   IFPACK_CHK_ERR(Inverse_->Compute());
00848 
00849   IsComputed_ = true; // need this here for Condest(Ifpack_Cheap)
00850   ++NumCompute_;
00851   ComputeTime_ += Time_->ElapsedTime();
00852 
00853 #ifdef IFPACK_FLOPCOUNTERS
00854   // sum up flops
00855   double partial = Inverse_->ComputeFlops();
00856    double total;
00857   Comm().SumAll(&partial, &total, 1);
00858   ComputeFlops_ += total;
00859 #endif
00860 
00861   // reset the Label
00862   string R = "";
00863   if (UseReordering_)
00864     R = ReorderingType_ + " reord, ";
00865 
00866   if (ComputeCondest_)
00867     Condest(Ifpack_Cheap);
00868   
00869   // add Condest() to label
00870   Label_ = "Ifpack_AdditiveSchwarz, ov = " + Ifpack_toString(OverlapLevel_)
00871     + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'"
00872     + "\n\t\t***** " + R + "Condition number estimate = "
00873     + Ifpack_toString(Condest_);
00874 
00875   return(0);
00876 }
00877 
00878 //==============================================================================
00879 template<typename T>
00880 int Ifpack_AdditiveSchwarz<T>::SetUseTranspose(bool UseTranspose_in)
00881 {
00882   // store the flag -- it will be set in Initialize() if Inverse_ does not
00883   // exist.
00884   UseTranspose_ = UseTranspose_in;
00885 
00886   // If Inverse_ exists, pass it right now.
00887   if (Inverse_!=Teuchos::null)
00888     IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
00889   return(0);
00890 }
00891 
00892 //==============================================================================
00893 template<typename T>
00894 int Ifpack_AdditiveSchwarz<T>::
00895 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00896 {
00897   IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
00898   return(0);
00899 }
00900 
00901 //==============================================================================
00902 template<typename T>
00903 double Ifpack_AdditiveSchwarz<T>::NormInf() const
00904 {
00905   return(-1.0);
00906 }
00907 
00908 //==============================================================================
00909 template<typename T>
00910 const char * Ifpack_AdditiveSchwarz<T>::Label() const
00911 {
00912   return(Label_.c_str());
00913 }
00914 
00915 //==============================================================================
00916 template<typename T>
00917 bool Ifpack_AdditiveSchwarz<T>::UseTranspose() const
00918 {
00919   return(UseTranspose_);
00920 }
00921 
00922 //==============================================================================
00923 template<typename T>
00924 bool Ifpack_AdditiveSchwarz<T>::HasNormInf() const
00925 {
00926   return(false);
00927 }
00928 
00929 //==============================================================================
00930 template<typename T>
00931 const Epetra_Comm & Ifpack_AdditiveSchwarz<T>::Comm() const
00932 {
00933   return(Matrix_->Comm());
00934 }
00935 
00936 //==============================================================================
00937 template<typename T>
00938 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorDomainMap() const
00939 {
00940   return(Matrix_->OperatorDomainMap());
00941 }
00942 
00943 //==============================================================================
00944 template<typename T>
00945 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorRangeMap() const
00946 {
00947   return(Matrix_->OperatorRangeMap());
00948 }
00949 
00950 //==============================================================================
00951 template<typename T>
00952 int Ifpack_AdditiveSchwarz<T>::
00953 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00954 {
00955   // compute the preconditioner is not done by the user
00956   if (!IsComputed())
00957     IFPACK_CHK_ERR(-3);
00958 
00959   int NumVectors = X.NumVectors();
00960 
00961   if (NumVectors != Y.NumVectors())
00962     IFPACK_CHK_ERR(-2); // wrong input
00963 
00964   Time_->ResetStartTime();
00965 
00966   Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
00967   Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
00968   Teuchos::RefCountPtr<Epetra_MultiVector> Xtmp;
00969 
00970   // for flop count, see bottom of this function
00971 #ifdef IFPACK_FLOPCOUNTERS
00972   double pre_partial = Inverse_->ApplyInverseFlops();
00973   double pre_total;
00974   Comm().SumAll(&pre_partial, &pre_total, 1);
00975 #endif
00976 
00977   // process overlap, may need to create vectors and import data
00978   if (IsOverlapping()) {
00979 #ifdef IFPACK_SUBCOMM_CODE
00980     if (OverlappingX == Teuchos::null) {
00981       OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), X.NumVectors()) );
00982       if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
00983     } else assert(OverlappingX->NumVectors() == X.NumVectors());
00984     if (OverlappingY == Teuchos::null) {
00985       OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(), Y.NumVectors()) );
00986       if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
00987     } else assert(OverlappingY->NumVectors() == Y.NumVectors());
00988 #else
00989 #   ifdef IFPACK_NODE_AWARE_CODE
00990     if (OverlappingX == Teuchos::null) {
00991       OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00992                                    X.NumVectors()) );
00993       if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
00994     } else assert(OverlappingX->NumVectors() == X.NumVectors());
00995     if (OverlappingY == Teuchos::null) {
00996       OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00997                                      Y.NumVectors()) );
00998       if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
00999     } else assert(OverlappingY->NumVectors() == Y.NumVectors());
01000 #else
01001     OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
01002                             X.NumVectors()) );
01003     OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
01004                             Y.NumVectors()) );
01005     if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
01006 #   endif
01007 #endif
01008     OverlappingY->PutScalar(0.0);
01009     OverlappingX->PutScalar(0.0);
01010     IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(X,*OverlappingX,Insert));
01011     // FIXME: this will not work with overlapping and non-zero starting
01012     // solutions. The same for other cases below.
01013     // IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(Y,*OverlappingY,Insert));
01014   }
01015   else {
01016     Xtmp = Teuchos::rcp( new Epetra_MultiVector(X) );
01017     OverlappingX = Xtmp;
01018     OverlappingY = Teuchos::rcp( &Y, false );
01019   }
01020 
01021   if (FilterSingletons_) {
01022     // process singleton filter
01023     Epetra_MultiVector ReducedX(SingletonFilter_->Map(),NumVectors);
01024     Epetra_MultiVector ReducedY(SingletonFilter_->Map(),NumVectors);
01025     IFPACK_CHK_ERR(SingletonFilter_->SolveSingletons(*OverlappingX,*OverlappingY));
01026     IFPACK_CHK_ERR(SingletonFilter_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX));
01027 
01028     // process reordering
01029     if (!UseReordering_) {
01030       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReducedX,ReducedY));
01031     }
01032     else {
01033       Epetra_MultiVector ReorderedX(ReducedX);
01034       Epetra_MultiVector ReorderedY(ReducedY);
01035       IFPACK_CHK_ERR(Reordering_->P(ReducedX,ReorderedX));
01036       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
01037       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,ReducedY));
01038     }
01039 
01040     // finish up with singletons
01041     IFPACK_CHK_ERR(SingletonFilter_->UpdateLHS(ReducedY,*OverlappingY));
01042   }
01043   else {
01044     // process reordering
01045     if (!UseReordering_) {
01046 #ifdef IFPACK_SUBCOMM_CODE
01047       tempX_.reset(&((*RangeVectorReindexer_)(*OverlappingX)), false);
01048       tempY_.reset(&((*DomainVectorReindexer_)(*OverlappingY)), false);
01049       IFPACK_CHK_ERR(Inverse_->ApplyInverse(*tempX_,*tempY_));
01050 #else
01051       IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
01052 #endif
01053     }
01054     else {
01055       Epetra_MultiVector ReorderedX(*OverlappingX);
01056       Epetra_MultiVector ReorderedY(*OverlappingY);
01057       IFPACK_CHK_ERR(Reordering_->P(*OverlappingX,ReorderedX));
01058       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
01059       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,*OverlappingY));
01060     }
01061   }
01062 
01063   if (IsOverlapping()) {
01064     IFPACK_CHK_ERR(OverlappingMatrix_->ExportMultiVector(*OverlappingY,Y,
01065                              CombineMode_));
01066   }
01067 
01068 #ifdef IFPACK_FLOPCOUNTERS
01069   // add flops. Note the we only have to add the newly counted
01070   // flops -- and each Inverse returns the cumulative sum
01071   double partial = Inverse_->ApplyInverseFlops();
01072   double total;
01073   Comm().SumAll(&partial, &total, 1);
01074   ApplyInverseFlops_ += total - pre_total;
01075 #endif
01076 
01077   // FIXME: right now I am skipping the overlap and singletons
01078   ++NumApplyInverse_;
01079   ApplyInverseTime_ += Time_->ElapsedTime();
01080 
01081   return(0);
01082  
01083 }
01084 
01085 //==============================================================================
01086 template<typename T>
01087 std::ostream& Ifpack_AdditiveSchwarz<T>::
01088 Print(std::ostream& os) const
01089 {
01090 #ifdef IFPACK_FLOPCOUNTERS
01091   double IF = InitializeFlops();
01092   double CF = ComputeFlops();
01093   double AF = ApplyInverseFlops();
01094 
01095   double IFT = 0.0, CFT = 0.0, AFT = 0.0;
01096   if (InitializeTime() != 0.0)
01097     IFT = IF / InitializeTime();
01098   if (ComputeTime() != 0.0)
01099     CFT = CF / ComputeTime();
01100   if (ApplyInverseTime() != 0.0)
01101     AFT = AF / ApplyInverseTime();
01102 #endif
01103 
01104   if (Matrix().Comm().MyPID())
01105     return(os);
01106 
01107   os << endl;
01108   os << "================================================================================" << endl;
01109   os << "Ifpack_AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
01110   if (CombineMode_ == Insert)
01111     os << "Combine mode                          = Insert" << endl;
01112   else if (CombineMode_ == Add)
01113     os << "Combine mode                          = Add" << endl;
01114   else if (CombineMode_ == Zero)
01115     os << "Combine mode                          = Zero" << endl;
01116   else if (CombineMode_ == Average)
01117     os << "Combine mode                          = Average" << endl;
01118   else if (CombineMode_ == AbsMax)
01119     os << "Combine mode                          = AbsMax" << endl;
01120 
01121   os << "Condition number estimate             = " << Condest_ << endl;
01122   os << "Global number of rows                 = " << Matrix_->NumGlobalRows() << endl;
01123 
01124 #ifdef IFPACK_SUBCOMM_CODE
01125   os << endl;
01126   os << "================================================================================" << endl;
01127   os << "Subcommunicator stats" << endl;
01128   os << "Number of MPI processes in simulation: " << NumMpiProcs_ << endl;
01129   os << "Number of subdomains: " << NumSubdomains_ << endl;
01130   os << "Number of MPI processes per subdomain: " << NumMpiProcsPerSubdomain_ << endl;
01131 #endif
01132 
01133   os << endl;
01134   os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
01135   os << "-----           -------   --------------       ------------     --------" << endl;
01136   os << "Initialize()    "   << std::setw(5) << NumInitialize()
01137      << "  " << std::setw(15) << InitializeTime() 
01138 #ifdef IFPACK_FLOPCOUNTERS
01139      << "  " << std::setw(15) << 1.0e-6 * IF 
01140      << "  " << std::setw(15) << 1.0e-6 * IFT
01141 #endif
01142      << endl;
01143   os << "Compute()       "   << std::setw(5) << NumCompute() 
01144      << "  " << std::setw(15) << ComputeTime()
01145 #ifdef IFPACK_FLOPCOUNTERS
01146      << "  " << std::setw(15) << 1.0e-6 * CF
01147      << "  " << std::setw(15) << 1.0e-6 * CFT
01148 #endif
01149      << endl;
01150   os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
01151      << "  " << std::setw(15) << ApplyInverseTime()
01152 #ifdef IFPACK_FLOPCOUNTERS
01153      << "  " << std::setw(15) << 1.0e-6 * AF
01154      << "  " << std::setw(15) << 1.0e-6 * AFT
01155 #endif
01156      << endl;
01157   os << "================================================================================" << endl;
01158   os << endl;
01159 
01160   return(os);
01161 }
01162 
01163 #include "Ifpack_Condest.h"
01164 //==============================================================================
01165 template<typename T>
01166 double Ifpack_AdditiveSchwarz<T>::
01167 Condest(const Ifpack_CondestType CT, const int MaxIters, 
01168         const double Tol, Epetra_RowMatrix* Matrix_in)
01169 {
01170   if (!IsComputed()) // cannot compute right now
01171     return(-1.0);
01172 
01173   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
01174 
01175   return(Condest_);
01176 }
01177 
01178 #endif // IFPACK_ADDITIVESCHWARZ_H
 All Classes Files Functions Variables Enumerations Friends