pgmres(4rheolef)						    rheolef-6.1 						  pgmres(4rheolef)

NAME

       pgmres -- generalized minimum residual method

SYNOPSIS

	  template <class Matrix, class Vector, class Preconditioner,
	    class SmallMatrix, class SmallVector, class Real, class Size>
	  int pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
	    SmallMatrix &H, const SmallVector& dummy,
	    Size m, Size &max_iter, Real &tol);

EXAMPLE

       The simplest call to pgmres has the folling form:

	       double tol = 1e-7;
	       size_t max_iter = 100;
	       size_t m = 6;
	       boost::numeric::ublas::matrix<double> H(m+1,m+1);
	       vec<double,sequential> dummy;
	       int status = pgmres (a, x, b, ic0(a), H, dummy, m, max_iter, tol);

DESCRIPTION

       pgmres solves the unsymmetric linear system Ax = b using the generalized minimum residual 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 In addition, M specifies a preconditioner, H specifies a matrix to hold the  coefficients	of
	      the upper Hessenberg matrix constructed by the pgmres iterations, m specifies the number of iterations for each restart.

       pgmres  requires  two  matrices as input, A and H.  The matrix A, which will typically be a sparse matrix) corresponds to the matrix in the
       linear system Ax=b.  The matrix H, which will be typically a dense matrix, corresponds to the upper Hessenberg matrix H that is constructed
       during  the  pgmres iterations. Within pgmres, H is used in a different way than A, so its class must supply different functionality.  That
       is, A is only accessed though its matrix-vector and transpose-matrix-vector multiplication functions.  On the other hand, pgmres  solves  a
       dense  upper  triangular  linear system of equations on H. Therefore, the class to which H belongs must provide H(i,j) operator for element
       acess.

NOTE

       It is important to remember that we use the convention that indices are 0-based. That is H(0,0) is the first component  of  the	matrix	H.
       Also,  the type of the matrix must be compatible with the type of single vector entry. That is, operations such as H(i,j)*x(j) must be able
       to be carried out.

       pgmres is an iterative template routine.

       pgmres follows the algorithm described on p. 20 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/templates.ps.  @end quotation

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

IMPLEMENTATION

       template <class SmallMatrix, class Vector, class SmallVector, class Size>
       void
       Update(Vector &x, Size k, SmallMatrix &h, SmallVector &s, Vector v[])
       {
	 SmallVector y = s;
	 // Back solve:
	 for (int i = k; i >= 0; i--) {
	   y(i) /= h(i,i);
	   for (int j = i - 1; j >= 0; j--)
	     y(j) -= h(j,i) * y(i);
	 }
	 for (Size j = 0; j <= k; j++) {
	   x += v[j] * y(j);
	 }
       }
       template<class Real>
       void GeneratePlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
       {
	 if (dy == Real(0)) {
	   cs = 1.0;
	   sn = 0.0;
	 } else if (abs(dy) > abs(dx)) {
	   Real temp = dx / dy;
	   sn = 1.0 / ::sqrt( 1.0 + temp*temp );
	   cs = temp * sn;
	 } else {
	   Real temp = dy / dx;
	   cs = 1.0 / ::sqrt( 1.0 + temp*temp );
	   sn = temp * cs;
	 }
       }
       template<class Real>
       void ApplyPlaneRotation(Real &dx, Real &dy, Real &cs, Real &sn)
       {
	 Real temp  =  cs * dx + sn * dy;
	 dy = -sn * dx + cs * dy;
	 dx = temp;
       }
       template <class Matrix, class Vector, class Preconditioner,
		 class SmallMatrix, class SmallVector, class Real, class Size>
       int
       pgmres (const Matrix &A, Vector &x, const Vector &b, const Preconditioner &M,
	     SmallMatrix &H, const SmallVector&, const Size &m, Size &max_iter, Real &tol)
       {
	 odiststream* p_derr = &derr;
	 std::string label = "pgmres";
	 if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
	 Vector w;
	 SmallVector s(m+1), cs(m+1), sn(m+1);
	 Size i;
	 Size j = 1;
	 Size k;
	 Real resid;
	 Real normb = norm(M.solve(b));
	 Vector r = M.solve(b - A * x);
	 Real beta = norm(r);
	 if (normb == Real(0)) {
	   normb = 1;
	 }
	 resid = norm(r) / normb;
	 if (resid  <= tol) {
	   tol = resid;
	   max_iter = 0;
	   return 0;
	 }
	 Vector *v = new Vector[m+1];
	 while (j <= max_iter) {
	   v[0] = r * (1.0 / beta);    // ??? r / beta
	   s = 0.0;
	   s(0) = beta;
	   for (i = 0; i < m && j <= max_iter; i++, j++) {
	     w = M.solve(A * v[i]);
	     for (k = 0; k <= i; k++) {
	       H(k, i) = dot(w, v[k]);
	       w -= H(k, i) * v[k];
	     }
	     H(i+1, i) = norm(w);
	     v[i+1] = w * (1.0 / H(i+1, i)); // ??? w / H(i+1, i)
	     for (k = 0; k < i; k++) {
	       ApplyPlaneRotation(H(k,i), H(k+1,i), cs(k), sn(k));
	     }
	     GeneratePlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
	     ApplyPlaneRotation(H(i,i), H(i+1,i), cs(i), sn(i));
	     ApplyPlaneRotation(s(i), s(i+1), cs(i), sn(i));
	     resid = abs(s(i+1)) / normb;
	     if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
	     if (resid	< tol) {
	       Update(x, i, H, s, v);
	       tol = resid;
	       max_iter = j;
	       delete [] v;
	       return 0;
	     }
	   }
	   Update(x, m - 1, H, s, v);
	   r = M.solve(b - A * x);
	   beta = norm(r);
	   resid = beta / normb;
	   if (p_derr) (*p_derr) << "[" << label << "] " << j << " " << resid << std::endl;
	   if (resid < tol) {
	     tol = resid;
	     max_iter = j;
	     delete [] v;
	     return 0;
	   }
	 }
	 tol = resid;
	 delete [] v;
	 return 1;
       }

rheolef-6.1							    rheolef-6.1 						  pgmres(4rheolef)