An Object Finite Element Library

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

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

#ifndef __BICG_H
#define __BICG_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 BiCG(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;
   T_ alpha=0, beta=0, rho_1=0, rho_2=0;
   Vect<T_> r(size), z(size), zz(size), p(size), pp(size), q(size), qq(size);

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

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

   for (it=1; it<=max_it; it++) {
      P.Solve(r,z);
      P.TransSolve(rr,zz);
      rho_1 = Dot(z,rr);
      if (rho_1 == T_(0)) {
         toler = r.getNorm2() / nrm;
         if (verbose > 1)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
      if (it == 1) {
         p = z;
         pp = zz;
      } else {
         beta = rho_1 / rho_2;
         p = z + beta*p;
         pp = zz + beta*pp;
      }
      A.Mult(p,q);
      A.TMult(pp,qq);
      alpha = rho_1/Dot(pp,q);
      Axpy( alpha,p,x);
      Axpy(-alpha,q,r);
      Axpy(-alpha,qq,rr);
      rho_2 = rho_1;

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


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<class T_>
int _BiCG(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();
   double res;
   T_ alpha=0, beta=0, rho_1=0, rho_2=0;
   Vect<T_> r(size), z(size), zz(size), p(size), pp(size), q(size), qq(size);

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

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

   for (it=1; it<=max_it; it++) {
      P->Solve(r,z);
      P->TransSolve(rr,zz);
      rho_1 = Dot(z,rr);
      if (rho_1 == T_(0)) {
         toler = r.getNorm2() / nrm;
         if (verbose > 1)
            cout << "Convergence after " << it << " iterations." << endl;
         return it;
      }
      if (it == 1) {
         p = z;
         pp = zz;
      } else {
         beta = rho_1 / rho_2;
         p = z + beta*p;
         pp = zz + beta*pp;
      }
      A->Mult(p,q);
      A->TMult(pp,qq);
      alpha = rho_1/Dot(pp,q);
      Axpy( alpha,p,x);
      Axpy(-alpha,q,r);
      Axpy(-alpha,qq,rr);
      rho_2 = rho_1;

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

} /* namespace OFELI */

#endif

Copyright © 1998-2007 Rachid Touzani