An Object Finite Element Library

BiCGStab.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 BiConjugate Gradient Stabilized Method
               BiCGStab is inspired from the SIAM Templates book.

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

#ifndef __BICG_STAB_H
#define __BICG_STAB_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 BiCGStab(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, omega=0;
   Vect<T_> p(size), pp(size), s(size), ss(size), t(size), v(size), r(size), rr(size);
   double res, nrm = b.getNorm2();
   A.Mult(x,r);
   r = b - r;
   rr = r;

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

   for (it=1; it<=max_it; it++) {
      rho_1 = Dot(rr,r);
      if (rho_1 == T_(0)) {
        toler = r.getNorm2() / nrm;
        return 2;
      }
      if (it == 1)
         p = r;
      else {
         beta = (rho_1/rho_2) * (alpha/omega);
         p = r + beta*(p - omega*v);
      }
      P.Solve(p,pp);
      A.Mult(pp,v);
      alpha = rho_1/Dot(rr,v);
      s = r - alpha * v;
      if ((res = s.getNorm2()/nrm) < toler) {
         x += alpha*pp;
         if (verbose > 1) {
            cout << "Iteration : " << setw(4) << it;
            cout << "  ... Residual : " << res << endl;
         }
         toler = res;
         return it;
      }
      P.Solve(s,ss);
      A.Mult(ss,t);
      omega = Dot(t,s)/Dot(t,t);
      x += alpha*pp + omega*ss;
      r = s - omega*t;

      rho_2 = rho_1;
      if (verbose > 1) {
         cout << "Iteration : " << setw(4) << it;
         cout << "  ... Residual : " << res << endl;
      }
      if (verbose > 2)
         cout << x;
      if ((res = r.getNorm2() / nrm) < toler) {
         toler = res;
         if (verbose)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
      if (omega == T_(0)) {
         toler = r.getNorm2() / nrm;
         return 3;
      }
   }
   if (verbose)
      cout << "Convergence after " << it << " iterations." << endl;
   toler = res;
   return -it;
}


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
int _BiCGStab(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, omega=0;
   Vect<T_> p(size), pp(size), s(size), ss(size), t(size), v(size), r(size), rr(size);
   double res, nrm = b.getNorm2();
   A->Mult(x,r);
   r = b - r;
   rr = r;

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

   for (it=1; it<=max_it; it++) {
      rho_1 = Dot(rr,r);
      if (rho_1 == T_(0)) {
        toler = r.getNorm2() / nrm;
        return 2;
      }
      if (it == 1)
         p = r;
      else {
         beta = (rho_1/rho_2) * (alpha/omega);
         p = r + beta*(p - omega*v);
      }
      P->Solve(p,pp);
      A->Mult(pp,v);
      alpha = rho_1/Dot(rr,v);
      s = r - alpha * v;
      if ((res = s.getNorm2()/nrm) < toler) {
         x += alpha*pp;
         if (verbose > 1) {
            cout << "Iteration : " << setw(4) << it;
            cout << "  ... Residual : " << res << endl;
         }
         toler = res;
         return it;
      }
      P->Solve(s,ss);
      A->Mult(ss,t);
      omega = Dot(t,s)/Dot(t,t);
      x += alpha*pp + omega*ss;
      r = s - omega*t;

      rho_2 = rho_1;
      if (verbose > 1) {
         cout << "Iteration : " << setw(4) << it;
         cout << "  ... Residual : " << res << endl;
      }
      if (verbose > 2)
         cout << x;
      if ((res = r.getNorm2() / nrm) < toler) {
         toler = res;
         if (verbose)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
      if (omega == T_(0)) {
         toler = r.getNorm2() / nrm;
         return 3;
      }
   }
   if (verbose)
      cout << "Convergence after " << it << " iterations." << endl;
   toler = res;
   return -it;
}
#endif /* DOXYGEN_SHOULD_SKIP_THIS */

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani