BiConjugateGradientStabilized.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: BiConjugateGradientStabilized.hpp,v 1.4 2007/08/29 00:31:04 delpinux Exp $

#ifndef BI_CONJUGATE_GRADIENT_STABILIZED_HPP
#define BI_CONJUGATE_GRADIENT_STABILIZED_HPP

#include <BiConjugateGradientStabilizedOptions.hpp>
#include <GetParameter.hpp>

class BiConjugateGradientStabilized
{
 private:
  real_t epsilon;
  int    max_iter;

  GetParameter<BiConjugateGradientStabilizedOptions> __options;

 public:
  template <typename VectorType,
            typename MatrixType,
            typename PreconditionerType>
  BiConjugateGradientStabilized(const VectorType& b,
                                const MatrixType& A, 
                                const PreconditionerType& P,
                                VectorType& x)
  {
    epsilon = __options.value().epsilon();
    max_iter = __options.value().maxiter();

    ffout(2) << "- bi-conjugate gradient stabilized\n";
    ffout(3) << "  epsilon = "
             << std::setiosflags(std::ios_base::scientific)
             << epsilon
             << std::resetiosflags(std::ios_base::scientific)
             << '\n';
    ffout(3) << "  maximum number of iterations: " << max_iter << '\n';

    VectorType r_k_1(b.size());
    A.timesX(x,r_k_1);
    r_k_1 = b - r_k_1;

    real_t residu = Norm(r_k_1);

    if(residu != 0) {
      real_t resid0 = residu;

      VectorType rTilda_0 = r_k_1;
      VectorType p_k = r_k_1;

      VectorType s_k(x.size());

      VectorType Ap_k(x.size());
      VectorType As_k(x.size());

      VectorType r_k(x.size());

      ffout(2) << "   initial residu: "
               << std::setiosflags(std::ios_base::scientific)
               << resid0
               << std::resetiosflags(std::ios_base::scientific)
               << '\n';
      for (int i=1; i<= max_iter; ++i) {
        ffout(3) << "  - iteration: "
                 << std::setw(6)
                 << i
                 << "\tresidu: "
                 << std::setiosflags(std::ios_base::scientific)
                 << residu/resid0;
        ffout(4) << "\tabsolute: " << residu;
        ffout(3) << std::resetiosflags(std::ios_base::scientific)
                 << '\n';
        A.timesX(p_k,Ap_k);
        real_t alpha_k = (r_k_1*rTilda_0) / (Ap_k*rTilda_0);

        s_k = r_k_1 - alpha_k*Ap_k;
        A.timesX(s_k,As_k);

        real_t w_k = (As_k * s_k) / (As_k * As_k);

        x += alpha_k * p_k + w_k * s_k;
        r_k = s_k - (w_k * As_k);

        real_t beta_k = (r_k * rTilda_0)/(r_k_1 * rTilda_0) * (alpha_k / w_k);

        p_k -= w_k*Ap_k;
        p_k *= beta_k;
        p_k += r_k;

        if ((residu = (Norm(r_k)))/resid0<epsilon) {
          break;
        }

        r_k_1 = r_k;
      }

      if (residu/resid0 > epsilon) {
        ffout(2) << "  bi_conjugate gradient stabilized: *NOT CONVERGED*\n";
        ffout(2) << "  - epsilon:          " << epsilon << '\n';
        ffout(2) << "  - relative residu : " << residu/resid0 << '\n';
        ffout(2) << "  - absolute residu : " << residu << '\n';
      }
    }

    if (StreamCenter::instance().getDebugLevel()>3) {
      A.timesX(x,r_k_1);
      r_k_1 -= b;
      ffout(4) << "- final *real* residu :   " << Norm(r_k_1) << '\n';
    }
    ffout(2) << "- bi-conjugate gradient stabilized: done\n";
  }
};
  
#endif // BI_CONJUGATE_GRADIENT_STABILIZED_HPP

---