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 #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 #ifdef IFPACK_FLOPCOUNTERS
00708   // count flops by summing up all the InitializeFlops() in each
00709   // Inverse. Each Inverse() can only give its flops -- it acts on one
00710   // process only
00711   double partial = Inverse_->InitializeFlops();
00712   double total;
00713   Comm().SumAll(&partial, &total, 1);
00714   InitializeFlops_ += total;
00715 #endif
00716 
00717   return(0);
00718 }
00719 
00720 //==============================================================================
00721 template<typename T>
00722 int Ifpack_AdditiveSchwarz<T>::Compute()
00723 {
00724 
00725   if (IsInitialized() == false)
00726     IFPACK_CHK_ERR(Initialize());
00727 
00728   Time_->ResetStartTime();
00729   IsComputed_ = false;
00730   Condest_ = -1.0;
00731   
00732   IFPACK_CHK_ERR(Inverse_->Compute());
00733 
00734   IsComputed_ = true; // need this here for Condest(Ifpack_Cheap)
00735   ++NumCompute_;
00736   ComputeTime_ += Time_->ElapsedTime();
00737 
00738 #ifdef IFPACK_FLOPCOUNTERS
00739   // sum up flops
00740   double partial = Inverse_->ComputeFlops();
00741    double total;
00742   Comm().SumAll(&partial, &total, 1);
00743   ComputeFlops_ += total;
00744 #endif
00745 
00746   // reset the Label
00747   string R = "";
00748   if (UseReordering_)
00749     R = ReorderingType_ + " reord, ";
00750 
00751   if (ComputeCondest_)
00752     Condest(Ifpack_Cheap);
00753   
00754   // add Condest() to label
00755   Label_ = "Ifpack_AdditiveSchwarz, ov = " + Ifpack_toString(OverlapLevel_)
00756     + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'"
00757     + "\n\t\t***** " + R + "Condition number estimate = "
00758     + Ifpack_toString(Condest());
00759 
00760   return(0);
00761 }
00762 
00763 //==============================================================================
00764 template<typename T>
00765 int Ifpack_AdditiveSchwarz<T>::SetUseTranspose(bool UseTranspose_in)
00766 {
00767   // store the flag -- it will be set in Initialize() if Inverse_ does not
00768   // exist.
00769   UseTranspose_ = UseTranspose_in;
00770 
00771   // If Inverse_ exists, pass it right now.
00772   if (Inverse_!=Teuchos::null)
00773     IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
00774   return(0);
00775 }
00776 
00777 //==============================================================================
00778 template<typename T>
00779 int Ifpack_AdditiveSchwarz<T>::
00780 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00781 {
00782   IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
00783   return(0);
00784 }
00785 
00786 //==============================================================================
00787 template<typename T>
00788 double Ifpack_AdditiveSchwarz<T>::NormInf() const
00789 {
00790   return(-1.0);
00791 }
00792 
00793 //==============================================================================
00794 template<typename T>
00795 const char * Ifpack_AdditiveSchwarz<T>::Label() const
00796 {
00797   return(Label_.c_str());
00798 }
00799 
00800 //==============================================================================
00801 template<typename T>
00802 bool Ifpack_AdditiveSchwarz<T>::UseTranspose() const
00803 {
00804   return(UseTranspose_);
00805 }
00806 
00807 //==============================================================================
00808 template<typename T>
00809 bool Ifpack_AdditiveSchwarz<T>::HasNormInf() const
00810 {
00811   return(false);
00812 }
00813 
00814 //==============================================================================
00815 template<typename T>
00816 const Epetra_Comm & Ifpack_AdditiveSchwarz<T>::Comm() const
00817 {
00818   return(Matrix_->Comm());
00819 }
00820 
00821 //==============================================================================
00822 template<typename T>
00823 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorDomainMap() const
00824 {
00825   return(Matrix_->OperatorDomainMap());
00826 }
00827 
00828 //==============================================================================
00829 template<typename T>
00830 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorRangeMap() const
00831 {
00832   return(Matrix_->OperatorRangeMap());
00833 }
00834 
00835 //==============================================================================
00836 template<typename T>
00837 int Ifpack_AdditiveSchwarz<T>::
00838 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00839 {
00840   // compute the preconditioner is not done by the user
00841   if (!IsComputed())
00842     IFPACK_CHK_ERR(-3);
00843 
00844   int NumVectors = X.NumVectors();
00845 
00846   if (NumVectors != Y.NumVectors())
00847     IFPACK_CHK_ERR(-2); // wrong input
00848 
00849   Time_->ResetStartTime();
00850 
00851   Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingX;
00852   Teuchos::RefCountPtr<Epetra_MultiVector> OverlappingY;
00853   Teuchos::RefCountPtr<Epetra_MultiVector> Xtmp;
00854 
00855   // for flop count, see bottom of this function
00856 #ifdef IFPACK_FLOPCOUNTERS
00857   double pre_partial = Inverse_->ApplyInverseFlops();
00858   double pre_total;
00859   Comm().SumAll(&pre_partial, &pre_total, 1);
00860 #endif
00861 
00862   // process overlap, may need to create vectors and import data
00863   if (IsOverlapping()) {
00864 #   ifdef IFPACK_NODE_AWARE_CODE
00865     if (OverlappingX == Teuchos::null) {
00866       OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00867                                    X.NumVectors()) );
00868       if (OverlappingX == Teuchos::null) IFPACK_CHK_ERR(-5);
00869     } else assert(OverlappingX->NumVectors() == X.NumVectors());
00870     if (OverlappingY == Teuchos::null) {
00871       OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00872                                      Y.NumVectors()) );
00873       if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
00874     } else assert(OverlappingY->NumVectors() == Y.NumVectors());
00875 #else
00876     OverlappingX = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00877                             X.NumVectors()) );
00878     OverlappingY = Teuchos::rcp( new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00879                             Y.NumVectors()) );
00880     if (OverlappingY == Teuchos::null) IFPACK_CHK_ERR(-5);
00881 #   endif
00882     OverlappingY->PutScalar(0.0);
00883     OverlappingX->PutScalar(0.0);
00884     IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(X,*OverlappingX,Insert));
00885     // FIXME: this will not work with overlapping and non-zero starting
00886     // solutions. The same for other cases below.
00887     // IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(Y,*OverlappingY,Insert));
00888   }
00889   else {
00890     Xtmp = Teuchos::rcp( new Epetra_MultiVector(X) );
00891     OverlappingX = Xtmp;
00892     OverlappingY = Teuchos::rcp( &Y, false );
00893   }
00894 
00895   if (FilterSingletons_) {
00896     // process singleton filter
00897     Epetra_MultiVector ReducedX(SingletonFilter_->Map(),NumVectors);
00898     Epetra_MultiVector ReducedY(SingletonFilter_->Map(),NumVectors);
00899     IFPACK_CHK_ERR(SingletonFilter_->SolveSingletons(*OverlappingX,*OverlappingY));
00900     IFPACK_CHK_ERR(SingletonFilter_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX));
00901 
00902     // process reordering
00903     if (!UseReordering_) {
00904       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReducedX,ReducedY));
00905     }
00906     else {
00907       Epetra_MultiVector ReorderedX(ReducedX);
00908       Epetra_MultiVector ReorderedY(ReducedY);
00909       IFPACK_CHK_ERR(Reordering_->P(ReducedX,ReorderedX));
00910       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
00911       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,ReducedY));
00912     }
00913 
00914     // finish up with singletons
00915     IFPACK_CHK_ERR(SingletonFilter_->UpdateLHS(ReducedY,*OverlappingY));
00916   }
00917   else {
00918     // process reordering
00919     if (!UseReordering_) {
00920       IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
00921     }
00922     else {
00923       Epetra_MultiVector ReorderedX(*OverlappingX);
00924       Epetra_MultiVector ReorderedY(*OverlappingY);
00925       IFPACK_CHK_ERR(Reordering_->P(*OverlappingX,ReorderedX));
00926       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
00927       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,*OverlappingY));
00928     }
00929   }
00930 
00931   if (IsOverlapping()) {
00932     IFPACK_CHK_ERR(OverlappingMatrix_->ExportMultiVector(*OverlappingY,Y,
00933                              CombineMode_));
00934   }
00935 
00936 #ifdef IFPACK_FLOPCOUNTERS
00937   // add flops. Note the we only have to add the newly counted
00938   // flops -- and each Inverse returns the cumulative sum
00939   double partial = Inverse_->ApplyInverseFlops();
00940   double total;
00941   Comm().SumAll(&partial, &total, 1);
00942   ApplyInverseFlops_ += total - pre_total;
00943 #endif
00944 
00945   // FIXME: right now I am skipping the overlap and singletons
00946   ++NumApplyInverse_;
00947   ApplyInverseTime_ += Time_->ElapsedTime();
00948 
00949   return(0);
00950  
00951 }
00952 
00953 //==============================================================================
00954 template<typename T>
00955 std::ostream& Ifpack_AdditiveSchwarz<T>::
00956 Print(std::ostream& os) const
00957 {
00958 #ifdef IFPACK_FLOPCOUNTERS
00959   double IF = InitializeFlops();
00960   double CF = ComputeFlops();
00961   double AF = ApplyInverseFlops();
00962 
00963   double IFT = 0.0, CFT = 0.0, AFT = 0.0;
00964   if (InitializeTime() != 0.0)
00965     IFT = IF / InitializeTime();
00966   if (ComputeTime() != 0.0)
00967     CFT = CF / ComputeTime();
00968   if (ApplyInverseTime() != 0.0)
00969     AFT = AF / ApplyInverseTime();
00970 #endif
00971 
00972   if (Matrix().Comm().MyPID())
00973     return(os);
00974 
00975   os << endl;
00976   os << "================================================================================" << endl;
00977   os << "Ifpack_AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
00978   if (CombineMode_ == Insert)
00979     os << "Combine mode                          = Insert" << endl;
00980   else if (CombineMode_ == Add)
00981     os << "Combine mode                          = Add" << endl;
00982   else if (CombineMode_ == Zero)
00983     os << "Combine mode                          = Zero" << endl;
00984   else if (CombineMode_ == Average)
00985     os << "Combine mode                          = Average" << endl;
00986   else if (CombineMode_ == AbsMax)
00987     os << "Combine mode                          = AbsMax" << endl;
00988 
00989   os << "Condition number estimate             = " << Condest_ << endl;
00990   os << "Global number of rows                 = " << Matrix_->NumGlobalRows() << endl;
00991 
00992   os << endl;
00993   os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00994   os << "-----           -------   --------------       ------------     --------" << endl;
00995   os << "Initialize()    "   << std::setw(5) << NumInitialize()
00996      << "  " << std::setw(15) << InitializeTime() 
00997 #ifdef IFPACK_FLOPCOUNTERS
00998      << "  " << std::setw(15) << 1.0e-6 * IF 
00999      << "  " << std::setw(15) << 1.0e-6 * IFT
01000 #endif
01001      << endl;
01002   os << "Compute()       "   << std::setw(5) << NumCompute() 
01003      << "  " << std::setw(15) << ComputeTime()
01004 #ifdef IFPACK_FLOPCOUNTERS
01005      << "  " << std::setw(15) << 1.0e-6 * CF
01006      << "  " << std::setw(15) << 1.0e-6 * CFT
01007 #endif
01008      << endl;
01009   os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
01010      << "  " << std::setw(15) << ApplyInverseTime()
01011 #ifdef IFPACK_FLOPCOUNTERS
01012      << "  " << std::setw(15) << 1.0e-6 * AF
01013      << "  " << std::setw(15) << 1.0e-6 * AFT
01014 #endif
01015      << endl;
01016   os << "================================================================================" << endl;
01017   os << endl;
01018 
01019   return(os);
01020 }
01021 
01022 #include "Ifpack_Condest.h"
01023 //==============================================================================
01024 template<typename T>
01025 double Ifpack_AdditiveSchwarz<T>::
01026 Condest(const Ifpack_CondestType CT, const int MaxIters, 
01027         const double Tol, Epetra_RowMatrix* Matrix_in)
01028 {
01029   if (!IsComputed()) // cannot compute right now
01030     return(-1.0);
01031 
01032   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
01033 
01034   return(Condest_);
01035 }
01036 
01037 #endif // IFPACK_ADDITIVESCHWARZ_H
 All Classes Files Functions Variables Enumerations Friends