ConjugateGradient.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: ConjugateGradient.hpp,v 1.5 2007/11/02 20:37:27 delpinux Exp $


#ifndef CONJUGATE_GRADIENT_HPP
#define CONJUGATE_GRADIENT_HPP

#include <ConjugateGradientOptions.hpp>
#include <GetParameter.hpp>

class ConjugateGradient
{
private:
  real_t __epsilon;
  int    __maxiter;

  GetParameter<ConjugateGradientOptions> __options;

public:
  template <typename VectorType,
            typename MatrixType,
            typename PreconditionerType>
  ConjugateGradient(const VectorType& f,
                    const MatrixType& A, 
                    const PreconditionerType& C,
                    VectorType& x,
                    const int maxiter = -1,
                    const real_t epsilon = -1.)
  {
    if (epsilon == -1) {
      __epsilon  = __options.value().epsilon();
    } else {
      __epsilon = epsilon;
    }

    if (maxiter == -1) {
      __maxiter = __options.value().maxiter();
    } else {
      __maxiter = maxiter;
    }

    ffout(2) << "- 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: " << __maxiter << '\n';

    VectorType h(x.size());
    VectorType b(f);

    VectorType g(b.size());
    VectorType cg = b;

    real_t gcg=0;
    real_t gcg0=1;

    real_t relativeEpsilon = __epsilon;

    for (int i=1; i<=__maxiter; ++i) {
      if (i==1) {
        A.timesX(x,h);

        cg -= h;

        C.computes(cg, g);

        gcg = g * cg;

        h=g;
      }

      A.timesX(h,b);

      real_t hAh = h*b;

      if (hAh==0) {
        hAh=1.;
      }
      real_t ro = gcg/hAh;
      cg -= ro*b;

      VectorType b2 = b;
      C.computes(b2,b);
      
      x+=ro*h;
      g-=ro*b;

      real_t gamma=gcg;
      gcg = g * cg;

      if ((i == 1)&&(gcg != 0)) {
        relativeEpsilon = __epsilon*gcg;
        gcg0=gcg;
        ffout(2) << "  initial residu: "
                 << std::setiosflags(std::ios_base::scientific)
                 << gcg
                 << std::resetiosflags(std::ios_base::scientific)
                 << '\n';
      }
      ffout(3) << "  - iteration "
               << std::setw(6)
               << i
               << "\tresidu: "
               << std::setiosflags(std::ios_base::scientific)
               << gcg/gcg0;
      ffout(4) << "\tabsolute: " << gcg;
      ffout(3) << std::resetiosflags(std::ios_base::scientific)
               << '\n';

      if (gcg<relativeEpsilon) {
        break;
      }

      gamma=gcg/gamma;

      h *= gamma;
      h += g;
    }

    if (gcg > relativeEpsilon) {
      ffout(2) << "  conjugate gradient: *NOT CONVERGED*\n";
      ffout(2) << "  - epsilon:          " << __epsilon << '\n';
      ffout(2) << "  - relative residu : " << gcg/gcg0 << '\n';
      ffout(2) << "  - absolute residu : " << gcg << '\n';
    }

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

---