BiConjugateGradient Class Reference

#include <BiConjugateGradient.hpp>

Collaboration diagram for BiConjugateGradient:

Collaboration graph
[legend]

List of all members.

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


Detailed Description

Bi Conjugate Gradient.

Resolve the linear system $ p\in T $#0 using the preconditionner P.
Author:
Stéphane Del Pino.

Definition at line 35 of file BiConjugateGradient.hpp.


Constructor & Destructor Documentation

template<typename VectorType, typename MatrixType, typename PreconditionerType>
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   }

Here is the call graph for this function:


Member Data Documentation

real_t BiConjugateGradient::epsilon [private]

Definition at line 38 of file BiConjugateGradient.hpp.

Referenced by BiConjugateGradient().

Definition at line 39 of file BiConjugateGradient.hpp.

Referenced by BiConjugateGradient().

Definition at line 41 of file BiConjugateGradient.hpp.

Referenced by BiConjugateGradient().


The documentation for this class was generated from the following file:

Generated on Wed Nov 19 00:04:06 2008 for FreeFEM3D (aka ff3d) by  doxygen 1.5.6