---
name: "cg"
section: "5rheolef"
platform: "linux"
package: "librheolef-dev"
version: "7.2-4"
release: "Debian 13 (Trixie)"
platforms: ["linux"]
url: "https://new.carta.tech/cg/5rheolef"
---

# cg(5rheolef)

conjugate gradient algorithm (rheolef-7.2)

## NAME

cg - conjugate gradient algorithm (rheolef-7.2)

## 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())
```

## 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(5)** one is a do-nothing preconditionner, for simplicity. Finally, the **solver\_option(4)** 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(4)** 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::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;
}
```

## AUTHOR

Pierre  Saramito  \<Pierre.Saramito@imag.fr\>

## COPYRIGHT

Copyright   (C)  2000-2018  Pierre  Saramito  \<Pierre.Saramito@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.