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

MLAPI_MultiLevelSA.h

Go to the documentation of this file.
00001 #ifndef MLAPI_MULTILEVEL_H
00002 #define MLAPI_MULTILEVEL_H
00003 
00013 /* ******************************************************************** */
00014 /* See the file COPYRIGHT for a complete copyright notice, contact      */
00015 /* person and disclaimer.                                               */        
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   // @{ \name Constructors and destructors
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   // @{ \name Set and Get methods
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   // @{ \name Mathematical methods
00132 
00134   void Compute() 
00135   {
00136     ResetTimer();
00137     StackPush();
00138     IsComputed_ = false;
00139 
00140     // get parameter from the input list
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     // build up the default null space
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;     // contains the next-level null space
00167 
00168     A_.resize(MaxLevels);
00169     R_.resize(MaxLevels);
00170     P_.resize(MaxLevels);
00171     S_.resize(MaxLevels);
00172 
00173     // work on increasing hierarchies only.
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       // only an alias
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       // load current level into database
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       // build smoothers
00254       S.Reshape(A, SmootherType, List_);
00255 
00256       // put operators and inverse in hierarchy
00257       R_[level    ] = R;
00258       P_[level    ] = P;
00259       A_[level + 1] = C;
00260       S_[level    ] = S;
00261 
00262       // break if coarse matrix is below specified tolerance
00263       if (C.GetNumGlobalRows() <= MaxCoarseSize) {
00264         ++level;
00265         break;
00266       }
00267     }
00268 
00269     // set coarse solver
00270     S.Reshape(A_[level], CoarseType, List_);
00271     S_[level] = S;
00272     MaxLevels_ = level + 1;
00273 
00274     // set the label
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     // FIXME: update flops!
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     // reset flop counter
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     // apply pre-smoother
00330     x_f = S(level) * b_f;
00331     // new residual
00332     r_f = b_f - A(level) * x_f;
00333     // restrict to coarse
00334     r_c = R(level) * r_f;
00335     // solve coarse problem
00336     SolveMultiLevelSA(r_c,z_c,level + 1);
00337     // prolongate back and add to solution
00338     x_f = x_f + P(level) * z_c;
00339     // apply post-smoother
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   // @{ \name Miscellaneous methods
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 } // namespace MLAPI
00398 
00399 #endif