GMRES.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: GMRES.hpp,v 1.4 2008/03/01 18:05:57 delpinux Exp $


#ifndef GMRES_HPP
#define GMRES_HPP

#include <GMRESOptions.hpp>
#include <GetParameter.hpp>


class GMRES
{
private:
  GetParameter<GMRESOptions>
  __options;                    
  const real_t __epsilon;       
  const int __maxiter;          
  const size_t __m;             
  void __givens(const real_t& beta,
                Vector<Vector<real_t> >& H,
                Vector<real_t>& y)
  {
    ffout(3) << std::setiosflags(std::ios_base::scientific);
    
    y = 0;
    Vector<real_t> b(H.size());
    b = 0;
    b[0]=beta;

    for (size_t i=0; i<H.size()-1; ++i) {
      const real_t h0 = H[i][i];
      const real_t h1 = H[i+1][i];
      const real_t r = std::sqrt(h0*h0 + h1*h1);
      const real_t c = h0/r;
      const real_t s = h1/r;

      for (size_t j=i; j<H.size()-1; ++j) {
        real_t& h_i_j   = H[i][j];
        real_t& h_ip1_j = H[i+1][j];
        
        const real_t newH_i_j   = c*h_i_j   + s*h_ip1_j;
        const real_t newH_ip1_j = c*h_ip1_j - s*h_i_j;

        h_i_j   = newH_i_j;
        h_ip1_j = newH_ip1_j;
      }

      const real_t newB_i   = c*b[i]   + s*b[i+1];
      const real_t newB_ip1 = c*b[i+1] - s*b[i];

      b[i]   = newB_i;
      b[i+1] = newB_ip1;
    }

    for (int i = H.size()-2; i>=0; --i) {
      real_t sum = b[i];
      for (int j=H.size()-2; j>i; --j) {
        sum -= H[i][j]*y[j];
      }
      y[i] = sum / H[i][i];
    }

    ffout(3) << std::resetiosflags(std::ios_base::scientific);
  }

public:
  template <typename VectorType,
            typename MatrixType,
            typename PreconditionerType>
  GMRES(const VectorType& b,
        const MatrixType& A, 
        const PreconditionerType& P,
        VectorType& x)
    : __epsilon(__options.value().epsilon()),
      __maxiter(__options.value().maxiter()),
      __m(__options.value().basis())
  {
    ffout(2) << "- general minimum residual method (m)\n";
    ffout(3) << "  epsilon = "
             << std::setiosflags(std::ios_base::scientific)
             << __epsilon
             << std::resetiosflags(std::ios_base::scientific)
             << '\n';
    ffout(3) << "  number of Krylov basis vectors: " << __m << '\n';
    ffout(3) << "  maximum number of iterations: " << __maxiter << '\n';


    VectorType w(b.size());
    VectorType z(b.size());

    Vector<VectorType> V(__m+1);
    for (size_t i=0; i<V.size(); ++i) {
      V[i].resize(b.size());
    }

    // create the Hessenberg matrix H
    Vector<Vector<real_t> > H;
    H.resize(__m+1);
    for (size_t k=0; k< H.size(); ++k) {
      H[k].resize(__m);
      H[k]=0;
    }

    real_t normb = 0;
    {
      VectorType P_1b(b.size());
      P.computes(b,P_1b);
      normb = Norm(P_1b);
    }

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

    VectorType r(b.size());

    P.computes(Pr,r);

    real_t normPr = Norm(Pr);

    ffout(2) << "  initial residu: "
             << std::setiosflags(std::ios_base::scientific)
             << normPr
             << std::resetiosflags(std::ios_base::scientific)
             << '\n';

    bool converged = (normPr<=__epsilon*normb);
    if (not converged) {
      for (int iteration=1; iteration<=__maxiter; ++iteration) {
        real_t beta = Norm(r);
        V[0] = (1./beta) * r;

        for (size_t j=0; j<__m; ++j) {
          A.timesX(V[j],z);
          P.computes(z, w);

          for (size_t i=0; i<=j; ++i) {
            H[i][j] = w * V[i];
            w -= H[i][j] * V[i];
          }
          H[j+1][j] = Norm(w);
          if (H[j+1][j] != 0)
            V[j+1] = (1./H[j+1][j]) * w;
          else {
            throw ErrorHandler(__FILE__,__LINE__,
                               "singular matrix!",
                               ErrorHandler::normal);
          }
        }

        Vector<real_t> y(__m);
        this->__givens(beta, H, y);

        for (size_t j=0; j<__m; ++j) {
          x += y[j]*V[j];
        }

        A.timesX(x,Ax);
        Pr = b;
        Pr -= Ax;

        normPr = Norm(Pr);

        ffout(3) << "  - iteration "
                 << std::setw(6)
                 << iteration
                 << "\tresidu: "
                 << std::setiosflags(std::ios_base::scientific)
                 << normPr/normb;
        ffout(4) << "\tabsolute: " << normPr;
        ffout(3) << std::resetiosflags(std::ios_base::scientific)
                 << '\n';

        if (normPr < normb * __epsilon) {
          converged = true;
          break;
        }

        P.computes(Pr,r);
      }
    }

    if (not converged) {
      ffout(2) << " general minimum residual : *NOT CONVERGED*\n";
      ffout(2) << "  - m:                " << __m << '\n';
      ffout(2) << "  - epsilon:          " << __epsilon << '\n';
      ffout(2) << "  - relative residu : " << normPr/normb << '\n';
      ffout(2) << "  - absolute residu : " << normPr << '\n';
    }

    if (StreamCenter::instance().getDebugLevel()>3) {
      A.timesX(x,Ax);
      Ax -= b;
      ffout(4) << "- final *real* residu :   " << Ax*Ax << '\n';
    }

    ffout(2) << "- general minimum residual method (m): done\n";
  }
};

#endif // GMRES_HPP

---