• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

ml_MatrixFreePreconditioner.h

Go to the documentation of this file.
00001 /* ******************************************************************** */
00002 /* See the file COPYRIGHT for a complete copyright notice, contact      */
00003 /* person and disclaimer.                                               */        
00004 /* ******************************************************************** */
00005 #ifndef ML_MATRIX_FREE_PRECONDITIONER
00006 #define ML_MATRIX_FREE_PRECONDITIONER
00007 
00021 #include "ml_include.h"
00022 #if defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS) && defined(HAVE_ML_EPETRAEXT) && defined(HAVE_ML_IFPACK)
00023 #include <string>
00024 #include "ml_epetra.h"
00025 #include "Epetra_Time.h"
00026 #include "Epetra_Operator.h"
00027 #include "Epetra_Comm.h"
00028 #include "Epetra_CrsMatrix.h"
00029 #include "Teuchos_ParameterList.hpp"
00030 #include "Teuchos_RefCountPtr.hpp"
00031 #include <vector>
00032 #include <map>
00033 #include "ml_MultiLevelPreconditioner.h"
00034 
00035 class Epetra_Map;
00036 class Epetra_BlockMap;
00037 class Epetra_CrsGraph;
00038 class Epetra_Vector;
00039 class Epetra_MultiVector;
00040 class Epetra_RowMatrix;
00041 class Epetra_FECrsMatrix;
00042 class Ifpack_Chebyshev;
00043 
00044 namespace ML_Epetra {
00045 
00046 class MultiLevelPreconditioner;
00047 
00064 class MatrixFreePreconditioner : public Epetra_Operator 
00065 {
00066   public:
00068     
00070     MatrixFreePreconditioner(const Epetra_Operator& Operator,
00071                              const Epetra_CrsGraph& Graph,
00072                              Epetra_MultiVector& NullSpace,
00073                              const Epetra_Vector& PointDiagonal,
00074                              Teuchos::ParameterList& List);
00075 
00077     virtual ~MatrixFreePreconditioner();
00078 
00079     // @}
00080     // @{ \name Query methods.
00081     
00083     int SetUseTranspose(bool UseTranspose);
00084 
00086     int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00087 
00089     int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00090 
00092     double NormInf() const
00093     {
00094       return(-1.0);
00095     }
00096 
00098     const char * Label() const
00099     {
00100       return(Label_.c_str());
00101     }
00102 
00104     bool UseTranspose() const
00105     {
00106       return(false);
00107     }
00108 
00110     bool HasNormInf() const
00111     {
00112       return(false);
00113     }
00114 
00116     const Epetra_Comm& Comm() const
00117     {
00118       return(Comm_);
00119     }
00120 
00122     const Epetra_Map & OperatorDomainMap() const
00123     {
00124       return(Operator_.OperatorDomainMap());
00125     }
00126 
00128     const Epetra_Map & OperatorRangeMap() const
00129     {
00130       return(Operator_.OperatorRangeMap());
00131     }
00132 
00134     const Epetra_RowMatrix& C() const
00135     {
00136       if (!IsComputed())
00137         throw(-1); // prec not computed yet
00138 
00139       return(*C_);
00140     }
00141 
00142     const MultiLevelPreconditioner& MLP() const
00143     {
00144       if (!IsComputed())
00145         throw(-1); // prec not computed yet
00146 
00147       return(*MLP_);
00148     }
00149     
00151     const Epetra_CrsMatrix& R() const
00152     {
00153       return(*R_);
00154     }
00155 
00157     ML_Comm* Comm_ML()
00158     {
00159       return(Comm_ML_);
00160     }
00161 
00163     inline int MyPID() const
00164     {
00165       return(Comm().MyPID());
00166     }
00167 
00169     inline int NumProc() const
00170     {
00171       return(Comm().NumProc());
00172     }
00173 
00175     bool IsComputed() const
00176     {
00177       return(IsComputed_);
00178     }
00179 
00181     double TotalCPUTime() const;
00182 
00183     bool CheckSPD(const Epetra_Operator& Op, 
00184                   const bool UseApply = true,
00185                   const int NumChecks = 1,
00186                   const int NumVectors = 1) const;
00187 
00188     // @}
00189     // @{ \name Construction methods.
00190     
00192     int Coarsen(ML_Operator* A, ML_Aggregate** aggr, ML_Operator** P, 
00193                 ML_Operator** R, ML_Operator** C, int NumPDEEqns = 1, 
00194                 int NullSpaceDim = 1, double* NullSpace = NULL);
00195 
00197     int GetBlockDiagonal(const Epetra_CrsGraph& Graph, std::string DiagonalColoringType);
00198 
00199   private:
00200 
00202     int Compute(const Epetra_CrsGraph& Graph, Epetra_MultiVector& NullSpace);
00203 
00204     // @}
00205     // @{ \name Basic smoothers
00206     
00208     int ApplyPreSmoother(Epetra_MultiVector& X) const;
00209 
00211     int ApplyPostSmoother(Epetra_MultiVector& X, const Epetra_MultiVector& Y,
00212                           Epetra_MultiVector& tmp) const;
00213 
00215     int ApplyJacobi(Epetra_MultiVector& X, const double omega) const;
00216 
00218     int ApplyJacobi(Epetra_MultiVector& X, const Epetra_MultiVector& B,
00219                     const double omega, Epetra_MultiVector& tmp) const;
00220 
00222     int ApplyBlockJacobi(Epetra_MultiVector& X, const double omega) const;
00223 
00225     int ApplyBlockJacobi(Epetra_MultiVector& X, const Epetra_MultiVector& B,
00226                     const double omega, Epetra_MultiVector& tmp) const;
00227 
00228     int ApplyInvBlockDiag(const double alpha, Epetra_MultiVector& X,
00229                           const double gamma, const Epetra_MultiVector& B) const;
00230 
00231     // @} 
00232     // @{ \name Timing
00233     
00234     inline void ResetStartTime() const
00235     {
00236       Time_->ResetStartTime();
00237     }
00238 
00239     inline void AddAndResetStartTime(const std::string& Label, const int print = false) const
00240     {
00241       TimeTable[Label] += Time_->ElapsedTime();
00242       Time_->ResetStartTime();
00243       if (print)
00244       {
00245         if (MyPID() == 0 && ML_Get_PrintLevel() > 5)
00246           std::cout << "Time for " << Label << " = " << TimeTable[Label] << " (s)" << std::endl;
00247       }
00248     }
00249 
00250     void PrintTimings() const
00251     {
00252       if (MyPID() == 0)
00253       {
00254         double Total = 0.0;
00255         std::cout << "Cumulative timing so far:" << std::endl;
00256         std::cout << "- for coarsening = " << TimeTable["coarsening"] << std::endl;
00257         std::cout << "- total time     = " << Total << std::endl;
00258       }
00259     }
00260 
00261     // @} 
00262     // @{ \name Private data
00263     
00265     bool verbose_;
00267     ML_Comm* Comm_ML_;
00269     const Epetra_Comm& Comm_;
00270 
00272     std::string Label_; 
00274     bool IsComputed_;
00276     int PrecType_;
00278     int SmootherType_;
00280     double omega_;
00282     Teuchos::ParameterList List_;
00283 
00285     const Epetra_Operator& Operator_;
00287     Teuchos::RefCountPtr<Epetra_Vector> InvPointDiagonal_;
00289     std::vector<double> InvBlockDiag_;
00290 
00292     Teuchos::RefCountPtr<Ifpack_Chebyshev> PreSmoother_;
00294     Teuchos::RefCountPtr<Ifpack_Chebyshev> PostSmoother_;
00296     Teuchos::RefCountPtr<Epetra_CrsMatrix> R_;
00298     Teuchos::RefCountPtr<Epetra_RowMatrix> C_;
00300     ML_Operator* C_ML_;
00302     Teuchos::RefCountPtr<MultiLevelPreconditioner> MLP_;
00303 
00305     int NumPDEEqns_;
00306     int NumMyBlockRows_;
00307 
00309     mutable Teuchos::RefCountPtr<Epetra_Time> Time_;
00310     mutable std::map<std::string, double> TimeTable;
00311     // @} 
00312 
00313 }; // class MatrixFreePreconditioner
00314 
00315 } // namespace ML_Epetra
00316 
00317 #endif // HAVE_ML_EPETRA
00318 #endif // ML_MATRIX_FREE_PRECONDITIONER