Provided by: librheolef-dev_5.93-2_i386 bug

NAME

       gmres -- generalized minimum residual method

SYNOPSIS

               template <class Operator, class Vector, class Preconditioner,
                   class Matrix, class Real, class Int>
               int gmres (const Operator &A, Vector &x, const Vector &b,
                   const Preconditioner &M, Matrix &H, Int m, Int &max_iter, Real &tol);

EXAMPLE

       The simplest call to gmres has the folling form:

               int m = 6;
               dns H(m+1,m+1);
               int status = gmres(a, x, b, EYE, H, m, 100, 1e-7);

DESCRIPTION

       gmres 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
       successful 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  gmres
              iterations, m  specifies  the  number  of  iterations  for  each
              restart.

       gmres  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 gmres iterations. Within gmres, 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,  gmres
       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.

       gmres is an iterative template routine.

       gmres  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 Matrix, class Vector, class Int >
       void
       Update(Vector &x, Int k, Matrix &h, Vector &s, Vector v[])
       {
         Vector y(s);

         // Backsolve:
         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 (Int j = 0; j <= k; j++)
           x += v[j] * y(j);
       }

       #ifdef TO_CLEAN
       template < class Real >
       Real
       abs(Real x)
       {
         return (x > Real(0) ? x : -x);
       }
       #endif // TO_CLEAN

       template < class Operator, class Vector, class Preconditioner,
                  class Matrix, class Real, class Int >
       int
       gmres(const Operator &A, Vector &x, const Vector &b,
             const Preconditioner &M, Matrix &H, const Int &m, Int &max_iter,
             Real &tol)
       {
         Real resid;
         Int i, j = 1, k;
         Vector s(m+1), cs(m+1), sn(m+1), w;

         Real normb = norm(M.solve(b));
         Vector r = M.solve(b - A * x);
         Real beta = norm(r);

         if (normb == Real(0))
           normb = 1;

         if ((resid = norm(r) / normb) <= 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));

             if ((resid = abs(s(i+1)) / normb) < 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);
           if ((resid = beta / normb) < tol) {
             tol = resid;
             max_iter = j;
             delete [] v;
             return 0;
           }
         }

         tol = resid;
         delete [] v;
         return 1;
       }
       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;
       }