#ifndef __RICHARDSON_H
#define __RICHARDSON_H

// Iterative template routine -- preconditioned Richardson
//
// RICHARDSON solves the linear system Ax=b using
// the preconditioned richardson iteration.
// The returned value indicates convergence within
// max_iter iterations (return value 0)
// or no convergence within max_iter iterations (return value 1)
// Upon successful return (0), the output arguments have the
// following values:
//        x: computed solution
// mat_iter: number of iterations to satisfy the stopping criterion
//      tol: residual after the final iteration

template <class MATRIX, class VECTOR, class PRECONDITIONER, class REAL>
int
RICHARDSON (const MATRIX & A, VECTOR & x, const VECTOR & b,
	    const PRECONDITIONER & M, int & max_iter, REAL & tol)
{
  REAL resid;
  VECTOR z(b.size ());
  REAL normb = norm (b);
  VECTOR r = b - A * x;

  if (normb == 0.0) normb = 1;
  resid = norm (r) / normb;
  
  if (resid <= tol)
    {
      tol = resid;
      max_iter = 0;
      return 0;
    }
  
  for (int i=1; i<max_iter; i++)
    {
      z = M.solve (r);
      x += z;
      r = b - A * x;
      resid = norm(r) / normb;
      
      if (resid <= tol)
	{
	  tol = resid;
	  max_iter = i;
	  return 0;
	}
    }
  
  tol = resid;
  return 1;
}



#endif // __RICHARDSON_H