---
name: "cg"
section: "4rheolef"
platform: "linux"
package: "librheolef-dev"
version: "7.0-2+b1"
release: "Debian 10 (Buster)"
platforms: ["linux"]
url: "https://new.carta.tech/cg/4rheolef"
---

# cg(4rheolef)

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

solver\_option(2)

## 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.