Ifpack_AdditiveSchwarz.h

00001 #ifndef IFPACK_ADDITIVESCHWARZ_H
00002 #define IFPACK_ADDITIVESCHWARZ_H
00003 
00004 #include "Ifpack_ConfigDefs.h"
00005 #if defined(HAVE_IFPACK_TEUCHOS)
00006 #include "Ifpack_Preconditioner.h"
00007 #include "Teuchos_ParameterList.hpp"
00008 #include "Ifpack_ConfigDefs.h"
00009 #include "Ifpack_Preconditioner.h"
00010 #include "Ifpack_Reordering.h"
00011 #include "Ifpack_RCMReordering.h"
00012 #include "Ifpack_METISReordering.h"
00013 #include "Ifpack_LocalFilter.h"
00014 #include "Ifpack_SingletonFilter.h"
00015 #include "Ifpack_ReorderFilter.h"
00016 #include "Ifpack_Utils.h"
00017 #include "Ifpack_OverlappingRowMatrix.h"
00018 #include "Epetra_CombineMode.h"
00019 #include "Epetra_MultiVector.h"
00020 #include "Epetra_Map.h"
00021 #include "Epetra_Comm.h"
00022 #include "Epetra_Time.h"
00023 #include "Epetra_LinearProblem.h"
00024 #include "Epetra_RowMatrix.h"
00025 #include "Epetra_CrsMatrix.h"
00026 #include "Teuchos_ParameterList.hpp"
00027 
00029 
00082 template<typename T>
00083 class Ifpack_AdditiveSchwarz : public virtual Ifpack_Preconditioner {
00084       
00085 public:
00086 
00088 
00089 
00100   Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix,
00101              int OverlapLevel = 0);
00102   
00104   virtual ~Ifpack_AdditiveSchwarz();
00106 
00108 
00110 
00119     virtual int SetUseTranspose(bool UseTranspose);
00121   
00123 
00125 
00135     virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00136 
00138 
00149     virtual int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00150 
00152     virtual double NormInf() const;
00154   
00156 
00158     virtual const char * Label() const;
00159 
00161     virtual bool UseTranspose() const;
00162 
00164     virtual bool HasNormInf() const;
00165 
00167     virtual const Epetra_Comm & Comm() const;
00168 
00170     virtual const Epetra_Map & OperatorDomainMap() const;
00171 
00173     virtual const Epetra_Map & OperatorRangeMap() const;
00175 
00177   virtual bool IsInitialized() const
00178   {
00179     return(IsInitialized_);
00180   }
00181 
00183   virtual bool IsComputed() const
00184   {
00185     return(IsComputed_);
00186   }
00187 
00189 
00208   virtual int SetParameters(Teuchos::ParameterList& List);
00209 
00210   // @}
00211 
00212   // @{ Query methods
00213   
00215   virtual int Initialize();
00216 
00218   virtual int Compute();
00219 
00221   virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
00222                          const int MaxIters = 1550,
00223                          const double Tol = 1e-9,
00224              Epetra_RowMatrix* Matrix = 0);
00225 
00227   virtual double Condest() const
00228   {
00229     return(Condest_);
00230   }
00231 
00233   virtual const Epetra_RowMatrix& Matrix() const
00234   {
00235     return(*Matrix_);
00236   }
00237 
00239   virtual bool IsOverlapping() const
00240   {
00241     return(IsOverlapping_);
00242   }
00243 
00245   virtual std::ostream& Print(std::ostream&) const;
00246   
00247   virtual const T* Inverse() const
00248   {
00249     return(Inverse_);
00250   }
00251 
00253   virtual int NumInitialize() const
00254   {
00255     return(NumInitialize_);
00256   }
00257 
00259   virtual int NumCompute() const
00260   {
00261     return(NumCompute_);
00262   }
00263 
00265   virtual int NumApplyInverse() const
00266   {
00267     return(NumApplyInverse_);
00268   }
00269 
00271   virtual double InitializeTime() const
00272   {
00273     return(InitializeTime_);
00274   }
00275 
00277   virtual double ComputeTime() const
00278   {
00279     return(ComputeTime_);
00280   }
00281 
00283   virtual double ApplyInverseTime() const
00284   {
00285     return(ApplyInverseTime_);
00286   }
00287 
00289   virtual double InitializeFlops() const
00290   {
00291     return(InitializeFlops_);
00292   }
00293 
00294   virtual double ComputeFlops() const
00295   {
00296     return(ComputeFlops_);
00297   }
00298 
00299   virtual double ApplyInverseFlops() const
00300   {
00301     return(ApplyInverseFlops_);
00302   }
00303 
00305   virtual int OverlapLevel() const
00306   {
00307     return(OverlapLevel_);
00308   }
00309 
00311   virtual const Teuchos::ParameterList& List() const
00312   {
00313     return(List_);
00314   }
00315 
00316 protected:
00317 
00318   // @}
00319 
00320   // @{ Internal merhods.
00321   
00323   Ifpack_AdditiveSchwarz(const Ifpack_AdditiveSchwarz& RHS)
00324   { }
00325 
00327   int Setup();
00328   
00330   void Destroy();
00331 
00332   // @}
00333 
00334   // @{ Internal data.
00335   
00337   const Epetra_RowMatrix* Matrix_;
00339   Ifpack_OverlappingRowMatrix* OverlappingMatrix_;
00341   Ifpack_LocalFilter* LocalizedMatrix_;
00343   string Label_;
00345   bool IsInitialized_;
00347   bool IsComputed_;
00349   T* Inverse_;
00351   bool UseTranspose_;
00353   bool IsOverlapping_;
00355   int OverlapLevel_;
00357   Teuchos::ParameterList List_;
00359   Epetra_CombineMode CombineMode_;
00361   double Condest_;
00363   bool ComputeCondest_;
00365   bool UseReordering_;
00367   string ReorderingType_;
00369   Ifpack_Reordering* Reordering_;
00371   Ifpack_ReorderFilter* ReorderedLocalizedMatrix_;
00373   bool FilterSingletons_;
00375   Ifpack_SingletonFilter* SingletonFilter_;
00377   int NumInitialize_;
00379   int NumCompute_;
00381   mutable int NumApplyInverse_;
00383   double InitializeTime_;
00385   double ComputeTime_;
00387   mutable double ApplyInverseTime_;
00389   double InitializeFlops_;
00391   double ComputeFlops_;
00393   mutable double ApplyInverseFlops_;
00395   Epetra_Time* Time_;
00396 
00397 }; // class Ifpack_AdditiveSchwarz<T>
00398 
00399 //==============================================================================
00400 template<typename T>
00401 Ifpack_AdditiveSchwarz<T>::
00402 Ifpack_AdditiveSchwarz(Epetra_RowMatrix* Matrix,
00403                int OverlapLevel) :
00404   Matrix_(Matrix),
00405   OverlappingMatrix_(0),
00406   LocalizedMatrix_(0),
00407   IsInitialized_(false),
00408   IsComputed_(false),
00409   Inverse_(0),
00410   UseTranspose_(false),
00411   IsOverlapping_(false),
00412   OverlapLevel_(OverlapLevel),
00413   CombineMode_(Zero),
00414   Condest_(-1.0),
00415   ComputeCondest_(true),
00416   UseReordering_(false),
00417   ReorderingType_("none"),
00418   Reordering_(0),
00419   ReorderedLocalizedMatrix_(0),
00420   FilterSingletons_(false),
00421   SingletonFilter_(0),
00422   NumInitialize_(0),
00423   NumCompute_(0),
00424   NumApplyInverse_(0),
00425   InitializeTime_(0.0),
00426   ComputeTime_(0.0),
00427   ApplyInverseTime_(0.0),
00428   InitializeFlops_(0.0),
00429   ComputeFlops_(0.0),
00430   ApplyInverseFlops_(0.0),
00431   Time_(0)
00432 {
00433   if (Matrix_->Comm().NumProc() == 1)
00434     OverlapLevel_ = 0;
00435 
00436   if ((OverlapLevel_ != 0) && (Matrix_->Comm().NumProc() > 1))
00437     IsOverlapping_ = true;
00438   // Sets parameters to default values
00439   Teuchos::ParameterList List;
00440   SetParameters(List);
00441 }
00442 
00443 //==============================================================================
00444 template<typename T>
00445 Ifpack_AdditiveSchwarz<T>::~Ifpack_AdditiveSchwarz()
00446 {
00447   Destroy();
00448 }
00449 
00450 //==============================================================================
00451 template<typename T>
00452 void Ifpack_AdditiveSchwarz<T>::Destroy() 
00453 {
00454   if (OverlappingMatrix_)
00455     delete OverlappingMatrix_;
00456   OverlappingMatrix_ = 0;
00457 
00458   if (Inverse_)
00459     delete Inverse_;
00460   Inverse_ = 0;
00461 
00462   if (LocalizedMatrix_)
00463     delete LocalizedMatrix_;
00464   LocalizedMatrix_ = 0;
00465 
00466   if (ReorderedLocalizedMatrix_)
00467     delete ReorderedLocalizedMatrix_;
00468   ReorderedLocalizedMatrix_ = 0;
00469 
00470   if (SingletonFilter_)
00471     delete SingletonFilter_;
00472   SingletonFilter_ = 0;
00473 
00474   if (Reordering_)
00475     delete Reordering_;
00476   Reordering_ = 0;
00477 
00478   if (Time_)
00479     delete Time_;
00480   Time_ = 0;
00481 }
00482 
00483 //==============================================================================
00484 template<typename T>
00485 int Ifpack_AdditiveSchwarz<T>::Setup()
00486 {
00487 
00488   Epetra_RowMatrix* MatrixPtr;
00489 
00490   if (OverlappingMatrix_)
00491     LocalizedMatrix_ = new Ifpack_LocalFilter(OverlappingMatrix_);
00492   else
00493     LocalizedMatrix_ = new Ifpack_LocalFilter(Matrix_);
00494 
00495   if (LocalizedMatrix_ == 0)
00496     IFPACK_CHK_ERR(-5);
00497 
00498   // users may want to skip singleton check
00499   if (FilterSingletons_) {
00500     SingletonFilter_ = new Ifpack_SingletonFilter(LocalizedMatrix_);
00501     MatrixPtr = SingletonFilter_;
00502   }
00503   else
00504     MatrixPtr = LocalizedMatrix_;
00505 
00506   if (UseReordering_) {
00507 
00508     // create reordeing and compute it
00509     if (ReorderingType_ == "rcm")
00510       Reordering_ = new Ifpack_RCMReordering();
00511     else if (ReorderingType_ == "metis")
00512       Reordering_ = new Ifpack_METISReordering();
00513     else {
00514       cerr << "reordering type not correct (" << ReorderingType_ << ")" << endl;
00515       exit(EXIT_FAILURE);
00516     }
00517     if (Reordering_ == 0) IFPACK_CHK_ERR(-5);
00518 
00519     IFPACK_CHK_ERR(Reordering_->SetParameters(List_));
00520     IFPACK_CHK_ERR(Reordering_->Compute(*MatrixPtr));
00521 
00522     // now create reordered localized matrix
00523     ReorderedLocalizedMatrix_ = 
00524       new Ifpack_ReorderFilter(MatrixPtr,Reordering_);
00525 
00526     if (ReorderedLocalizedMatrix_ == 0) IFPACK_CHK_ERR(-5);
00527 
00528     MatrixPtr = ReorderedLocalizedMatrix_;
00529   }
00530 
00531   Inverse_ = new T(MatrixPtr);
00532 
00533   if (Inverse_ == 0)
00534     IFPACK_CHK_ERR(-5);
00535 
00536   return(0);
00537 }
00538 
00539 //==============================================================================
00540 template<typename T>
00541 int Ifpack_AdditiveSchwarz<T>::SetParameters(Teuchos::ParameterList& List)
00542 {
00543  
00544   // compute the condition number each time Compute() is invoked.
00545   ComputeCondest_ = List.get("schwarz: compute condest", ComputeCondest_);
00546   // combine mode
00547   if( Teuchos::ParameterEntry *combineModeEntry = List.getEntryPtr("schwarz: combine mode") )
00548   {
00549     if( typeid(std::string) == combineModeEntry->getAny().type() )
00550     {
00551       std::string mode = List.get("schwarz: combine mode", "Add");
00552       if (mode == "Add")
00553         CombineMode_ = Add;
00554       else if (mode == "Zero")
00555         CombineMode_ = Zero;
00556       else if (mode == "Insert")
00557         CombineMode_ = Insert;
00558       else if (mode == "InsertAdd")
00559         CombineMode_ = InsertAdd;
00560       else if (mode == "Average")
00561         CombineMode_ = Average;
00562       else if (mode == "AbsMax")
00563         CombineMode_ = AbsMax;
00564       else
00565       {
00566         TEST_FOR_EXCEPTION(
00567           true,std::logic_error
00568           ,"Error, The (Epetra) combine mode of \""<<mode<<"\" is not valid!  Only the values"
00569           " \"Add\", \"Zero\", \"Insert\", \"InsertAdd\", \"Average\", and \"AbsMax\" are accepted!"
00570           );
00571       }
00572     }
00573     else if ( typeid(Epetra_CombineMode) == combineModeEntry->getAny().type() )
00574     {
00575       CombineMode_ = Teuchos::any_cast<Epetra_CombineMode>(combineModeEntry->getAny());
00576     }
00577     else
00578     {
00579       // Throw exception with good error message!
00580       Teuchos::getParameter<std::string>(List,"schwarz: combine mode");
00581     }
00582   }
00583   else
00584   {
00585     // Make the default be a string to be consistent with the valid parameters!
00586     List.get("schwarz: combine mode","Zero");
00587   }
00588   // type of reordering
00589   ReorderingType_ = List.get("schwarz: reordering type", ReorderingType_);
00590   if (ReorderingType_ == "none")
00591     UseReordering_ = false;
00592   else 
00593     UseReordering_ = true;
00594   // if true, filter singletons. NOTE: the filtered matrix can still have
00595   // singletons! A simple example: upper triangular matrix, if I remove
00596   // the lower node, I still get a matrix with a singleton! However, filter
00597   // singletons should help for PDE problems with Dirichlet BCs.
00598   FilterSingletons_ = List.get("schwarz: filter singletons", FilterSingletons_);
00599 
00600   // This copy may be needed by Amesos or other preconditioners.
00601   List_ = List;
00602 
00603   return(0);
00604 }
00605 
00606 //==============================================================================
00607 template<typename T>
00608 int Ifpack_AdditiveSchwarz<T>::Initialize()
00609 {
00610   IsInitialized_ = false;
00611   IsComputed_ = false; // values required
00612   Condest_ = -1.0; // zero-out condest
00613 
00614   Destroy();
00615 
00616   if (Time_ == 0)
00617     Time_ = new Epetra_Time(Comm());
00618 
00619   Time_->ResetStartTime();
00620 
00621   // compute the overlapping matrix if necessary
00622   if (IsOverlapping_) {
00623     OverlappingMatrix_ = 
00624       new Ifpack_OverlappingRowMatrix(Matrix_, OverlapLevel_);
00625     if (OverlappingMatrix_ == 0)
00626       IFPACK_CHK_ERR(-5);
00627   }
00628 
00629   IFPACK_CHK_ERR(Setup());
00630 
00631   if (Inverse_ == 0)
00632     IFPACK_CHK_ERR(-5);
00633 
00634   if (LocalizedMatrix_ == 0)
00635     IFPACK_CHK_ERR(-5);
00636 
00637   IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose()));
00638   IFPACK_CHK_ERR(Inverse_->SetParameters(List_));
00639   IFPACK_CHK_ERR(Inverse_->Initialize());
00640 
00641   // Label is for Aztec-like solvers
00642   Label_ = "Ifpack_AdditiveSchwarz, ";
00643   if (UseTranspose())
00644     Label_ += ", transp";
00645   Label_ += ", ov = " + Ifpack_toString(OverlapLevel_)
00646     + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'";
00647 
00648   IsInitialized_ = true;
00649   ++NumInitialize_;
00650   InitializeTime_ += Time_->ElapsedTime();
00651 
00652   // count flops by summing up all the InitializeFlops() in each
00653   // Inverse. Each Inverse() can only give its flops -- it acts on one
00654   // process only
00655   double partial = Inverse_->InitializeFlops();
00656   double total;
00657   Comm().SumAll(&partial, &total, 1);
00658   InitializeFlops_ += total;
00659 
00660   return(0);
00661 }
00662 
00663 //==============================================================================
00664 template<typename T>
00665 int Ifpack_AdditiveSchwarz<T>::Compute()
00666 {
00667 
00668   if (IsInitialized() == false)
00669     IFPACK_CHK_ERR(Initialize());
00670 
00671   Time_->ResetStartTime();
00672   IsComputed_ = false;
00673   Condest_ = -1.0;
00674   
00675   IFPACK_CHK_ERR(Inverse_->Compute());
00676 
00677   IsComputed_ = true; // need this here for Condest(Ifpack_Cheap)
00678   ++NumCompute_;
00679   ComputeTime_ += Time_->ElapsedTime();
00680 
00681   // sum up flops
00682   double partial = Inverse_->ComputeFlops();
00683   double total;
00684   Comm().SumAll(&partial, &total, 1);
00685   ComputeFlops_ += total;
00686 
00687   // reset the Label
00688   string R = "";
00689   if (UseReordering_)
00690     R = ReorderingType_ + " reord, ";
00691 
00692   if (ComputeCondest_)
00693     Condest(Ifpack_Cheap);
00694   
00695   // add Condest() to label
00696   Label_ = "Ifpack_AdditiveSchwarz, ov = " + Ifpack_toString(OverlapLevel_)
00697     + ", local solver = \n\t\t***** `" + string(Inverse_->Label()) + "'"
00698     + "\n\t\t***** " + R + "Condition number estimate = "
00699     + Ifpack_toString(Condest());
00700 
00701   return(0);
00702 }
00703 
00704 //==============================================================================
00705 template<typename T>
00706 int Ifpack_AdditiveSchwarz<T>::SetUseTranspose(bool UseTranspose)
00707 {
00708   // store the flag -- it will be set in Initialize() if Inverse_ does not
00709   // exist.
00710   UseTranspose_ = UseTranspose;
00711 
00712   // If Inverse_ exists, pass it right now.
00713   if (Inverse_)
00714     IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose));
00715   return(0);
00716 }
00717 
00718 //==============================================================================
00719 template<typename T>
00720 int Ifpack_AdditiveSchwarz<T>::
00721 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00722 {
00723   IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
00724   return(0);
00725 }
00726 
00727 //==============================================================================
00728 template<typename T>
00729 double Ifpack_AdditiveSchwarz<T>::NormInf() const
00730 {
00731   return(-1.0);
00732 }
00733 
00734 //==============================================================================
00735 template<typename T>
00736 const char * Ifpack_AdditiveSchwarz<T>::Label() const
00737 {
00738   return(Label_.c_str());
00739 }
00740 
00741 //==============================================================================
00742 template<typename T>
00743 bool Ifpack_AdditiveSchwarz<T>::UseTranspose() const
00744 {
00745   return(UseTranspose_);
00746 }
00747 
00748 //==============================================================================
00749 template<typename T>
00750 bool Ifpack_AdditiveSchwarz<T>::HasNormInf() const
00751 {
00752   return(false);
00753 }
00754 
00755 //==============================================================================
00756 template<typename T>
00757 const Epetra_Comm & Ifpack_AdditiveSchwarz<T>::Comm() const
00758 {
00759   return(Matrix_->Comm());
00760 }
00761 
00762 //==============================================================================
00763 template<typename T>
00764 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorDomainMap() const
00765 {
00766   return(Matrix_->OperatorDomainMap());
00767 }
00768 
00769 //==============================================================================
00770 template<typename T>
00771 const Epetra_Map & Ifpack_AdditiveSchwarz<T>::OperatorRangeMap() const
00772 {
00773   return(Matrix_->OperatorRangeMap());
00774 }
00775 
00776 //==============================================================================
00777 template<typename T>
00778 int Ifpack_AdditiveSchwarz<T>::
00779 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00780 {
00781   // compute the preconditioner is not done by the user
00782   if (!IsComputed())
00783     IFPACK_CHK_ERR(-3);
00784   
00785   int NumVectors = X.NumVectors();
00786 
00787   if (NumVectors != Y.NumVectors())
00788     IFPACK_CHK_ERR(-2); // wrong input
00789 
00790   Time_->ResetStartTime();
00791 
00792   Epetra_MultiVector* OverlappingX;
00793   Epetra_MultiVector* OverlappingY;
00794   Epetra_MultiVector* Xtmp = 0;
00795 
00796   // for flop count, see bottom of this function
00797   double pre_partial = Inverse_->ApplyInverseFlops();
00798   double pre_total;
00799   Comm().SumAll(&pre_partial, &pre_total, 1);
00800 
00801   // process overlap, may need to create vectors and import data
00802   if (IsOverlapping()) {
00803     OverlappingX = new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00804                       X.NumVectors());
00805     OverlappingY = new Epetra_MultiVector(OverlappingMatrix_->RowMatrixRowMap(),
00806                                           Y.NumVectors());
00807     if (OverlappingY == 0) IFPACK_CHK_ERR(-5);
00808 
00809     OverlappingY->PutScalar(0.0);
00810     OverlappingX->PutScalar(0.0);
00811     IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(X,*OverlappingX,Insert));
00812     // FIXME: this will not work with overlapping and non-zero starting
00813     // solutions. The same for other cases below.
00814     // IFPACK_CHK_ERR(OverlappingMatrix_->ImportMultiVector(Y,*OverlappingY,Insert));
00815   }
00816   else {
00817     Xtmp = new Epetra_MultiVector(X);
00818     OverlappingX = (Epetra_MultiVector*)Xtmp;
00819     OverlappingY = &Y;
00820   }
00821 
00822   if (FilterSingletons_) {
00823     // process singleton filter
00824     Epetra_MultiVector ReducedX(SingletonFilter_->Map(),NumVectors);
00825     Epetra_MultiVector ReducedY(SingletonFilter_->Map(),NumVectors);
00826     IFPACK_CHK_ERR(SingletonFilter_->SolveSingletons(*OverlappingX,*OverlappingY));
00827     IFPACK_CHK_ERR(SingletonFilter_->CreateReducedRHS(*OverlappingY,*OverlappingX,ReducedX));
00828 
00829     // process reordering
00830     if (!UseReordering_) {
00831       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReducedX,ReducedY));
00832     }
00833     else {
00834       Epetra_MultiVector ReorderedX(ReducedX);
00835       Epetra_MultiVector ReorderedY(ReducedY);
00836       IFPACK_CHK_ERR(Reordering_->P(ReducedX,ReorderedX));
00837       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
00838       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,ReducedY));
00839     }
00840 
00841     // finish up with singletons
00842     IFPACK_CHK_ERR(SingletonFilter_->UpdateLHS(ReducedY,*OverlappingY));
00843   }
00844   else {
00845     // process reordering
00846     if (!UseReordering_) {
00847       IFPACK_CHK_ERR(Inverse_->ApplyInverse(*OverlappingX,*OverlappingY));
00848     }
00849     else {
00850       Epetra_MultiVector ReorderedX(*OverlappingX);
00851       Epetra_MultiVector ReorderedY(*OverlappingY);
00852       IFPACK_CHK_ERR(Reordering_->P(*OverlappingX,ReorderedX));
00853       IFPACK_CHK_ERR(Inverse_->ApplyInverse(ReorderedX,ReorderedY));
00854       IFPACK_CHK_ERR(Reordering_->Pinv(ReorderedY,*OverlappingY));
00855     }
00856   }
00857 
00858   if (IsOverlapping()) {
00859     IFPACK_CHK_ERR(OverlappingMatrix_->ExportMultiVector(*OverlappingY,Y,
00860                                                         CombineMode_));
00861 
00862     delete OverlappingX;
00863     delete OverlappingY;
00864   }
00865 
00866   if (Xtmp)
00867     delete Xtmp;
00868 
00869   // add flops. Note the we only have to add the newly counted
00870   // flops -- and each Inverse returns the cumulative sum
00871   double partial = Inverse_->ApplyInverseFlops();
00872   double total;
00873   Comm().SumAll(&partial, &total, 1);
00874   ApplyInverseFlops_ += total - pre_total;
00875 
00876   // FIXME: right now I am skipping the overlap and singletons
00877   ++NumApplyInverse_;
00878   ApplyInverseTime_ += Time_->ElapsedTime();
00879 
00880   return(0);
00881  
00882 }
00883 
00884 //==============================================================================
00885 template<typename T>
00886 std::ostream& Ifpack_AdditiveSchwarz<T>::
00887 Print(std::ostream& os) const
00888 {
00889   double IF = InitializeFlops();
00890   double CF = ComputeFlops();
00891   double AF = ApplyInverseFlops();
00892 
00893   double IFT = 0.0, CFT = 0.0, AFT = 0.0;
00894   if (InitializeTime() != 0.0)
00895     IFT = IF / InitializeTime();
00896   if (ComputeTime() != 0.0)
00897     CFT = CF / ComputeTime();
00898   if (ApplyInverseTime() != 0.0)
00899     AFT = AF / ApplyInverseTime();
00900 
00901   if (Matrix().Comm().MyPID())
00902     return(os);
00903 
00904   os << endl;
00905   os << "================================================================================" << endl;
00906   os << "Ifpack_AdditiveSchwarz, overlap level = " << OverlapLevel_ << endl;
00907   if (CombineMode_ == Insert)
00908     os << "Combine mode                          = Insert" << endl;
00909   else if (CombineMode_ == Add)
00910     os << "Combine mode                          = Add" << endl;
00911   else if (CombineMode_ == Zero)
00912     os << "Combine mode                          = Zero" << endl;
00913   else if (CombineMode_ == Average)
00914     os << "Combine mode                          = Average" << endl;
00915   else if (CombineMode_ == AbsMax)
00916     os << "Combine mode                          = AbsMax" << endl;
00917 
00918   os << "Condition number estimate             = " << Condest_ << endl;
00919   os << "Global number of rows                 = " << Matrix_->NumGlobalRows() << endl;
00920 
00921   os << endl;
00922   os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00923   os << "-----           -------   --------------       ------------     --------" << endl;
00924   os << "Initialize()    "   << std::setw(5) << NumInitialize()
00925      << "  " << std::setw(15) << InitializeTime() 
00926      << "  " << std::setw(15) << 1.0e-6 * IF 
00927      << "  " << std::setw(15) << 1.0e-6 * IFT << endl;
00928   os << "Compute()       "   << std::setw(5) << NumCompute() 
00929      << "  " << std::setw(15) << ComputeTime()
00930      << "  " << std::setw(15) << 1.0e-6 * CF
00931      << "  " << std::setw(15) << 1.0e-6 * CFT << endl;
00932   os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00933      << "  " << std::setw(15) << ApplyInverseTime()
00934      << "  " << std::setw(15) << 1.0e-6 * AF
00935      << "  " << std::setw(15) << 1.0e-6 * AFT << endl;
00936   os << "================================================================================" << endl;
00937   os << endl;
00938 
00939   return(os);
00940 }
00941 
00942 #include "Ifpack_Condest.h"
00943 //==============================================================================
00944 template<typename T>
00945 double Ifpack_AdditiveSchwarz<T>::
00946 Condest(const Ifpack_CondestType CT, const int MaxIters, 
00947         const double Tol, Epetra_RowMatrix* Matrix)
00948 {
00949   if (!IsComputed()) // cannot compute right now
00950     return(-1.0);
00951 
00952   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);
00953 
00954   return(Condest_);
00955 }
00956 
00957 #endif // HAVE_IFPACK_TEUCHOS
00958 #endif // IFPACK_ADDITIVESCHWARZ_H

Generated on Thu Sep 18 12:37:07 2008 for IFPACK by doxygen 1.3.9.1