#include <BiConjugateGradientStabilized.hpp>
Public Member Functions | |
| template<typename VectorType, typename MatrixType, typename PreconditionerType> | |
| BiConjugateGradientStabilized (const VectorType &b, const MatrixType &A, const PreconditionerType &P, VectorType &x) | |
Private Attributes | |
| real_t | epsilon |
| int | max_iter |
| GetParameter < BiConjugateGradientStabilizedOptions > | __options |
Definition at line 35 of file BiConjugateGradientStabilized.hpp.
| BiConjugateGradientStabilized::BiConjugateGradientStabilized | ( | const VectorType & | b, | |
| const MatrixType & | A, | |||
| const PreconditionerType & | P, | |||
| VectorType & | x | |||
| ) | [inline] |
Definition at line 47 of file BiConjugateGradientStabilized.hpp.
References __options, BiConjugateGradientStabilizedOptions::epsilon(), epsilon, ffout(), StaticBase< StreamCenter >::instance(), max_iter, BiConjugateGradientStabilizedOptions::maxiter(), and GetParameter< T >::value().
00051 { 00052 epsilon = __options.value().epsilon(); 00053 max_iter = __options.value().maxiter(); 00054 00055 ffout(2) << "- bi-conjugate gradient stabilized\n"; 00056 ffout(3) << " epsilon = " 00057 << std::setiosflags(std::ios_base::scientific) 00058 << epsilon 00059 << std::resetiosflags(std::ios_base::scientific) 00060 << '\n'; 00061 ffout(3) << " maximum number of iterations: " << max_iter << '\n'; 00062 00063 VectorType r_k_1(b.size()); 00064 A.timesX(x,r_k_1); 00065 r_k_1 = b - r_k_1; 00066 00067 real_t residu = Norm(r_k_1); 00068 00069 if(residu != 0) { 00070 real_t resid0 = residu; 00071 00072 VectorType rTilda_0 = r_k_1; 00073 VectorType p_k = r_k_1; 00074 00075 VectorType s_k(x.size()); 00076 00077 VectorType Ap_k(x.size()); 00078 VectorType As_k(x.size()); 00079 00080 VectorType r_k(x.size()); 00081 00082 ffout(2) << " initial residu: " 00083 << std::setiosflags(std::ios_base::scientific) 00084 << resid0 00085 << std::resetiosflags(std::ios_base::scientific) 00086 << '\n'; 00087 for (int i=1; i<= max_iter; ++i) { 00088 ffout(3) << " - iteration: " 00089 << std::setw(6) 00090 << i 00091 << "\tresidu: " 00092 << std::setiosflags(std::ios_base::scientific) 00093 << residu/resid0; 00094 ffout(4) << "\tabsolute: " << residu; 00095 ffout(3) << std::resetiosflags(std::ios_base::scientific) 00096 << '\n'; 00097 A.timesX(p_k,Ap_k); 00098 real_t alpha_k = (r_k_1*rTilda_0) / (Ap_k*rTilda_0); 00099 00100 s_k = r_k_1 - alpha_k*Ap_k; 00101 A.timesX(s_k,As_k); 00102 00103 real_t w_k = (As_k * s_k) / (As_k * As_k); 00104 00105 x += alpha_k * p_k + w_k * s_k; 00106 r_k = s_k - (w_k * As_k); 00107 00108 real_t beta_k = (r_k * rTilda_0)/(r_k_1 * rTilda_0) * (alpha_k / w_k); 00109 00110 p_k -= w_k*Ap_k; 00111 p_k *= beta_k; 00112 p_k += r_k; 00113 00114 if ((residu = (Norm(r_k)))/resid0<epsilon) { 00115 break; 00116 } 00117 00118 r_k_1 = r_k; 00119 } 00120 00121 if (residu/resid0 > epsilon) { 00122 ffout(2) << " bi_conjugate gradient stabilized: *NOT CONVERGED*\n"; 00123 ffout(2) << " - epsilon: " << epsilon << '\n'; 00124 ffout(2) << " - relative residu : " << residu/resid0 << '\n'; 00125 ffout(2) << " - absolute residu : " << residu << '\n'; 00126 } 00127 } 00128 00129 if (StreamCenter::instance().getDebugLevel()>3) { 00130 A.timesX(x,r_k_1); 00131 r_k_1 -= b; 00132 ffout(4) << "- final *real* residu : " << Norm(r_k_1) << '\n'; 00133 } 00134 ffout(2) << "- bi-conjugate gradient stabilized: done\n"; 00135 }
real_t BiConjugateGradientStabilized::epsilon [private] |
Definition at line 38 of file BiConjugateGradientStabilized.hpp.
Referenced by BiConjugateGradientStabilized().
int BiConjugateGradientStabilized::max_iter [private] |
Definition at line 39 of file BiConjugateGradientStabilized.hpp.
Referenced by BiConjugateGradientStabilized().
GetParameter<BiConjugateGradientStabilizedOptions> BiConjugateGradientStabilized::__options [private] |
Definition at line 41 of file BiConjugateGradientStabilized.hpp.
Referenced by BiConjugateGradientStabilized().
1.5.6