puzawa(4rheolef)						    rheolef-6.1 						  puzawa(4rheolef)

NAME

       puzawa -- Uzawa algorithm.

SYNOPSIS

	   template <class Matrix, class Vector, class Preconditioner, class Real>
	   int puzawa (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
	     int &max_iter, Real &tol, const Real& rho, odiststream *p_derr=0);

EXAMPLE

       The simplest call to 'puzawa' has the folling form:

	   size_t max_iter = 100;
	   double tol = 1e-7;
	   int status = puzawa(A, x, b, EYE, max_iter, tol, 1.0, &derr);

DESCRIPTION

       puzawa solves the linear system A*x=b using the Uzawa method. The Uzawa method is a descent method in the direction opposite to the  gradi-
       ent, with a constant step length 'rho'. The convergence is assured when the step length 'rho' is small enough.  If matrix  A  is  symmetric
       positive definite, please uses 'pcg' that computes automatically the optimal descdnt step length at each iteration.

       The return value indicates convergence within max_iter (input) iterations(0), or no convergence within max_iter iterations(1).  Upon suc-
       cessful 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

IMPLEMENTATION

       template < class Matrix, class Vector, class Preconditioner, class Real, class Size>
       int puzawa(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
	   Size &max_iter, Real &tol, const Real& rho,
	   odiststream *p_derr, std::string label)
       {
	   Vector b = M.solve(Mb);
	   Real norm2_b = dot(Mb,b);
	   Real norm2_r = norm2_b;
	   if (norm2_b == Real(0)) norm2_b = 1;
	   if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
	   for (Size n = 0; n <= max_iter; n++) {
	       Vector Mr = A*x - Mb;
	       Vector r = M.solve(Mr);
	       norm2_r = dot(Mr, r);
	       if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << sqrt(norm2_r/norm2_b) << std::endl;
	       if (norm2_r <= sqr(tol)*norm2_b) {
		 tol = sqrt(norm2_r/norm2_b);
		 max_iter = n;
		 return 0;
	       }
	       x  -= rho*r;
	   }
	   tol = sqrt(norm2_r/norm2_b);
	   return 1;
       }

rheolef-6.1							    rheolef-6.1 						  puzawa(4rheolef)