00001 #ifndef MLAPI_MULTILEVEL_H
00002 #define MLAPI_MULTILEVEL_H
00003
00013
00014
00015
00016
00017
00018 #include "ml_common.h"
00019 #include "ml_agg_genP.h"
00020 #include "MLAPI_Error.h"
00021 #include "MLAPI_CompObject.h"
00022 #include "MLAPI_TimeObject.h"
00023 #include "MLAPI_Operator.h"
00024 #include "MLAPI_Operator_Utils.h"
00025 #include "MLAPI_MultiVector.h"
00026 #include "MLAPI_InverseOperator.h"
00027 #include "MLAPI_Expressions.h"
00028 #include "MLAPI_BaseOperator.h"
00029 #include "MLAPI_Workspace.h"
00030 #include "MLAPI_Aggregation.h"
00031 #include "MLAPI_Eig.h"
00032 #include <vector>
00033
00034 namespace MLAPI {
00035
00047 class MultiLevelSA : public BaseOperator, public CompObject, public TimeObject {
00048
00049 public:
00050
00051
00052
00054 MultiLevelSA(const Operator FineMatrix, Teuchos::ParameterList& List,
00055 const bool ConstructNow = true) :
00056 IsComputed_(false)
00057 {
00058 FineMatrix_ = FineMatrix;
00059 List_ = List;
00060 if (ConstructNow) Compute();
00061 }
00062
00064 virtual ~MultiLevelSA()
00065 { }
00066
00067
00068
00069
00071 const Space GetOperatorDomainSpace() const
00072 {
00073 return(FineMatrix_.GetDomainSpace());
00074 }
00075
00077 const Space GetOperatorRangeSpace() const
00078 {
00079 return(FineMatrix_.GetRangeSpace());
00080 }
00081
00083 inline const Space GetDomainSpace() const
00084 {
00085 return(FineMatrix_.GetDomainSpace());
00086 }
00087
00089 inline const Space GetRangeSpace() const
00090 {
00091 return(FineMatrix_.GetRangeSpace());
00092 }
00093
00095 inline const Operator& R(const int i) const
00096 {
00097 return(R_[i]);
00098 }
00099
00101 inline const Operator& A(const int i) const
00102 {
00103 return(A_[i]);
00104 }
00105
00107 inline const Operator& P(const int i) const
00108 {
00109 return(P_[i]);
00110 }
00111
00113 inline const InverseOperator& S(const int i) const
00114 {
00115 return(S_[i]);
00116 }
00117
00119 inline int GetMaxLevels() const
00120 {
00121 return(MaxLevels_);
00122 }
00123
00125 inline bool IsComputed() const
00126 {
00127 return(IsComputed_);
00128 }
00129
00130
00131
00132
00134 void Compute()
00135 {
00136 ResetTimer();
00137 StackPush();
00138 IsComputed_ = false;
00139
00140
00141 int MaxLevels = List_.get("max levels", 10);
00142 double Damping = List_.get("aggregation: damping factor", 1.3333);
00143 string EigenAnalysis = List_.get("eigen-analysis: type", "Anorm");
00144 int MaxCoarseSize = List_.get("coarse: max size", 32);
00145 MultiVector EmptySpace;
00146 MultiVector ThisNS = List_.get("aggregation: null space", EmptySpace);
00147 int NumPDEEqns = List_.get("PDE equations", 1);
00148 string SmootherType = List_.get("smoother: type", "symmetric Gauss-Seidel");
00149 string CoarseType = List_.get("coarse: type", "Amesos-KLU");
00150
00151
00152 if (ThisNS.GetNumVectors() == 0) {
00153 ThisNS.Reshape(FineMatrix_.GetDomainSpace(),NumPDEEqns);
00154 if (NumPDEEqns == 1)
00155 ThisNS = 1.0;
00156 else
00157 {
00158 ThisNS = 0.0;
00159 for (int i = 0 ; i < ThisNS.GetMyLength() ; ++i)
00160 for (int j = 0 ; j < NumPDEEqns ;++j)
00161 if (i % NumPDEEqns == j)
00162 ThisNS(i,j) = 1.0;
00163 }
00164 }
00165
00166 MultiVector NextNS;
00167
00168 A_.resize(MaxLevels);
00169 R_.resize(MaxLevels);
00170 P_.resize(MaxLevels);
00171 S_.resize(MaxLevels);
00172
00173
00174 A_[0] = FineMatrix_;
00175
00176 double LambdaMax;
00177 Operator A;
00178 Operator C;
00179 Operator R;
00180 Operator P;
00181 Operator Ptent;
00182 Operator IminusA;
00183 InverseOperator S;
00184
00185 int level;
00186
00187 for (level = 0 ; level < MaxLevels - 1 ; ++level) {
00188
00189
00190 A = A_[level];
00191
00192 if (level)
00193 List_.set("PDE equations", ThisNS.GetNumVectors());
00194
00195 if (GetPrintLevel()) {
00196 ML_print_line("-", 80);
00197 cout << "current working level = " << level << endl;
00198 cout << "number of global rows = " << A.GetNumGlobalRows() << endl;
00199 cout << "number of global nnz = " << A.GetNumGlobalNonzeros() << endl;
00200 cout << "threshold = " << List_.get("aggregation: threshold", 0.0) << endl;
00201 cout << "number of PDE equations = " << NumPDEEqns << endl;
00202 cout << "null space dimension = " << ThisNS.GetNumVectors() << endl;
00203 }
00204
00205
00206 List_.set("workspace: current level", level);
00207
00208 GetPtent(A, List_, ThisNS, Ptent, NextNS);
00209 ThisNS = NextNS;
00210
00211 if (Damping) {
00212
00213 if (EigenAnalysis == "Anorm")
00214 LambdaMax = MaxEigAnorm(A,true);
00215 else if (EigenAnalysis == "cg")
00216 LambdaMax = MaxEigCG(A,true);
00217 else if (EigenAnalysis == "power-method")
00218 LambdaMax = MaxEigPowerMethod(A,true);
00219 else
00220 ML_THROW("incorrect parameter (" + EigenAnalysis + ")", -1);
00221
00222 #if 0
00223 MultiVector Diag = GetDiagonal(A);
00224 Diag.Reciprocal();
00225 Diag.Scale(Damping / LambdaMax);
00226 Operator Dinv = GetDiagonal(Diag);
00227 Operator DinvA = Dinv * A;
00228 Operator I = GetIdentity(A.GetDomainSpace(),A.GetRangeSpace());
00229 Operator IminusA = I - DinvA;
00230 #else
00231 IminusA = GetJacobiIterationOperator(A,Damping / LambdaMax);
00232 #endif
00233 P = IminusA * Ptent;
00234 }
00235 else {
00236 P = Ptent;
00237 LambdaMax = -1.0;
00238 }
00239
00240 if (GetPrintLevel()) {
00241 cout << "omega = " << Damping << endl;
00242 if (LambdaMax != -1.0) {
00243 cout << "lambda max = " << LambdaMax << endl;
00244 cout << "damping factor = " << Damping / LambdaMax << endl;
00245 }
00246 cout << "smoother type = " << SmootherType << endl;
00247 cout << "relaxation sweeps = " << List_.get("smoother: sweeps", 1) << endl;
00248 cout << "smoother damping = " << List_.get("smoother: damping factor", 0.67) << endl;
00249 }
00250
00251 R = GetTranspose(P);
00252 C = GetRAP(R,A,P);
00253
00254 S.Reshape(A, SmootherType, List_);
00255
00256
00257 R_[level ] = R;
00258 P_[level ] = P;
00259 A_[level + 1] = C;
00260 S_[level ] = S;
00261
00262
00263 if (C.GetNumGlobalRows() <= MaxCoarseSize) {
00264 ++level;
00265 break;
00266 }
00267 }
00268
00269
00270 S.Reshape(A_[level], CoarseType, List_);
00271 S_[level] = S;
00272 MaxLevels_ = level + 1;
00273
00274
00275 SetLabel("SA, L = " + GetString(MaxLevels_) +
00276 ", smoother = " + SmootherType);
00277
00278 if (GetPrintLevel()) {
00279 ML_print_line("-", 80);
00280 cout << "final level = " << level << endl;
00281 cout << "number of global rows = " << A_[level].GetNumGlobalRows() << endl;
00282 cout << "number of global nnz = " << A_[level].GetNumGlobalNonzeros() << endl;
00283 cout << "coarse solver = " << CoarseType << endl;
00284 cout << "time spent in constr. = " << GetTime() << " (s)" << endl;
00285 ML_print_line("-", 80);
00286 }
00287
00288 IsComputed_ = true;
00289 StackPop();
00290
00291
00292 UpdateTime();
00293
00294 }
00295
00297 int Apply(const MultiVector& b_f, MultiVector& x_f) const
00298 {
00299 ResetTimer();
00300 StackPush();
00301
00302 if (IsComputed() == false)
00303 ML_THROW("Method Compute() must be called before Apply()", -1);
00304 SolveMultiLevelSA(b_f,x_f,0);
00305 UpdateTime();
00306
00307 StackPop();
00308 return(0);
00309 }
00310
00312 int SolveMultiLevelSA(const MultiVector& b_f,MultiVector& x_f, int level) const
00313 {
00314 if (level == MaxLevels_ - 1) {
00315 x_f = S(level) * b_f;
00316 return(0);
00317 }
00318
00319 MultiVector r_f(P(level).GetRangeSpace());
00320 MultiVector r_c(P(level).GetDomainSpace());
00321 MultiVector z_c(P(level).GetDomainSpace());
00322
00323
00324 S(level).SetFlops(0.0);
00325 A(level).SetFlops(0.0);
00326 R(level).SetFlops(0.0);
00327 P(level).SetFlops(0.0);
00328
00329
00330 x_f = S(level) * b_f;
00331
00332 r_f = b_f - A(level) * x_f;
00333
00334 r_c = R(level) * r_f;
00335
00336 SolveMultiLevelSA(r_c,z_c,level + 1);
00337
00338 x_f = x_f + P(level) * z_c;
00339
00340 S(level).Apply(b_f,x_f);
00341
00342 UpdateFlops(2.0 * S(level).GetFlops());
00343 UpdateFlops(A(level).GetFlops());
00344 UpdateFlops(R(level).GetFlops());
00345 UpdateFlops(P(level).GetFlops());
00346 UpdateFlops(2.0 * x_f.GetGlobalLength());
00347
00348 return(0);
00349 }
00350
00351
00352
00353
00355 std::ostream& Print(std::ostream& os,
00356 const bool verbose = true) const
00357 {
00358 if (GetMyPID() == 0) {
00359 os << endl;
00360 os << "*** MLAPI::MultiLevelSA, label = `" << GetLabel() << "'" << endl;
00361 os << endl;
00362 os << "Number of levels = " << GetMaxLevels() << endl;
00363 os << "Flop count = " << GetFlops() << endl;
00364 os << "Cumulative time = " << GetTime() << endl;
00365 if (GetTime() != 0.0)
00366 os << "MFlops rate = " << 1.0e-6 * GetFlops() / GetTime() << endl;
00367 else
00368 os << "MFlops rate = 0.0" << endl;
00369 os << endl;
00370 }
00371 return(os);
00372 }
00373
00374
00375
00376 private:
00377
00379 int MaxLevels_;
00381 Operator FineMatrix_;
00383 std::vector<Operator> A_;
00385 std::vector<Operator> R_;
00387 std::vector<Operator> P_;
00389 std::vector<InverseOperator> S_;
00391 Teuchos::ParameterList List_;
00393 bool IsComputed_;
00394
00395 };
00396
00397 }
00398
00399 #endif