BiConjugateGradient.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: BiConjugateGradient.hpp,v 1.5 2008/02/24 20:45:17 delpinux Exp $

#ifndef BI_CONJUGATE_GRADIENT_HPP
#define BI_CONJUGATE_GRADIENT_HPP

#include <BiConjugateGradientOptions.hpp>
#include <GetParameter.hpp>

class BiConjugateGradient
{
 private:
  real_t epsilon;
  int    max_iter;

  GetParameter<BiConjugateGradientOptions> __options;

 public:

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

    ffout(2) << "- bi-conjugate gradient\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';

    real_t alpha, beta;
    real_t rho_1 = 0;
    real_t rho_2 = 1;

    VectorType z(b.size());
    VectorType zTilda(b.size());
    VectorType p(b.size());
    VectorType pTilda(b.size());
    VectorType q(b.size());
    VectorType qTilda(b.size());
    z = 0;
    p = 0;

    VectorType r = b;
    A.timesX(u,q);
    r -= q;
    VectorType rTilda = r;

    real_t residu = Norm(r);
    if (residu==0)
      residu = 1.;
    real_t resid0 = residu;

    ffout(2) << "   initial residu: "
             << std::setiosflags(std::ios_base::scientific)
             << resid0
             << std::resetiosflags(std::ios_base::scientific)
             << '\n';
    for (size_t i = 1; i <= (size_t)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';
      // z = P^{-1}(r)
      P.computes(r,z);

      // zTilda = P^{-T}(rTilda) Here we only use symetric preconditionners.
      P.computesTransposed(rTilda,zTilda);

      rho_1 = z * rTilda;
      if (rho_1 == 0) { 
        break;
      }

      if (i == 1) {
        p = z;
        pTilda = zTilda;
      } else {
        beta = rho_1/rho_2;
        // p = z + beta * p;
        p *= beta;
        p += z;

        // pTilda = zTilda + beta * pTilda;
        pTilda *= beta;
        pTilda += zTilda;
      }

      A.transposedTimesX(pTilda,qTilda);
      A.timesX(p,q);
      alpha = rho_1/(pTilda*q);
      u += alpha * p;
      r -= alpha * q;
      rTilda -= alpha * qTilda;

      rho_2 = rho_1;
      if ((residu = (Norm(r)))/resid0<epsilon) {
        break;
      }
    }

    if (residu/resid0 > epsilon) {
      ffout(2) << "  bi_conjugate gradient: *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(u,r);
      r -= b;
      ffout(4) << "- final *real* residu :   " << Norm(r) << '\n';
    }
    ffout(2) << "- bi-conjugate gradient: finished\n";
  }

};
  
#endif // BI_CONJUGATE_GRADIENT_HPP

---