BiCgStabl.cxx

// Copyright (C) 2003-2009 Marc Duruflé
// Copyright (C) 2001-2009 Vivien Mallet
//
// This file is part of the linear-algebra library Seldon,
// http://seldon.sourceforge.net/.
//
// Seldon is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 2.1 of the License, or (at your option)
// any later version.
//
// Seldon 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 Lesser General Public License for
// more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Seldon. If not, see http://www.gnu.org/licenses/.
 
 
#ifndef SELDON_FILE_ITERATIVE_BICGSTABL_CXX
 
namespace Seldon
{
 
 
#ifdef SELDON_WITH_VIRTUAL
  template<class T, class Vector1>
  int BiCgStabl(const VirtualMatrix<T>& A, Vector1& x, const Vector1& b,
                Preconditioner_Base<T>& M,
                Iteration<typename ClassComplexType<T>::Treal>& iter)
#else
  template <class Titer, class Matrix1, class Vector1, class Preconditioner>
  int BiCgStabl(const Matrix1& A, Vector1& x, const Vector1& b,
                Preconditioner& M, Iteration<Titer> & iter)
#endif
  {
    const int N = A.GetM();
    if (N <= 0)
      return 0;
 
    typedef typename Vector1::value_type Complexe;
    int m = iter.GetRestart();
    int l = m;
    Complexe rho_0, rho_1, alpha, beta, omega, sigma, zero, one;
    SetComplexZero(zero);
    SetComplexOne(one);
 
    // q temporary vector before preconditioning, r0 initial residual
    Vector1 q(b), r0(b);
    Vector<Complexe> gamma(l+1), gamma_prime(l+1), gamma_twice(l+1);
    Matrix<Complexe, General, RowMajor> tau(l+1,l+1);
    // history of u and residual r
    std::vector<Vector1> r(l+1, b), u(l+1, b);
    for (int i = 0; i <= l; i++)
      {
        r[i].Fill(zero);
        u[i].Fill(zero);
      }
    
    tau.Fill(zero); gamma.Fill(zero);
    gamma_prime.Fill(zero); gamma_twice.Fill(zero);
 
    // we compute the residual r = (b - Ax)
    Copy(b, r[0]);
    if (!iter.IsInitGuess_Null())
      iter.MltAdd(-one, A, x, one, r[0]);
    else
      x.Fill(zero);
 
    // we initialize iter
    int success_init = iter.Init(b);
    if (success_init !=0 )
      return iter.ErrorCode();
 
    Copy(r[0], r0); // we keep the first residual
 
    // we initialize constants
    rho_0 = one; alpha = zero; omega = one;
    tau.Fill(zero);
 
    iter.SetNumberIteration(0);
    // Loop until the stopping criteria are satisfied
    while (! iter.Finished(r[0]))
      {
        rho_0 *= -omega;
 
        // Bi-CG Part
        for (int j = 0; j < l; j++)
          {
            rho_1 = DotProd(r[j], r0);
            if (rho_0 == zero)
              {
                iter.Fail(1, "Bicgstabl breakdown #1");
                break;
              }
            beta = alpha*(rho_1/rho_0);
            rho_0 = rho_1;
            for (int i = 0; i <= j; i++)
              {
                Mlt(-beta, u[i]);
                Add(one, r[i], u[i]);
              }
            
            M.Solve(A, u[j], q); // preconditioning
            iter.Mlt(A, q, u[j+1]); // product Matrix Vector
 
            ++iter;
            sigma = DotProd(u[j+1], r0);
            if (sigma == zero)
              {
                iter.Fail(2, "Bicgstabl Breakdown #2");
                break;
              }
            alpha = rho_1/sigma;
            Add(alpha, u[0], x);
            for (int i = 0; i <= j; i++)
              Add(-alpha, u[i+1], r[i]);
 
            M.Solve(A, r[j], q); // preconditioning
            iter.Mlt(A, q, r[j+1]); // product matrix vector
 
            ++iter;
          }
 
        // MR Part  modified Gram-Schmidt
        for (int j = 1; j <= l; j++)
          {
            for (int i = 1; i < j; i++)
              {
                if (gamma(i) != zero)
                  {
                    tau(i,j) = DotProd(r[j], r[i])/gamma(i);
                    Add(-tau(i,j), r[i], r[j]);
                  }
              }
            gamma(j) = DotProd(r[j], r[j]);
            if (gamma(j) != zero)
              gamma_prime(j) = DotProd(r[0], r[j])/gamma(j);
          }
 
        // gamma = tau-1 * gamma_prime
        gamma(l) = gamma_prime(l); omega = gamma(l);
        for (int j = l-1; j >= 1; j--)
          {
            sigma = zero;
            for (int i = j+1; i <= l; i++)
              sigma += tau(j,i)*gamma(i);
 
            gamma(j) = gamma_prime(j)-sigma;
          }
 
        // gamma_twice=T*S*gamma
        for (int j = 1; j <= l-1; j++)
          {
            sigma = zero;
            for (int i = j+1; i <= l-1; i++)
              sigma += tau(j,i)*gamma(i+1);
 
            gamma_twice(j) = gamma(j+1)+sigma;
          }
 
        // update
        Add(gamma(1), r[0], x);
        Add(-gamma_prime(l), r[l], r[0]);
        Add(-gamma(l), u[l], u[0]);
        for (int j = 1;j <= l-1; j++)
          {
            Add(-gamma(j), u[j], u[0]);
            Add(gamma_twice(j), r[j], x);
            Add(-gamma_prime(j), r[j], r[0]);
          }
      }
    // change of coordinates (right preconditioning)
    Copy(x,q); M.Solve(A, q, x);
    return iter.ErrorCode();
  }
 
}
 
#define SELDON_FILE_ITERATIVE_BICGSTABL_CXX
#endif