Provided by: librheolef-dev_5.93-2_i386

NAME

```       pcg -- conjugate gradient algorithm.

```

SYNOPSIS

```           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, std::ostream *p_cerr=0);

```

EXAMPLE

```       The simplest call to 'pcg' has the folling form:

size_t max_iter = 100;
double tol = 1e-7;
int status = pcg(a, x, b, EYE, max_iter, tol, &cerr);

```

DESCRIPTION

```       pcg solves the symmetric positive definite linear system Ax=b using the

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

max_iter
the number of iterations  performed  before  the  tolerance  was
reached

tol    the residual after the final iteration

```

NOTE

```       pcg is an iterative template routine.

pcg follows the algorithm described on p. 15 in

@quotation  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.  @end quotation

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

```

IMPLEMENTATION

```       template < class Matrix, class Vector, class Preconditioner, class Real, class Size>
int pcg(const Matrix &A, Vector &x, const Vector &Mb, const Preconditioner &M,
Size &max_iter, Real &tol, std::ostream *p_cerr = 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_cerr) (*p_cerr) << "[" << 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_cerr) (*p_cerr) << "[" << 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;
}
```