.TH "cg" 5rheolef "Version 7.2" "rheolef" \" -*- nroff -*-
.ad l
.nh
.SH NAME
cg \- conjugate gradient algorithm (rheolef-7\&.2)
.PP
  
.SH "SYNOPSIS"
.PP
.PP
.nf
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())
.fi
.PP
.SH "EXAMPLE"
.PP
.PP
.nf
    solver_option sopt;
    sopt\&.max_iter = 100;
    sopt\&.tol = 1e-7;
    int status = cg (A, x, b, eye(), sopt);
.fi
.PP
 
.SH "DESCRIPTION"
.PP
This function solves the \fIsymmetric positive definite\fP linear system \fCA*x=b\fP with the conjugate gradient method\&. The \fCcg\fP function follows the algorithm described on p\&. 15 in 
.PP
.nf
    R\&. Barrett, M\&. Berry, T\&. F\&. Chan, J\&. Demmel, J\&. Donato,
    J\&. Dongarra, V\&. Eijkhout, R\&. Pozo, C\&. Romine and H\&. Van der Vorst,
    Templates for the solution of linear systems: building blocks for iterative methods, 
    SIAM, 1994\&.

.fi
.PP
 The fourth argument of \fCcg\fP is a preconditionner: here, the \fBeye(5)\fP one is a do-nothing preconditionner, for simplicity\&. Finally, the \fBsolver_option(4)\fP variable \fCsopt\fP transmits the stopping criterion with \fCsopt\&.tol\fP and \fCsopt\&.max_iter\fP\&.
.PP
On return, the \fCsopt\&.residue\fP and \fCsopt\&.n_iter\fP indicate the reached residue and the number of iterations effectively performed\&. The return \fCstatus\fP is zero when the prescribed tolerance \fCtol\fP has been obtained, and non-zero otherwise\&. Also, the \fCx\fP variable contains the approximate solution\&. See also the \fBsolver_option(4)\fP for more controls upon the stopping criterion\&.
.SH "IMPLEMENTATION"
.PP
This documentation has been generated from file linalg/lib/cg\&.h
.PP
The present template implementation is inspired from the \fCIML++ 1\&.2\fP iterative method library, http://math.nist.gov/iml++
.PP
.PP
.nf
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())
.fi
.PP
 
.PP
.nf
{
  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;
}

.fi
.PP
 
.SH AUTHOR
Pierre  Saramito  <Pierre.Saramito@imag.fr>
.SH 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.