.\"
.de Id
..
.de Sp
.if n .sp
.if t .sp 0.4
..
.TH pcg 4rheolef "rheolef-6.7" "rheolef-6.7" "rheolef-6.7"
.\" label: /*D:pcg
.SH NAME
\fBpcg\fP -- conjugate gradient algorithm.
.\" skip: @findex pcg
.\" skip: @cindex conjugate gradient algorithm
.\" skip: @cindex iterative solver
.\" skip: @cindex preconditioner
.SH SYNOPSIS

.\" begin_example
.Sp
.nf
    template <class Matrix, class Vector, class Preconditioner, class Real>
    int pcg (const Matrix &A, Vector &x, const Vector &b,
      const Preconditioner &M, int &max_iter, Real &tol, odiststream *p_derr=0);
.Sp
.fi
.\" end_example
.PP
.SH EXAMPLE

The simplest call to 'pcg' has the folling form:
.\" begin_example
.Sp
.nf
    size_t max_iter = 100;
    double tol = 1e-7;
    int status = pcg(a, x, b, EYE, max_iter, tol, &derr);
.Sp
.fi
.\" end_example
.PP
.SH DESCRIPTION

\fBpcg\fP solves the symmetric positive definite linear
system Ax=b using the Conjugate Gradient method.
.PP
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:
.\" begin table
.\" start item
.TP
.B x
approximate solution to Ax = b

.\" start item
.TP
.B max_iter
the number of iterations performed before the tolerance was reached

.\" start item
.TP
.B tol
the residual after the final iteration
.\" end table
.PP
.SH NOTE

.PP
\fBpcg\fP is an iterative template routine.
.PP
\fBpcg\fP follows the algorithm described on p. 15 in
.PP
\fITemplates for the solution of linear systems: building blocks for iterative methods\fP,
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,
\fBftp.netlib.org/templates/templates.ps\fP.
.PP
The present implementation is inspired from
\fBIML++ 1.2\fP iterative method library,
\fBhttp://math.nist.gov/iml++\fP.
.PP
.\" skip start:AUTHOR:

.\" skip start:DATE:

.\" skip start:METHODS:
.\" END
.SH IMPLEMENTATION
.\" begin_example
.Sp
.nf
template <class Matrix, class Vector, class Vector2, class Preconditioner, class Real, class Size>
int pcg(const Matrix &A, Vector &x, const Vector2 &Mb, const Preconditioner &M,
        Size &max_iter, Real &tol, odiststream *p_derr = 0, std::string label = "cg")
{
    Vector b = M.solve(Mb);
    Real norm2_b = dot(Mb,b);
    if (norm2_b == Real(0)) norm2_b = 1;
    Vector Mr = Mb - A*x;
    Real  norm2_r = 0;
    if (p_derr) (*p_derr) << "[" << label << "] #iteration residue" << std::endl;
    Vector p;
    for (Size n = 0; n <= max_iter; n++) {
        Vector r = M.solve(Mr);
        Real prev_norm2_r = norm2_r;
        norm2_r = dot(Mr, r);
        if (p_derr) (*p_derr) << "[" << label << "] " << n << " " << sqrt(norm2_r/norm2_b) << std::endl;
        if (norm2_r <= sqr(tol)*norm2_b) {
          tol = sqrt(norm2_r/norm2_b);
          max_iter = n;
          return 0;
        }
        if (n == 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;
    }
    tol = sqrt(norm2_r/norm2_b);
    return 1;
}
.Sp
.fi
.\" end_example