An Object Finite Element Library

CG.h

Go to the documentation of this file.

/*==============================================================================

                                    O  F  E  L  I

                           Object  Finite  Element  Library

  ==============================================================================

   Copyright (C) 1998 - 2008 Rachid Touzani

   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; Version 2 of the License.

   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

  ==============================================================================

                 Template function for Conjugate Gradient Method
                   CG is inspired from the SIAM Templates book

  ==============================================================================*/

#ifndef __CG_H
#define __CG_H

#include <iostream>
using std::ostream;
using std::endl;

#include <iomanip>
using std::setw;
using std::ios;
using std::setprecision;


#include "OFELI_Config.h"
#include "Vect.h"
#include "Precond.h"
#include "util.h"
#include "output.h"

namespace OFELI {

template<class T_, class M_, class P_ >
int CG(const M_ &A, const P_ &P, const Vect<T_> &b, Vect<T_> &x, int max_it,
       double &toler, int verbose)
{
   int it;
   size_t size = x.getSize();
   double res, nrm=b.getNorm2();
   T_ rho=0, rho_1=T_(1), beta=0;

   if (nrm == 0)
      nrm = 1;
   Vect<T_> r(size), p(size), q(size), z(size);
   A.Mult(x,r);
   r = b - r;
   if ((res = r.getNorm2()/nrm) <= toler) {
      if (verbose > 1)
         cout << "Convergence after 0 iterations." << endl;
      return 0;
   }

   for (it=1; it<=max_it; it++) {
       P.Solve(r,z);
       rho = Dot(r,z);
       if (it==1)
          p = z;
       else {
          beta = rho/rho_1;
          p = z + beta * p;
       }
       A.Mult(p,q);
       T_ alpha = rho/Dot(p,q);
       x += alpha * p;
       r -= alpha * q;
       res = r.getNorm2()/nrm;

       if (verbose > 1) {
          cout << "Iteration: " << setw(4) << it;
          cout << "  ... Residual: " << res << endl;
       }
       if (verbose > 2)
          cout << x;
       if (res <= toler) {
          toler = res;
          if (verbose)
             cout << "Convergence after " << it << " iterations." << endl;
          return it;
       }
       rho_1 = rho;
   }
   if (verbose)
      cout << "No Convergence after " << it << " iterations." << endl;
   return -it;
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
int _CG(Matrix<T_> *A, Prec<T_> *P, const Vect<T_> &b,
        Vect<T_> &x, int max_it, double toler, int verbose)
{
   int it=0;
   size_t size = x.getSize();
   double nrm=b.getNorm2(), rho=0, res=0, beta=0, rho_1=0;

   if (nrm == 0)
      nrm = 1;
   Vect<T_> r(size), z(size), p(size), q(size);
   A->Mult(x,r);
   r = b - r;
   if ((res = r.getNorm2()/nrm) <= toler) {
      if (verbose > 1)
         cout << "Convergence after " << it << " iterations." << endl;
      return it;
   }

   for (it=1; it<=max_it; it++) {
      P->Solve(r,z);
      rho = Dot(r,z);
      if (it==1)
         p = z;
      else {
         beta = rho/rho_1;
         p = z + beta * p;
      }
      A->Mult(p,q);
      double alpha = rho/Dot(p,q);
      x += alpha * p;
      r -= alpha * q;
      res = r.getNorm2()/nrm;

      if (verbose > 1) {
         cout << "Iteration : " << setw(4) << it;
         cout << "  ... Residual : " << res << endl;
      }
      if (verbose > 2)
         cout << x;
      if (res <= toler) {
         if (verbose)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
      rho_1 = rho;
   }
   if (verbose)
      cout << "No Convergence after " << it << " iterations." << endl;
   return -it;
}


template<class T_, class M_, class P_>
int CG(const P_ &P, const Vect<T_> &b, Vect<T_> &x, int max_it, double toler,
       int verbose)
{
   int it=0;
   size_t size = x.Dim();
   double nrm=b.getNorm2(), rho=0, res=0, beta=0, rho_1=0;
   if (nrm == 0)
      nrm = 1;
   Vect<T_> r(size), z(size), p(size), q(size);
   r = b;
   r.MultAdd(Mxv(x),-1);

   if ((res = r.getNorm2() / nrm) <= toler) {
      if (verbose)
         cout << "Convergence after " << it << " iterations." << endl;
      return it;
   }

   for (it=1; it<=max_it; it++) {
      P.Solve(r,z);
      rho = Dot(r,z);
      if (it==1)
         p = z;
      else {
         beta = rho/rho_1;
         p *= beta;
         p.MultAdd(z);
      }

      q = Mxv(p);
      double alpha = rho/Dot(p,q);
      x.MultAdd(p, alpha);
      r.MultAdd(q,-alpha);
      res = r.Norm2()/nrm;

      if (verbose > 1)
         cout << "\nIteration : " << it << "  ..." << endl;
      if (verbose > 2)
         cout << x;
      if (res <= toler) {
         if (verbose)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
      rho_1 = rho;
   }
   if (verbose)
      cout << "No Convergence after " << it << " iterations." << endl;
   return -it;
}


template<class T_>
int _CG(const Prec<T_> *P, const Vect<T_> &b, Vect<T_> &x, int max_it,
        double toler, int verbose)
{
   int it=0;
   size_t size = x.Dim();
   double nrm=b.Norm2(), rho=0, res=0, beta=0, rho_1=0;
   if (nrm == 0)
      nrm = 1;
   Vect<T_> r(size), z(size), p(size), q(size);
   r = b;
   r.MultAdd(Mxv(x),-1);

   if ((res = r.Norm2() / nrm) <= toler) {
      if (verbose)
         cout << "Convergence after " << it << " iterations." << endl;
      return it;
   }

   for (it=1; it<=max_it; it++) {
      P->Solve(r,z);
      rho = Dot(r,z);
      if (it==1)
         p = z;
      else {
         beta = rho/rho_1;
         p *= beta;
         p.MultAdd(z);
      }

      q = Mxv(p);
      double alpha = rho/Dot(p,q);
      x.MultAdd(p, alpha);
      r.MultAdd(q,-alpha);
      res = r.Norm2()/nrm;

      if (verbose > 1)
         cout << "\nIteration : " << it << "  ..." << endl;
      if (verbose > 2)
         cout << x;
      if (res <= toler) {
         if (verbose)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
      rho_1 = rho;
   }
   if (verbose)
      cout << "No Convergence after " << it << " iterations." << endl;
   return -it;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani