GNU.WIKI: The GNU/Linux Knowledge Base

  [HOME] [PHP Manual] [HowTo] [ABS] [MAN1] [MAN2] [MAN3] [MAN4] [MAN5] [MAN6] [MAN7] [MAN8] [MAN9]

  [0-9] [Aa] [Bb] [Cc] [Dd] [Ee] [Ff] [Gg] [Hh] [Ii] [Jj] [Kk] [Ll] [Mm] [Nn] [Oo] [Pp] [Qq] [Rr] [Ss] [Tt] [Uu] [Vv] [Ww] [Xx] [Yy] [Zz]


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
       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  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 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.

       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;
       }



  All copyrights belong to their respective owners. Other content (c) 2014-2018, GNU.WIKI. Please report site errors to webmaster@gnu.wiki.
Page load time: 0.128 seconds. Last modified: November 04 2018 12:49:43.