BiConjugateGradientStabilized Class Reference

#include <BiConjugateGradientStabilized.hpp>

Collaboration diagram for BiConjugateGradientStabilized:

Collaboration graph
[legend]

List of all members.

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


Detailed Description

Bi-Conjugate Gradient Stabilized.

Resolve the linear system #0 using the preconditionner P.
Author:
Stéphane Del Pino.

Definition at line 35 of file BiConjugateGradientStabilized.hpp.


Constructor & Destructor Documentation

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

Here is the call graph for this function:


Member Data Documentation

Definition at line 38 of file BiConjugateGradientStabilized.hpp.

Referenced by BiConjugateGradientStabilized().

Definition at line 39 of file BiConjugateGradientStabilized.hpp.

Referenced by BiConjugateGradientStabilized().

Definition at line 41 of file BiConjugateGradientStabilized.hpp.

Referenced by BiConjugateGradientStabilized().


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

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