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())
int cg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M, const solver_option &sopt=solver_option())
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::float_type Real;
std::string label = (sopt.label != "" ? sopt.label : "cg");
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;
for (sopt.n_iter = 0; sopt.n_iter <= sopt.max_iter; sopt.n_iter++) {
Real prev_norm2_r = norm2_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) {
} else {
Real
beta = norm2_r/prev_norm2_r;
}
}
return 1;
}
Float alpha[pmax+1][pmax+1]
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
DATA::size_type size_type