---
name: "pcg"
section: "4rheolef"
platform: "linux"
package: "librheolef-dev"
version: "6.7-1+b4"
release: "Debian 9 (Stretch)"
platforms: ["linux"]
url: "https://new.carta.tech/pcg"
---

# pcg(4rheolef)

## 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, odiststream *p_derr=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, &derr);

```

## DESCRIPTION

**pcg** 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

- **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

*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, 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;
}

```