00001
00002
00003
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
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);
00138
00139 return(*C_);
00140 }
00141
00142 const MultiLevelPreconditioner& MLP() const
00143 {
00144 if (!IsComputed())
00145 throw(-1);
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
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
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
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
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 };
00314
00315 }
00316
00317 #endif // HAVE_ML_EPETRA
00318 #endif // ML_MATRIX_FREE_PRECONDITIONER