GitHub

Conjugate Gradients (CG)

Conjugate Gradients solves $Ax = b$ approximately for $x$ where $A$ is a symmetric, positive-definite linear operator and $b$ the right-hand side vector. The method uses short recurrences and therefore has fixed memory costs and fixed computational costs per iteration.

Usage

IterativeSolvers.cg — Function
cg(A, b; kwargs...) -> x, [history]

Same as cg!, but allocates a solution vector x initialized with zeros.

source
IterativeSolvers.cg! — Function
cg!(x, A, b; kwargs...) -> x, [history]

Arguments

  • x: Initial guess, will be updated in-place;
  • A: linear operator;
  • b: right-hand side.

Keywords

  • statevars::CGStateVariables: Has 3 arrays similar to x to hold intermediate results;
  • initially_zero::Bool: If true assumes that iszero(x) so that one matrix-vector product can be saved when computing the initial residual vector;
  • Pl = Identity(): left preconditioner of the method. Should be symmetric, positive-definite like A;
  • abstol::Real = zero(real(eltype(b))), reltol::Real = sqrt(eps(real(eltype(b)))): absolute and relative tolerance for the stopping condition |r_k| ≤ max(reltol * |r_0|, abstol), where r_k ≈ A * x_k - b is approximately the residual in the kth iteration.
    Note

    The true residual norm is never explicitly computed during the iterations for performance reasons; it may accumulate rounding errors.

  • maxiter::Int = size(A,2): maximum number of iterations;
  • verbose::Bool = false: print method information;
  • log::Bool = false: keep track of the residual norm in each iteration.

Output

if log is false

  • x: approximated solution.

if log is true

  • x: approximated solution.
  • ch: convergence history.

ConvergenceHistory keys

  • :tol => ::Real: stopping tolerance.
  • :resnom => ::Vector: residual norm at each iteration.
source

On the GPU

The method should work fine on the GPU. As a minimal working example, consider:

using LinearAlgebra, CuArrays, IterativeSolvers

n = 100
A = cu(rand(n, n))
A = A + A' + 2*n*I
b = cu(rand(n))
x = cg(A, b)
Note

Make sure that all state vectors are stored on the GPU. For instance when calling cg!(x, A, b), one might have an issue when x is stored on the GPU, while b is stored on the CPU – IterativeSolvers.jl does not copy the vectors to the same device.

Implementation details

The current implementation follows a rather standard approach. Note that preconditioned CG (or PCG) is slightly different from ordinary CG, because the former must compute the residual explicitly, while it is available as byproduct in the latter. Our implementation of CG ensures the minimal number of vector operations.

Tip

CG can be used as an iterator.