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 #include "Ifpack_SingletonFilter.h"
00014 #include "Ifpack_ReorderFilter.h"
00015 #include "Ifpack_Utils.h"
00016 #include "Ifpack_OverlappingRowMatrix.h"
00017 #include "Epetra_CombineMode.h"
00018 #include "Epetra_MultiVector.h"
00019 #include "Epetra_Map.h"
00020 #include "Epetra_Comm.h"
00021 #include "Epetra_Time.h"
00022 #include "Epetra_LinearProblem.h"
00023 #include "Epetra_RowMatrix.h"
00024 #include "Epetra_CrsMatrix.h"
00025 #include "Teuchos_ParameterList.hpp"
00026 #include "Teuchos_RefCountPtr.hpp"
00027 
00028 #ifdef IFPACK_NODE_AWARE_CODE
00029 #include "EpetraExt_OperatorOut.h"
00030 #include "EpetraExt_RowMatrixOut.h"
00031 #include "EpetraExt_BlockMapOut.h"
00032 #endif
00033 
00034 #ifdef HAVE_IFPACK_AMESOS
00035   #include "Ifpack_AMDReordering.h"
00036 #endif
00037 
00038 
00040 
00093 template<typename T>
00094 class Ifpack_AdditiveSchwarz : public virtual Ifpack_Preconditioner {
00095       
00096 public:
00097 
00099 
00100 
00111   Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
00112              int OverlapLevel_in = 0);
00113   
00115   virtual ~Ifpack_AdditiveSchwarz() {};
00117 
00119 
00121 
00130     virtual int SetUseTranspose(bool UseTranspose_in);
00132   
00134 
00136 
00146     virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00147 
00149 
00160     virtual int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00161 
00163     virtual double NormInf() const;
00165   
00167 
00169     virtual const char * Label() const;
00170 
00172     virtual bool UseTranspose() const;
00173 
00175     virtual bool HasNormInf() const;
00176 
00178     virtual const Epetra_Comm & Comm() const;
00179 
00181     virtual const Epetra_Map & OperatorDomainMap() const;
00182 
00184     virtual const Epetra_Map & OperatorRangeMap() const;
00186 
00188   virtual bool IsInitialized() const
00189   {
00190     return(IsInitialized_);
00191   }
00192 
00194   virtual bool IsComputed() const
00195   {
00196     return(IsComputed_);
00197   }
00198 
00200 
00219   virtual int SetParameters(Teuchos::ParameterList& List);
00220 
00221   // @}
00222 
00223   // @{ Query methods
00224   
00226   virtual int Initialize();
00227 
00229   virtual int Compute();
00230 
00232   virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
00233                          const int MaxIters = 1550,
00234                          const double Tol = 1e-9,
00235              Epetra_RowMatrix* Matrix_in = 0);
00236 
00238   virtual double Condest() const
00239   {
00240     return(Condest_);
00241   }
00242 
00244   virtual const Epetra_RowMatrix& Matrix() const
00245   {
00246     return(*Matrix_);
00247   }
00248 
00250   virtual bool IsOverlapping() const
00251   {
00252     return(IsOverlapping_);
00253   }
00254 
00256   virtual std::ostream& Print(std::ostream&) const;
00257   
00258   virtual const T* Inverse() const
00259   {
00260     return(&*Inverse_);
00261   }
00262 
00264   virtual int NumInitialize() const
00265   {
00266     return(NumInitialize_);
00267   }
00268 
00270   virtual int NumCompute() const
00271   {
00272     return(NumCompute_);
00273   }
00274 
00276   virtual int NumApplyInverse() const
00277   {
00278     return(NumApplyInverse_);
00279   }
00280 
00282   virtual double InitializeTime() const
00283   {
00284     return(InitializeTime_);
00285   }
00286 
00288   virtual double ComputeTime() const
00289   {
00290     return(ComputeTime_);
00291   }
00292 
00294   virtual double ApplyInverseTime() const
00295   {
00296     return(ApplyInverseTime_);
00297   }
00298 
00300   virtual double InitializeFlops() const
00301   {
00302     return(InitializeFlops_);
00303   }
00304 
00305   virtual double ComputeFlops() const
00306   {
00307     return(ComputeFlops_);
00308   }
00309 
00310   virtual double ApplyInverseFlops() const
00311   {
00312     return(ApplyInverseFlops_);
00313   }
00314 
00316   virtual int OverlapLevel() const
00317   {
00318     return(OverlapLevel_);
00319   }
00320 
00322   virtual const Teuchos::ParameterList& List() const
00323   {
00324     return(List_);
00325   }
00326 
00327 protected:
00328 
00329   // @}
00330 
00331   // @{ Internal merhods.
00332   
00334   Ifpack_AdditiveSchwarz(const Ifpack_AdditiveSchwarz& RHS)
00335   { }
00336 
00338   int Setup();
00339   
00340   // @}
00341 
00342   // @{ Internal data.
00343   
00345   Teuchos::RefCountPtr<const Epetra_RowMatrix> Matrix_;
00347   Teuchos::RefCountPtr<Ifpack_OverlappingRowMatrix> OverlappingMatrix_;
00349 /*
00350   //TODO if we choose to expose the node aware code, i.e., no ifdefs,
00351   //TODO then we should switch to this definition.
00352   Teuchos::RefCountPtr<Epetra_RowMatrix> LocalizedMatrix_;
00353 */
00354 # ifdef IFPACK_NODE_AWARE_CODE
00355   Teuchos::RefCountPtr<Ifpack_NodeFilter> LocalizedMatrix_;
00356 # else
00357   Teuchos::RefCountPtr<Ifpack_LocalFilter> LocalizedMatrix_;
00358 # endif
00359 
00360   string Label_;
00362   bool IsInitialized_;
00364   bool IsComputed_;
00366   bool UseTranspose_;
00368   bool IsOverlapping_;
00370   int OverlapLevel_;
00372   Teuchos::ParameterList List_;
00374   Epetra_CombineMode CombineMode_;
00376   double Condest_;
00378   bool ComputeCondest_;
00380   bool UseReordering_;
00382   string ReorderingType_;
00384   Teuchos::RefCountPtr<Ifpack_Reordering> Reordering_;
00386   Teuchos::RefCountPtr<Ifpack_ReorderFilter> ReorderedLocalizedMatrix_;
00388   bool FilterSingletons_;
00390   Teuchos::RefCountPtr<Ifpack_SingletonFilter> SingletonFilter_;
00392   int NumInitialize_;
00394   int NumCompute_;
00396   mutable int NumApplyInverse_;
00398   double InitializeTime_;
00400   double ComputeTime_;
00402   mutable double ApplyInverseTime_;
00404   double InitializeFlops_;
00406   double ComputeFlops_;
00408   mutable double ApplyInverseFlops_;
00410   Teuchos::RefCountPtr<Epetra_Time> Time_;
00412   Teuchos::RefCountPtr<T> Inverse_;
00414 # ifdef IFPACK_NODE_AWARE_CODE
00415   mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
00416   mutable Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
00417 #endif
00418 
00419 }; // class Ifpack_AdditiveSchwarz<T>
00420 
00421 //==============================================================================
00422 template<typename T>
00423 Ifpack_AdditiveSchwarz<T>::
00424 Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix_in,
00425                int OverlapLevel_in) :
00426   IsInitialized_(false),
00427   IsComputed_(false),
00428   UseTranspose_(false),
00429   IsOverlapping_(false),
00430   OverlapLevel_(OverlapLevel_in),
00431   CombineMode_(Zero),
00432   Condest_(-1.0),
00433   ComputeCondest_(true),
00434   UseReordering_(false),
00435   ReorderingType_("none"),
00436   FilterSingletons_(false),
00437   NumInitialize_(0),
00438   NumCompute_(0),
00439   NumApplyInverse_(0),
00440   InitializeTime_(0.0),
00441   ComputeTime_(0.0),
00442   ApplyInverseTime_(0.0),
00443   InitializeFlops_(0.0),
00444   ComputeFlops_(0.0),
00445   ApplyInverseFlops_(0.0)
00446 {
00447   // Construct a reference-counted pointer with the input matrix, don't manage the memory.
00448   Matrix_ = Teuchos::rcp( Matrix_in, false );
00449 
00450   if (Matrix_->Comm().NumProc() == 1)
00451     OverlapLevel_ = 0;
00452 
00453   if ((OverlapLevel_ != 0) && (Matrix_->Comm().NumProc() > 1))
00454     IsOverlapping_ = true;
00455   // Sets parameters to default values
00456   Teuchos::ParameterList List_in;
00457   SetParameters(List_in);
00458 }
00459 
00460 #ifdef IFPACK_NODE_AWARE_CODE
00461 extern int ML_NODE_ID;
00462 #endif
00463 
00464 //==============================================================================
00465 template<typename T>
00466 int Ifpack_AdditiveSchwarz<T>::Setup()
00467 {
00468 
00469   Epetra_RowMatrix* MatrixPtr;
00470 
00471 # ifdef IFPACK_NODE_AWARE_CODE
00472 /*
00473   sleep(3);
00474   if (Comm().MyPID() == 0) cout << "Printing out ovArowmap" << endl;
00475   Comm().Barrier();
00476 
00477   EpetraExt::BlockMapToMatrixMarketFile("ovArowmap",OverlappingMatrix_->RowMatrixRowMap());
00478   if (Comm().MyPID() == 0) cout << "Printing out ovAcolmap" << endl;
00479   Comm().Barrier();
00480   EpetraExt::BlockMapToMatrixMarketFile("ovAcolmap",OverlappingMatrix_->RowMatrixColMap());
00481   Comm().Barrier();
00482 */
00483 /*
00484   EpetraExt::RowMatrixToMatlabFile("ovA",*OverlappingMatrix_);
00485   fprintf(stderr,"p %d n %d matrix file done\n",Comm().MyPID(),ML_NODE_ID);
00486   Comm().Barrier();
00487 */
00488   int nodeID;
00489   try{ nodeID = List_.get("ML node id",0);}
00490   catch(...){fprintf(stderr,"%s","Ifpack_AdditiveSchwarz<T>::Setup(): no parameter \"ML node id\"\n\n");
00491              cout << List_ << endl;}
00492 # endif
00493 
00494   try{
00495   if (OverlappingMatrix_ != Teuchos::null)
00496   {
00497 #   ifdef IFPACK_NODE_AWARE_CODE
00498     Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(OverlappingMatrix_,nodeID); //FIXME
00499     LocalizedMatrix_ = Teuchos::rcp(tt);
00500     //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
00501 #   else
00502     LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(OverlappingMatrix_) );
00503 #   endif
00504   }
00505   else
00506   {
00507 #   ifdef IFPACK_NODE_AWARE_CODE
00508     Ifpack_NodeFilter *tt = new Ifpack_NodeFilter(Matrix_,nodeID); //FIXME
00509     LocalizedMatrix_ = Teuchos::rcp(tt);
00510     //LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
00511 #   else
00512     LocalizedMatrix_ = Teuchos::rcp( new Ifpack_LocalFilter(Matrix_) );
00513 #   endif
00514   }
00515   }
00516   catch(...) {
00517      fprintf(stderr,"%s","AdditiveSchwarz Setup: problem creating local filter matrix.\n");
00518   }
00519 
00520   if (LocalizedMatrix_ == Teuchos::null)
00521     IFPACK_CHK_ERR(-5);
00522 
00523   // users may want to skip singleton check
00524   if (FilterSingletons_) {
00525     SingletonFilter_ = Teuchos::rcp( new Ifpack_SingletonFilter(LocalizedMatrix_) );
00526     MatrixPtr = &*SingletonFilter_;
00527   }
00528   else
00529     MatrixPtr = &*LocalizedMatrix_;
00530 
00531   if (UseReordering_) {
00532 
00533     // create reordering and compute it
00534     if (ReorderingType_ == "rcm")
00535       Reordering_ = Teuchos::rcp( new Ifpack_RCMReordering() );
00536     else if (ReorderingType_ == "metis")
00537       Reordering_ = Teuchos::rcp( new Ifpack_METISReordering() );
00538 #ifdef HAVE_IFPACK_AMESOS   
00539     else if (ReorderingType_ == "amd" )
00540       Reordering_ = Teuchos::rcp( new Ifpack_AMDReordering() );
00541 #endif
00542     else {
00543       cerr << "reordering type not correct (" << ReorderingType_ << ")" << endl;
00544       exit(EXIT_FAILURE);
00545     }
00546     if (Reordering_ == Teuchos::null) IFPACK_CHK_ERR(-5);
00547 
00548     IFPACK_CHK_ERR(Reordering_->SetParameters(List_));
00549     IFPACK_CHK_ERR(Reordering_->Compute(*MatrixPtr));
00550 
00551     // now create reordered localized matrix
00552     ReorderedLocalizedMatrix_ = 
00553       Teuchos::rcp( new Ifpack_ReorderFilter(Teuchos::rcp( MatrixPtr, false ), Reordering_) );
00554 
00555     if (ReorderedLocalizedMatrix_ == Teuchos::null) IFPACK_CHK_ERR(-5);
00556 
00557     MatrixPtr = &*ReorderedLocalizedMatrix_;
00558   }
00559 
00560   Inverse_ = Teuchos::rcp( new T(MatrixPtr) );
00561 
00562   if (Inverse_ == Teuchos::null)
00563     IFPACK_CHK_ERR(-5);
00564 
00565   return(0);
00566 }
00567 
00568 //==============================================================================
00569 template<typename T>
00570 int Ifpack_AdditiveSchwarz<T>::SetParameters(Teuchos::ParameterList& List_in)
00571 {
00572  
00573   // compute the condition number each time Compute() is invoked.
00574   ComputeCondest_ = List_in.get("schwarz: compute condest", ComputeCondest_);
00575   // combine mode
00576   if( Teuchos::ParameterEntry *combineModeEntry = List_in.getEntryPtr("schwarz: combine mode") )
00577   {
00578     if( typeid(std::string) == combineModeEntry->getAny().type() )
00579     {
00580       std::string mode = List_in.get("schwarz: combine mode", "Add");
00581       if (mode == "Add")
00582         CombineMode_ = Add;
00583       else if (mode == "Zero")
00584         CombineMode_ = Zero;
00585       else if (mode == "Insert")
00586         CombineMode_ = Insert;
00587       else if (mode == "InsertAdd")
00588         CombineMode_ = InsertAdd;
00589       else if (mode == "Average")
00590         CombineMode_ = Average;
00591       else if (mode == "AbsMax")
00592         CombineMode_ = AbsMax;
00593       else
00594       {
00595         TEST_FOR_EXCEPTION(
00596           true,std::logic_error
00597           ,"Error, The (Epetra) combine mode of \""<<mode<<"\" is not valid!  Only the values"
00598           " \"Add\", \"Zero\", \"Insert\", \"InsertAdd\", \"Average\", and \"AbsMax\" are accepted!"
00599           );
00600       }
00601     }
00602     else if ( typeid(Epetra_CombineMode) == combineModeEntry->getAny().type() )
00603     {
00604       CombineMode_ = Teuchos::any_cast<Epetra_CombineMode>(combineModeEntry->getAny());
00605     }
00606     else
00607     {
00608       // Throw exception with good error message!
00609       Teuchos::getParameter<std::string>(List_in,"schwarz: combine mode");
00610     }
00611   }
00612   else
00613   {
00614     // Make the default be a string to be consistent with the valid parameters!
00615     List_in.get("schwarz: combine mode","Zero");
00616   }
00617   // type of reordering
00618   ReorderingType_ = List_in.get("schwarz: reordering type", ReorderingType_);
00619   if (ReorderingType_ == "none")
00620     UseReordering_ = false;
00621   else 
00622     UseReordering_ = true;
00623   // if true, filter singletons. NOTE: the filtered matrix can still have
00624   // singletons! A simple example: upper triangular matrix, if I remove
00625   // the lower node, I still get a matrix with a singleton! However, filter
00626   // singletons should help for PDE problems with Dirichlet BCs.
00627   FilterSingletons_ = List_in.get("schwarz: filter singletons", FilterSingletons_);
00628 
00629   // This copy may be needed by Amesos or other preconditioners.
00630   List_ = List_in;
00631 
00632   return(0);
00633 }
00634 
00635 //==============================================================================
00636 template<typename T>
00637 int Ifpack_AdditiveSchwarz<T>::Initialize()
00638 {
00639   IsInitialized_ = false;
00640   IsComputed_ = false; // values required
00641   Condest_ = -1.0; // zero-out condest
00642 
00643   if (Time_ == Teuchos::null)
00644     Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
00645 
00646   Time_->ResetStartTime();
00647 
00648   // compute the overlapping matrix if necessary
00649   if (IsOverlapping_) {
00650 #     ifdef IFPACK_NODE_AWARE_CODE
00651       int myNodeID;
00652       try{ myNodeID = List_.get("ML node id",-1);}
00653       catch(...){fprintf(stderr,"pid %d: no such entry (returned %d)\n",Comm().MyPID(),myNodeID);}
00654 /*
00655       cout << "pid " << Comm().MyPID()
00656            << ": calling Ifpack_OverlappingRowMatrix with myNodeID = "
00657            << myNodeID << ", OverlapLevel_ = " << OverlapLevel_ << endl;
00658 */
00659       OverlappingMatrix_ = Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_, myNodeID) );
00660 #   else
00661       OverlappingMatrix_ =
00662         Teuchos::rcp( new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_) );
00663 #   endif
00664 
00665     if (OverlappingMatrix_ == Teuchos::null) {
00666       IFPACK_CHK_ERR(-5);
00667     } 
00668   }
00669 
00670 # ifdef IFPACK_NODE_AWARE_CODE
00671 /*
00672   sleep(1);
00673   Comm().Barrier();
00674 */
00675 # endif
00676 
00677   IFPACK_CHK_ERR(Setup());
00678 
00679 # ifdef IFPACK_NODE_AWARE_CODE
00680 /*
00681   sleep(1);
00682   Comm().Barrier();
00683 */
00684 #endif
00685 
00686   if (Inverse_ == Teuchos::null)
00687     IFPACK_CHK_ERR(-5);
00688 
00689   if (LocalizedMatrix_ == Teuchos::null)
00690     IFPACK_CHK_ERR(-5);
00691 
00692   IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose()));
00693   IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
00694   IFPACK_CHK_ERR(Inverse_->Initialize());
00695 
00696   // Label is for Aztec-like solvers
00697   Label_ = "Ifpack_AdditiveSchwarz, ";
00698   if (UseTranspose())
00699     Label_ += ", transp";
00700   Label_ += ", ov = " + Ifpack_toString(OverlapLevel_)
00701     + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'";
00702 
00703   IsInitialized_ = true;
00704   ++NumInitialize_;
00705   InitializeTime_ += Time_->ElapsedTime();
00706 
00707   // count flops by summing up all the InitializeFlops() in each
00708   // Inverse. Each Inverse() can only give its flops -- it acts on one
00709   // process only
00710   double partial = Inverse_->InitializeFlops();
00711   double total;
00712   Comm().SumAll(&partial, &total, 1);
00713   InitializeFlops_ += total;
00714 
00715   return(0);
00716 }
00717 
00718 //==============================================================================
00719 template<typename T>
00720 int Ifpack_AdditiveSchwarz<T>::Compute()
00721 {
00722 
00723   if (IsInitialized() == false)
00724     IFPACK_CHK_ERR(Initialize());
00725 
00726   Time_->ResetStartTime();
00727   IsComputed_ = false;
00728   Condest_ = -1.0;
00729   
00730   IFPACK_CHK_ERR(Inverse_->Compute());
00731 
00732   IsComputed_ = true; // need this here for Condest(Ifpack_Cheap)
00733   ++NumCompute_;
00734   ComputeTime_ += Time_->ElapsedTime();
00735 
00736   // sum up flops
00737   double partial = Inverse_->ComputeFlops();
00738    double total;
00739   Comm().SumAll(&partial, &total, 1);
00740   ComputeFlops_ += total;
00741 
00742   // reset the Label
00743   string R = "";
00744   if (UseReordering_)
00745     R = ReorderingType_ + " reord, ";
00746 
00747   if (ComputeCondest_)
00748     Condest(Ifpack_Cheap);
00749   
00750   // add Condest() to label
00751   Label_ = "Ifpack_AdditiveSchwarz, ov = " + Ifpack_toString(OverlapLevel_)
00752     + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'"
00753     + "\n\t\t***** " + R + "Condition number estimate = "
00754     + Ifpack_toString(Condest());
00755 
00756   return(0);
00757 }
00758 
00759 //==============================================================================
00760 template<typename T>
00761 int Ifpack_AdditiveSchwarz<T>::SetUseTranspose(bool UseTranspose_in)
00762 {
00763   // store the flag -- it will be set in Initialize() if Inverse_ does not
00764   // exist.
00765   UseTranspose_ = UseTranspose_in;
00766 
00767   // If Inverse_ exists, pass it right now.
00768   if (Inverse_!=Teuchos::null)
00769     IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
00770   return(0);
00771 }
00772 
00773 //==============================================================================
00774 template<typename T>
00775 int Ifpack_AdditiveSchwarz<T>::
00776 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00777 {
00778   IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
00779   return(0);
00780 }
00781 
00782 //==============================================================================
00783 template<typename T>
00784 double Ifpack_AdditiveSchwarz<T>::NormInf() const
00785 {
00786   return(-1.0);
00787 }
00788 
00789 //==============================================================================
00790 template<typename T>
00791 const char * Ifpack_AdditiveSchwarz<T>::Label() const
00792 {
00793   return(Label_.c_str());
00794 }
00795 
00796 //==============================================================================
00797 template<typename T>
00798 bool Ifpack_AdditiveSchwarz<T>::UseTranspose() const
00799 {
00800   return(UseTranspose_);
00801 }
00802 
00803 //==============================================================================
00804 template<typename T>
00805 bool Ifpack_AdditiveSchwarz<T>::HasNormInf() const
00806 {
00807   return(false);
00808 }
00809 
00810 //==============================================================================
00811 template<typename T>
00812 const Epetra_Comm & Ifpack_AdditiveSchwarz<T>::Comm() const
00813 {
00814   return(Matrix_->Comm());
00815 }
00816 
00817 //==============================================================================
00818 template<typename T>
00819 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorDomainMap() const
00820 {
00821   return(Matrix_->OperatorDomainMap());
00822 }
00823 
00824 //==============================================================================
00825 template<typename T>
00826 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorRangeMap() const
00827 {
00828   return(Matrix_->OperatorRangeMap());
00829 }
00830 
00831 //==============================================================================
00832 template<typename T>
00833 int Ifpack_AdditiveSchwarz<T>::
00834 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00835 {
00836   // compute the preconditioner is not done by the user
00837   if (!IsComputed())
00838     IFPACK_CHK_ERR(-3);
00839 
00840   int NumVectors = X.NumVectors();
00841 
00842   if (NumVectors != Y.NumVectors())
00843     IFPACK_CHK_ERR(-2); // wrong input
00844 
00845   Time_->ResetStartTime();
00846 
00847   Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
00848   Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
00849   Teuchos::RefCountPtr<Epetra_MultiVector> Xtmp;
00850 
00851   // for flop count, see bottom of this function
00852   double pre_partial = Inverse_->ApplyInverseFlops();
00853   double pre_total;
00854   Comm().SumAll(&pre_partial, &pre_total, 1);
00855 
00856   // process overlap, may need to create vectors and import data
00857   if (IsOverlapping()) {
00858 #   ifdef IFPACK_NODE_AWARE_CODE
00859     if (OverlappingX == Teuchos::null) {
00860       OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00861                                    X.NumVectors()) );
00862       if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
00863     } else assert(OverlappingX->NumVectors() == X.NumVectors());
00864     if (OverlappingY == Teuchos::null) {
00865       OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00866                                      Y.NumVectors()) );
00867       if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
00868     } else assert(OverlappingY->NumVectors() == Y.NumVectors());
00869 #else
00870     OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00871                             X.NumVectors()) );
00872     OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00873                             Y.NumVectors()) );
00874     if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
00875 #   endif
00876     OverlappingY->PutScalar(0.0);
00877     OverlappingX->PutScalar(0.0);
00878     IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(X,*OverlappingX,Insert));
00879     // FIXME: this will not work with overlapping and non-zero starting
00880     // solutions. The same for other cases below.
00881     // IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(Y,*OverlappingY,Insert));
00882   }
00883   else {
00884     Xtmp = Teuchos::rcp( new Epetra_MultiVector(X) );
00885     OverlappingX = Xtmp;
00886     OverlappingY = Teuchos::rcp( &Y, false );
00887   }
00888 
00889   if (FilterSingletons_) {
00890     // process singleton filter
00891     Epetra_MultiVector ReducedX(SingletonFilter_->Map(),NumVectors);
00892     Epetra_MultiVector ReducedY(SingletonFilter_->Map(),NumVectors);
00893     IFPACK_CHK_ERR(SingletonFilter_->SolveSingletons(*OverlappingX,*OverlappingY));
00894     IFPACK_CHK_ERR(SingletonFilter_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX));
00895 
00896     // process reordering
00897     if (!UseReordering_) {
00898       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReducedX,ReducedY));
00899     }
00900     else {
00901       Epetra_MultiVector ReorderedX(ReducedX);
00902       Epetra_MultiVector ReorderedY(ReducedY);
00903       IFPACK_CHK_ERR(Reordering_->P(ReducedX,ReorderedX));
00904       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
00905       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,ReducedY));
00906     }
00907 
00908     // finish up with singletons
00909     IFPACK_CHK_ERR(SingletonFilter_->UpdateLHS(ReducedY,*OverlappingY));
00910   }
00911   else {
00912     // process reordering
00913     if (!UseReordering_) {
00914       IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
00915     }
00916     else {
00917       Epetra_MultiVector ReorderedX(*OverlappingX);
00918       Epetra_MultiVector ReorderedY(*OverlappingY);
00919       IFPACK_CHK_ERR(Reordering_->P(*OverlappingX,ReorderedX));
00920       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
00921       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,*OverlappingY));
00922     }
00923   }
00924 
00925   if (IsOverlapping()) {
00926     IFPACK_CHK_ERR(OverlappingMatrix_->ExportMultiVector(*OverlappingY,Y,
00927                              CombineMode_));
00928   }
00929 
00930   // add flops. Note the we only have to add the newly counted
00931   // flops -- and each Inverse returns the cumulative sum
00932   double partial = Inverse_->ApplyInverseFlops();
00933   double total;
00934   Comm().SumAll(&partial, &total, 1);
00935   ApplyInverseFlops_ += total - pre_total;
00936 
00937   // FIXME: right now I am skipping the overlap and singletons
00938   ++NumApplyInverse_;
00939   ApplyInverseTime_ += Time_->ElapsedTime();
00940 
00941   return(0);
00942  
00943 }
00944 
00945 //==============================================================================
00946 template<typename T>
00947 std::ostream& Ifpack_AdditiveSchwarz<T>::
00948 Print(std::ostream& os) const
00949 {
00950   double IF = InitializeFlops();
00951   double CF = ComputeFlops();
00952   double AF = ApplyInverseFlops();
00953 
00954   double IFT = 0.0, CFT = 0.0, AFT = 0.0;
00955   if (InitializeTime() != 0.0)
00956     IFT = IF / InitializeTime();
00957   if (ComputeTime() != 0.0)
00958     CFT = CF / ComputeTime();
00959   if (ApplyInverseTime() != 0.0)
00960     AFT = AF / ApplyInverseTime();
00961 
00962   if (Matrix().Comm().MyPID())
00963     return(os);
00964 
00965   os << endl;
00966   os << "================================================================================" << endl;
00967   os << "Ifpack_AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
00968   if (CombineMode_ == Insert)
00969     os << "Combine mode                          = Insert" << endl;
00970   else if (CombineMode_ == Add)
00971     os << "Combine mode                          = Add" << endl;
00972   else if (CombineMode_ == Zero)
00973     os << "Combine mode                          = Zero" << endl;
00974   else if (CombineMode_ == Average)
00975     os << "Combine mode                          = Average" << endl;
00976   else if (CombineMode_ == AbsMax)
00977     os << "Combine mode                          = AbsMax" << endl;
00978 
00979   os << "Condition number estimate             = " << Condest_ << endl;
00980   os << "Global number of rows                 = " << Matrix_->NumGlobalRows() << endl;
00981 
00982   os << endl;
00983   os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00984   os << "-----           -------   --------------       ------------     --------" << endl;
00985   os << "Initialize()    "   << std::setw(5) << NumInitialize()
00986      << "  " << std::setw(15) << InitializeTime() 
00987      << "  " << std::setw(15) << 1.0e-6 * IF 
00988      << "  " << std::setw(15) << 1.0e-6 * IFT << endl;
00989   os << "Compute()       "   << std::setw(5) << NumCompute() 
00990      << "  " << std::setw(15) << ComputeTime()
00991      << "  " << std::setw(15) << 1.0e-6 * CF
00992      << "  " << std::setw(15) << 1.0e-6 * CFT << endl;
00993   os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00994      << "  " << std::setw(15) << ApplyInverseTime()
00995      << "  " << std::setw(15) << 1.0e-6 * AF
00996      << "  " << std::setw(15) << 1.0e-6 * AFT << endl;
00997   os << "================================================================================" << endl;
00998   os << endl;
00999 
01000   return(os);
01001 }
01002 
01003 #include "Ifpack_Condest.h"
01004 //==============================================================================
01005 template<typename T>
01006 double Ifpack_AdditiveSchwarz<T>::
01007 Condest(const Ifpack_CondestType CT, const int MaxIters, 
01008         const double Tol, Epetra_RowMatrix* Matrix_in)
01009 {
01010   if (!IsComputed()) // cannot compute right now
01011     return(-1.0);
01012 
01013   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
01014 
01015   return(Condest_);
01016 }
01017 
01018 #endif // IFPACK_ADDITIVESCHWARZ_H
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:22 2011 for IFPACK by  doxygen 1.6.3