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

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, std::ostream *p_cerr=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, &cerr);

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

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/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 Preconditioner, class Real, class Size>
       int pcg(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
               Size &max_iter, Real &tol, std::ostream *p_cerr = 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_cerr) (*p_cerr) << "[" << 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_cerr) (*p_cerr) << "[" << 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;
       }