#ifndef __RICHARDSON_H
#define __RICHARDSON_H

// Iterative template routine -- preconditioned Richardson
//
// RICHARDSON solves the linear system Ax=b using
// the Richardson iteration with damping paramter tau.
// 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 REAL>
int
RICHARDSON (const MATRIX & A, VECTOR & x, const VECTOR & b,
	    REAL tau, int & max_iter, REAL & tol)
{
  assert (x.size() == b.size());

  REAL resid;
  VECTOR z(b.size ());
  REAL normb = b.l2norm ();
  VECTOR r = b - A * x;

  if (normb == 0.0) normb = 1;
  resid = r.l2norm() / normb;

  // IF YOU WANT OUTPUT
  cout << "RICHARDSON " << 0 << ": " << resid; cout.flush ();
  
  if (resid <= tol)
    {
      tol = resid;
      max_iter = 0;
      return 0;
    }
  
  for (int i=1; i<max_iter; i++)
    {
      z = tau * r;
      x += z;
      r = b - A * x;
      resid = r.l2norm() / normb;
      
      // IF YOU WANT OUTPUT
      cout << "\r" << "RICHARDSON " << i << ": " << resid << "      "; cout.flush ();
      
      if (resid <= tol)
	{
	  tol = resid;
	  max_iter = i;
	  cout << endl;
	  return 0;
	}
    }
  
  tol = resid;
  cout << endl;
  return 1;
}



#endif // __RICHARDSON_H