`cg`

conjugate gradient algorithm

# Synopsis

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())
rheolef::cg
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
Definition: cg.h:94
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416

# Example

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

# Description

This function solves the symmetric positive definite linear system `A*x=b` with the conjugate gradient method. The `cg` function follows the algorithm described on p. 15 in

        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.

The fourth argument of `cg` is a preconditionner: here, the `eye` one is a do-nothing preconditionner, for simplicity. Finally, the `solver_option` variable `sopt` transmits the stopping criterion with `sopt.tol` and `sopt.max_iter`.

On return, the `sopt.residue` and `sopt.n_iter` indicate the reached residue and the number of iterations effectively performed. The return `status` is zero when the prescribed tolerance `tol` has been obtained, and non-zero otherwise. Also, the `x` variable contains the approximate solution. See also the `solver_option` for more controls upon the stopping criterion.

# Implementation

This documentation has been generated from file linalg/lib/cg.h

The present template implementation is inspired from the `IML++ 1.2` iterative method library, http://math.nist.gov/iml++

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;
}
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
p
Definition: sphere.icc:25
rheolef::Vector::size_type
DATA::size_type size_type
Definition: Vector.h:188