/***************************************************************************
    templates_bicg.h, modified by Joris Koster, 2002, jhhkoste@math.uu.nl
 ***************************************************************************/

//*****************************************************************
// Iterative template routine -- BiCG
//
// BiCG solves the unsymmetric linear system Ax = b
// using the Preconditioned BiConjugate Gradient method
//
// BiCG follows the algorithm described on p. 22 of the
// SIAM Templates book.
//
// The return value indicates convergence within max_iter (input)
// iterations (0), or no convergence within max_iter iterations (1).
//
// Upon successful return, output arguments have the following values:
//
//        x  --  approximate solution to Ax = b
// max_iter  --  the number of iterations performed before the
//               tolerance was reached
//      tol  --  the residual after the final iteration
//
//*****************************************************************

#ifndef __TEMPLATES_BICG_H
#define __TEMPLATES_BICG_H


template < class Matrix, class Vector, class Preconditioner, class Real >
int
BiCG(const Matrix &A, Vector &x, const Vector &b,
     const Preconditioner &M, int &max_iter, Real &tol)
{
  Real resid;
  Real rho_1, rho_2, alpha, beta;

  Vector    z(x.M()),
            ztilde(x.M()),
            p(x.M()),
            ptilde(x.M()),
            q(x.M()),
            qtilde(x.M()),
            r((x.M())),
            rtilde((x.M()));

  Real normb = norm(b);
  rtilde = A *x;
  r = b - rtilde;
  rtilde = r;

  if (normb == 0.0)
    normb = 1;

  if ((resid = norm(r) / normb) <= tol) {
    tol = resid;
    max_iter = 0;
    return 0;
  }

  cprintf( "0: %.16e\n", resid);
  for (int i = 1; i <= max_iter; i++) {
    z = M.solve(r);
    ztilde = M.trans_solve(rtilde);
    rho_1 = dot(z, rtilde);
    if (rho_1 == 0) {
      tol = norm(r) / normb;
      max_iter = i;
      return 2;
    }
    if (i == 1) {
      p = z;
      ptilde = ztilde;
    } else {
      beta = rho_1 / rho_2;
      p = z + beta * p;
      ptilde = ztilde + beta * ptilde;
    }
    q = A * p;
    qtilde = A.trans_mult(ptilde);
    alpha = rho_1 / dot(ptilde, q);
    x += alpha * p;
    r -= alpha * q;
    rtilde -= alpha * qtilde;

    rho_2 = rho_1;
    if ((resid = norm(r) / normb) < tol) {
      cprintf( "** %d: %.16e\n", i, resid);
      tol = resid;
      max_iter = i;
      return 0;
    }
    cprintf( "%d: %.16e\n", i, resid);
  }

  tol = resid;
  return 1;
}

#endif // end ifndef __TEMPLATES_BICG_H