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
COPYRIGHT
Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito [AT] 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.