ConjugateGradient Class Reference
#include <ConjugateGradient.hpp>
Collaboration diagram for ConjugateGradient:
Public Member Functions | |
| template<typename VectorType, typename MatrixType, typename PreconditionerType> | |
| ConjugateGradient (const VectorType &f, const MatrixType &A, const PreconditionerType &C, VectorType &x, const int maxiter=-1, const real_t epsilon=-1.) | |
Private Attributes | |
| real_t | __epsilon |
| int | __maxiter |
| GetParameter < ConjugateGradientOptions > | __options |
Detailed Description
Conjugate Gradient.
- Resolve the linear system
#0 using the preconditionner P.
Definition at line 36 of file ConjugateGradient.hpp.
Constructor & Destructor Documentation
template<typename VectorType, typename MatrixType, typename PreconditionerType>
| ConjugateGradient::ConjugateGradient | ( | const VectorType & | f, | |
| const MatrixType & | A, | |||
| const PreconditionerType & | C, | |||
| VectorType & | x, | |||
| const int | maxiter = -1, |
|||
| const real_t | epsilon = -1. | |||
| ) | [inline] |
Definition at line 48 of file ConjugateGradient.hpp.
References __epsilon, __maxiter, __options, ConjugateGradientOptions::epsilon(), ffout(), StaticBase< StreamCenter >::instance(), ConjugateGradientOptions::maxiter(), and GetParameter< T >::value().
00054 { 00055 if (epsilon == -1) { 00056 __epsilon = __options.value().epsilon(); 00057 } else { 00058 __epsilon = epsilon; 00059 } 00060 00061 if (maxiter == -1) { 00062 __maxiter = __options.value().maxiter(); 00063 } else { 00064 __maxiter = maxiter; 00065 } 00066 00067 ffout(2) << "- conjugate gradient\n"; 00068 ffout(3) << " epsilon = " 00069 << std::setiosflags(std::ios_base::scientific) 00070 << __epsilon 00071 << std::resetiosflags(std::ios_base::scientific) 00072 << '\n'; 00073 ffout(3) << " maximum number of iterations: " << __maxiter << '\n'; 00074 00075 VectorType h(x.size()); 00076 VectorType b(f); 00077 00078 VectorType g(b.size()); 00079 VectorType cg = b; 00080 00081 real_t gcg=0; 00082 real_t gcg0=1; 00083 00084 real_t relativeEpsilon = __epsilon; 00085 00086 for (int i=1; i<=__maxiter; ++i) { 00087 if (i==1) { 00088 A.timesX(x,h); 00089 00090 cg -= h; 00091 00092 C.computes(cg, g); 00093 00094 gcg = g * cg; 00095 00096 h=g; 00097 } 00098 00099 A.timesX(h,b); 00100 00101 real_t hAh = h*b; 00102 00103 if (hAh==0) { 00104 hAh=1.; 00105 } 00106 real_t ro = gcg/hAh; 00107 cg -= ro*b; 00108 00109 VectorType b2 = b; 00110 C.computes(b2,b); 00111 00112 x+=ro*h; 00113 g-=ro*b; 00114 00115 real_t gamma=gcg; 00116 gcg = g * cg; 00117 00118 if ((i == 1)&&(gcg != 0)) { 00119 relativeEpsilon = __epsilon*gcg; 00120 gcg0=gcg; 00121 ffout(2) << " initial residu: " 00122 << std::setiosflags(std::ios_base::scientific) 00123 << gcg 00124 << std::resetiosflags(std::ios_base::scientific) 00125 << '\n'; 00126 } 00127 ffout(3) << " - iteration " 00128 << std::setw(6) 00129 << i 00130 << "\tresidu: " 00131 << std::setiosflags(std::ios_base::scientific) 00132 << gcg/gcg0; 00133 ffout(4) << "\tabsolute: " << gcg; 00134 ffout(3) << std::resetiosflags(std::ios_base::scientific) 00135 << '\n'; 00136 00137 if (gcg<relativeEpsilon) { 00138 break; 00139 } 00140 00141 gamma=gcg/gamma; 00142 00143 h *= gamma; 00144 h += g; 00145 } 00146 00147 if (gcg > relativeEpsilon) { 00148 ffout(2) << " conjugate gradient: *NOT CONVERGED*\n"; 00149 ffout(2) << " - epsilon: " << __epsilon << '\n'; 00150 ffout(2) << " - relative residu : " << gcg/gcg0 << '\n'; 00151 ffout(2) << " - absolute residu : " << gcg << '\n'; 00152 } 00153 00154 if (StreamCenter::instance().getDebugLevel()>3) { 00155 A.timesX(x,h); 00156 h -= b; 00157 ffout(4) << "- final *real* residu : " << h*h << '\n'; 00158 } 00159 ffout(2) << "- conjugate gradient: done\n"; 00160 }
Here is the call graph for this function:
Member Data Documentation
real_t ConjugateGradient::__epsilon [private] |
int ConjugateGradient::__maxiter [private] |
The documentation for this class was generated from the following file:
