AFEPack
|
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