NAME

       cg -- conjugate gradient algorithm.

SYNOPSIS

         template <class Matrix, class Vector, class Preconditioner, class Real>
         int cg (const Matrix &A, Vector &x, const Vector &b,
           const solver_option& sopt = solver_option())

EXAMPLE

       The simplest call to cg has the folling form:

           solver_option sopt;
           sopt.max_iter = 100;
           sopt.tol = 1e-7;
           int status = cg(a, x, b, eye(), sopt);

DESCRIPTION

       cg  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

       sopt.n_iter
              the number of iterations performed before the tolerance was reached

       sopt.residue
              the residual after the final iteration See also the solver_option(2).

NOTE

       cg is an iterative template routine.

       cg follows the algorithm described on p. 15 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 Matrix, class Vector, class Vector2, class Preconditioner>
       int cg (const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M,
               const solver_option& sopt = solver_option())
       {
         typedef typename Vector::size_type  Size;
         typedef typename Vector::float_type Real;
         std::string label = (sopt.label != "" ? sopt.label : "cg");
         Vector b = M.solve(Mb);
         Real norm2_b = dot(Mb,b);
         if (sopt.absolute_stopping || norm2_b == Real(0)) norm2_b = 1;
         Vector Mr = Mb - A*x;
         Real  norm2_r = 0;
         if (sopt.p_err) (*sopt.p_err) << "[" << label << "] #iteration residue" << std::endl;
         Vector p;
         for (sopt.n_iter = 0; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
           Vector r = M.solve(Mr);
           Real prev_norm2_r = norm2_r;
           norm2_r = dot(Mr, r);
           sopt.residue = sqrt(norm2_r/norm2_b);
           if (sopt.p_err) (*sopt.p_err) << "[" << label << "] " << sopt.n_iter << " " << sopt.residue << std::endl;
           if (sopt.residue <= sopt.tol) return 0;
           if (sopt.n_iter == 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;
         }
         return 1;
       }

SEE ALSO

       solver_option(2)

COPYRIGHT

       Copyright  (C) 2000-2018 Pierre Saramito <Pierre.Saramito@imag.fr> GPLv3+: GNU GPL version
       3 or later <http://gnu.org/licenses/gpl.html>.  This is free software:  you  are  free  to
       change and redistribute it.  There is NO WARRANTY, to the extent permitted by law.