AFEPack
AMGSolver.h
浏览该文件的文档。
00001 
00011 #ifndef _AMGSolver_h_
00012 #define _AMGSolver_h_
00013 
00014 #include <iostream>
00015 #include <fstream>
00016 #include <cmath>
00017 #include <vector>
00018 #include <list>
00019 #include <base/exceptions.h>
00020 #include <lac/vector.h>
00021 #include <lac/full_matrix.h>
00022 #include <lac/sparsity_pattern.h>
00023 #include <lac/sparse_matrix.h>
00024 
00025 #include "Miscellaneous.h"
00026 
00028 
00064 class AMGSolver {
00065  public:
00066   typedef SparseMatrix<double> Matrix;
00067  private:
00068   bool is_initialized; 
00069   u_int n_project; 
00070   std::vector<Matrix *> project_matrix; 
00071   std::vector<Matrix *> project_matrix_r; 
00072   std::vector<const Matrix *> projected_matrix; 
00073   bool is_most_project_full; 
00074   bool is_solve_most_project_exactly; 
00075   FullMatrix<double> M_n; 
00076   double toler; 
00077   double alpha; 
00078   u_int smooth_step; 
00079   u_int n_most_project_order; 
00080   double n_most_project_sparse_degree; 
00081  public:
00082   AMGSolver(); 
00083   AMGSolver(const Matrix&, 
00084             double tol = 1.0e-12, 
00085             u_int s = 3, 
00086             u_int nmpo = 50,
00087             double nmpsd = 0.382,
00088             double alp = 0.25);
00089   virtual ~AMGSolver(); 
00090   void clear(); 
00091 
00092   double tolerence() const {return toler;}; 
00093   double& tolerence() {return toler;};
00094   u_int smoothStep() const {return smooth_step;};
00095   u_int& smoothStep() {return smooth_step;};
00096   bool& isSolveMostProjectExactly() {return is_solve_most_project_exactly;}
00097   bool isSolveMostProjectExactly() const {return is_solve_most_project_exactly;}
00098 
00099   void reinit(const Matrix&); 
00100 
00110   void solve(Vector<double>& x, 
00111              const Vector<double>& r, 
00112              double tol = 0.0, 
00113              u_int step = 20,
00114              int mode = 0) const; 
00115 
00116   void lazyReinit(const Matrix&); 
00117 
00118  private:
00119   void lazyInit(const Matrix&);
00120   void lazyProject(const Matrix& M, 
00121                    Matrix *& P, 
00122                    Matrix *& PMPt,
00123                    Matrix *& Pt);
00124 
00125   void Project(const Matrix& M, 
00126                Matrix *& P, 
00127                Matrix *& PMPt);
00128 
00129   void GaussSidel(const Matrix& M, 
00130                   Vector<double>& x, 
00131                   const Vector<double>& r, 
00132                   const int& s) const;
00133   Matrix * getPMPt(const Matrix & P, 
00134                    const Matrix & M,
00135                    const Matrix & Pt) const;
00136 };
00137 
00138 
00139 
00141 
00167 class AMGPreconditioner : public AMGSolver
00168 {
00169  private:
00170   u_int iterate_step; 
00171  public:
00172   AMGPreconditioner() {};
00173   AMGPreconditioner(const Matrix& M, 
00174                     u_int is = 5, 
00175                     u_int s = 3) : 
00176     iterate_step(is) {
00177     lazyReinit(M);
00178     smoothStep() = s;
00179   }
00180   ~AMGPreconditioner() {}
00181  public:
00182   void reinit(const Matrix& M,
00183               u_int is = 5) {
00184     AMGSolver::lazyReinit(M);
00185     iterate_step = is;
00186   }
00187   u_int& iterateStep() {return iterate_step;}
00188   const u_int& iterateStep() const {return iterate_step;}
00189   void vmult(Vector<double>& x, 
00190              const Vector<double>& b) const
00191     {solve(x, b, 0, iterate_step, 1);}
00192 };
00193 
00194 
00195 
00196 #endif
00197 
00198 //
00199 // end of file
00201