#include <BiConjugateGradient.hpp>
Public Member Functions | |
| template<typename VectorType, typename MatrixType, typename PreconditionerType> | |
| BiConjugateGradient (const VectorType &b, const MatrixType &A, PreconditionerType &P, VectorType &u) | |
Private Attributes | |
| real_t | epsilon |
| int | max_iter |
| GetParameter < BiConjugateGradientOptions > | __options |
#0 using the preconditionner P.Definition at line 35 of file BiConjugateGradient.hpp.
| BiConjugateGradient::BiConjugateGradient | ( | const VectorType & | b, | |
| const MatrixType & | A, | |||
| PreconditionerType & | P, | |||
| VectorType & | u | |||
| ) | [inline] |
Definition at line 48 of file BiConjugateGradient.hpp.
References __options, BiConjugateGradientOptions::epsilon(), epsilon, ffout(), StaticBase< StreamCenter >::instance(), max_iter, BiConjugateGradientOptions::maxiter(), and GetParameter< T >::value().
00052 { 00053 epsilon = __options.value().epsilon(); 00054 max_iter = __options.value().maxiter(); 00055 00056 ffout(2) << "- bi-conjugate gradient\n"; 00057 ffout(3) << " epsilon = " 00058 << std::setiosflags(std::ios_base::scientific) 00059 << epsilon 00060 << std::resetiosflags(std::ios_base::scientific) 00061 << '\n'; 00062 ffout(3) << " maximum number of iterations: " << max_iter << '\n'; 00063 00064 real_t alpha, beta; 00065 real_t rho_1 = 0; 00066 real_t rho_2 = 1; 00067 00068 VectorType z(b.size()); 00069 VectorType zTilda(b.size()); 00070 VectorType p(b.size()); 00071 VectorType pTilda(b.size()); 00072 VectorType q(b.size()); 00073 VectorType qTilda(b.size()); 00074 z = 0; 00075 p = 0; 00076 00077 VectorType r = b; 00078 A.timesX(u,q); 00079 r -= q; 00080 VectorType rTilda = r; 00081 00082 real_t residu = Norm(r); 00083 if (residu==0) 00084 residu = 1.; 00085 real_t resid0 = residu; 00086 00087 ffout(2) << " initial residu: " 00088 << std::setiosflags(std::ios_base::scientific) 00089 << resid0 00090 << std::resetiosflags(std::ios_base::scientific) 00091 << '\n'; 00092 for (size_t i = 1; i <= (size_t)max_iter; i++) { 00093 ffout(3) << " - iteration: " 00094 << std::setw(6) 00095 << i 00096 << "\tresidu: " 00097 << std::setiosflags(std::ios_base::scientific) 00098 << residu/resid0; 00099 ffout(4) << "\tabsolute: " << residu; 00100 ffout(3) << std::resetiosflags(std::ios_base::scientific) 00101 << '\n'; 00102 // z = P^{-1}(r) 00103 P.computes(r,z); 00104 00105 // zTilda = P^{-T}(rTilda) Here we only use symetric preconditionners. 00106 P.computesTransposed(rTilda,zTilda); 00107 00108 rho_1 = z * rTilda; 00109 if (rho_1 == 0) { 00110 break; 00111 } 00112 00113 if (i == 1) { 00114 p = z; 00115 pTilda = zTilda; 00116 } else { 00117 beta = rho_1/rho_2; 00118 // p = z + beta * p; 00119 p *= beta; 00120 p += z; 00121 00122 // pTilda = zTilda + beta * pTilda; 00123 pTilda *= beta; 00124 pTilda += zTilda; 00125 } 00126 00127 A.transposedTimesX(pTilda,qTilda); 00128 A.timesX(p,q); 00129 alpha = rho_1/(pTilda*q); 00130 u += alpha * p; 00131 r -= alpha * q; 00132 rTilda -= alpha * qTilda; 00133 00134 rho_2 = rho_1; 00135 if ((residu = (Norm(r)))/resid0<epsilon) { 00136 break; 00137 } 00138 } 00139 00140 if (residu/resid0 > epsilon) { 00141 ffout(2) << " bi_conjugate gradient: *NOT CONVERGED*\n"; 00142 ffout(2) << " - epsilon: " << epsilon << '\n'; 00143 ffout(2) << " - relative residu : " << residu/resid0 << '\n'; 00144 ffout(2) << " - absolute residu : " << residu << '\n'; 00145 } 00146 00147 if (StreamCenter::instance().getDebugLevel()>3) { 00148 A.timesX(u,r); 00149 r -= b; 00150 ffout(4) << "- final *real* residu : " << Norm(r) << '\n'; 00151 } 00152 ffout(2) << "- bi-conjugate gradient: finished\n"; 00153 }
real_t BiConjugateGradient::epsilon [private] |
int BiConjugateGradient::max_iter [private] |
1.5.6