Provided by: librheolef-dev_5.93-2_amd64 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;
       }