Conjugate Gradients (CG)
Usage
#
IterativeSolvers.cg
— Function
cg(A, b; kwargs...) -> x, [history]
Same as cg!
, but allocates a solution vector x
initialized with zeros.
#
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 tox
to hold intermediate results; -
initially_zero::Bool
: Iftrue
assumes thatiszero(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 likeA
; -
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)
, wherer_k ≈ A * x_k - b
is approximately the residual in thek
th iteration.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.
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)
Make sure that all state vectors are stored on the GPU. For instance when calling |
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.
CG can be used as an iterator. |