puzawa(4rheolef) [debian man page]

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)

pcg(4rheolef)							    rheolef-6.1 						     pcg(4rheolef)

NAME

       pcg -- conjugate gradient algorithm.

SYNOPSIS

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

EXAMPLE

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

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

DESCRIPTION

       pcg solves the symmetric positive definite linear system Ax=b using the Conjugate Gradient method.

       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

NOTE

       pcg is an iterative template routine.

       pcg follows the algorithm described on p. 15 in

       @quotation Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition, R. Barrett, M.  Berry,  T.	F.
       Chan,  J.  Demmel,  J.  Donato,	J.  Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. Van der Vorst, SIAM, 1994, ftp.netlib.org/templates/tem-
       plates.ps.  @end quotation

       The present implementation is inspired from IML++ 1.2 iterative method library, http://math.nist.gov/iml++.

IMPLEMENTATION

       template <class Matrix, class Vector, class Vector2, class Preconditioner, class Real, class Size>
       int pcg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M,
	       Size &max_iter, Real &tol, odiststream *p_derr = 0, std::string label = "cg")
       {
	   Vector b = M.solve(Mb);
	   Real norm2_b = dot(Mb,b);
	   if (norm2_b == Real(0)) norm2_b = 1;
	   Vector Mr = Mb - A*x;
	   Real  norm2_r = 0;
	   if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
	   Vector p;
	   for (Size n = 0; n <= max_iter; n++) {
	       Vector r = M.solve(Mr);
	       Real prev_norm2_r = norm2_r;
	       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;
	       }
	       if (n == 0) {
		 p = r;
	       } else {
		 Real beta = norm2_r/prev_norm2_r;
		 p = r + beta*p;
	       }
	       Vector Mq = A*p;
	       Real alpha = norm2_r/dot(Mq, p);
	       x  += alpha*p;
	       Mr -= alpha*Mq;
	   }
	   tol = ::sqrt(norm2_r/norm2_b);
	   return 1;
       }

rheolef-6.1							    rheolef-6.1 						     pcg(4rheolef)