An Object Finite Element Library

CGS.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 Squared Method
                 CGS is inspired from the SIAM Templates book

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

#ifndef __CGS_H
#define __CGS_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 CGS(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();
   T_ rho_1=0, rho_2=0, alpha=0, beta=0;
   double res;
   Vect<T_> r(size), p(size), y(size), q(size), z(size), v(size), u(size), w(size);

   double nrm = b.getNorm2();
   A.Mult(x,r);
   Vect<T_> s = b - r;
   r = s;

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

   for (it=1; it<=max_it; it++) {
      rho_1 = Dot(s,r);
      if (rho_1 == T_(0)) {
         toler = r.getNorm2() / nrm;
         if (verbose)
            cout << "Convergence after 2 iterations." << endl;
         return 2;
      }
      if (it == 1) {
         u = r;
         p = u;
      } else {
         beta = rho_1 / rho_2;
         u = r + beta * q;
         p = u + beta * (q + beta * p);
      }
      y = p;
      P.Solve(y);
      A.Mult(y,v);
      alpha = rho_1 / Dot(s,v);
      w = u - alpha*v;
      w += u;
      P.Solve(w);
      Axpy(alpha,w,x);
      A.Mult(w,z);
      Axpy(-alpha,z,r);
      rho_2 = rho_1;
      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;
      }
   }
   toler = res;
   if (verbose)
      cout << "No Convergence after " << it << " iterations." << endl;
   return -max_it;
}

#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
int _CGS(Matrix<T_> *A, Prec<T_> *P, const Vect<T_> &b, Vect<T_> &x, int max_it,
         double &toler, int verbose)
{
   int it;
   size_t size=x.getSize();
   T_ rho_1=0, rho_2=0, alpha=0, beta=0;
   double res;
   Vect<T_> r(size), p(size), y(size), q(size), z(size), v(size), u(size), w(size);

   double nrm = b.getNorm2();
   A->Mult(x,r);
   Vect<T_> s = b - r;
   r = s;

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

   for (it=1; it<=max_it; it++) {
      rho_1 = Dot(s,r);
      if (rho_1 == T_(0)) {
         toler = r.getNorm2() / nrm;
         if (verbose)
            cout << "Convergence after 2 iterations." << endl;
         return 2;
      }
      if (it == 1) {
         u = r;
         p = u;
      } else {
         beta = rho_1 / rho_2;
         u = r + beta * q;
         p = u + beta * (q + beta * p);
      }
      y = p;
      P->Solve(y);
      A->Mult(y,v);
      alpha = rho_1 / Dot(s,v);
      w = u - alpha*v;
      w += u;
      P->Solve(w);
      Axpy(alpha,w,x);
      A->Mult(w,z);
      Axpy(-alpha,z,r);
      rho_2 = rho_1;
      res = r.getNorm2()/nrm;
      if (verbose > 1)
         cout << "Iteration: " << setw(4) << it << "  ... Residual: " << res << endl;
      if (verbose > 2)
         cout << x;
      if (res < toler) {
         toler = res;
         if (verbose)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
   }
   toler = res;
   if (verbose)
      cout << "No Convergence after " << it << " iterations." << endl;
   return -max_it;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani