cPCG (version 1.0)

pcgsolve: Preconditioned conjugate gradient method

Description

Preconditioned conjugate gradient method for solving system of linear equations Ax = b, where A is symmetric and positive definite, b is a column vector.

Usage

pcgsolve(A, b, preconditioner = "Jacobi", tol = 1e-6, maxIter = 1000)

Arguments

A

matrix, symmetric and positive definite.

b

vector, with same dimension as number of rows of A.

preconditioner

string, method for preconditioning: "Jacobi" (default), "SSOR", or "ICC".

tol

numeric, threshold for convergence, default is 1e-6.

maxIter

numeric, maximum iteration, default is 1000.

Value

Returns a vector representing solution x.

Preconditioners

Jacobi: The Jacobi preconditioner is the diagonal of the matrix A, with an assumption that all diagonal elements are non-zero.

SSOR: The symmetric successive over-relaxation preconditioner, implemented as \(M = (D+L) D^{-1} (D+L)^T\). [1]

ICC: The incomplete Cholesky factorization preconditioner. [2]

Warning

Users need to check that input matrix A is symmetric and positive definite before applying the function.

Details

When the condition number for \(A\) is large, the conjugate gradient (CG) method may fail to converge in a reasonable number of iterations. The Preconditioned Conjugate Gradient (PCG) Method applies a precondition matrix \(C\) and approaches the problem by solving: $${C}^{-1} A x = {C}^{-1} b$$ where the symmetric and positive-definite matrix \(C\) approximates \(A\) and \({C}^{-1} A \) improves the condition number of \(A\).

Common choices for the preconditioner include: Jacobi preconditioning, symmetric successive over-relaxation (SSOR), and incomplete Cholesky factorization [2].

References

[1] David Young. <U+201C>Iterative methods for solving partial difference equations of elliptic type<U+201D>. In: Transactions of the American Mathematical Society 76.1 (1954), pp. 92<U+2013>111.

[2] David S Kershaw. <U+201C>The incomplete Cholesky<U+2014>conjugate gradient method for the iter- ative solution of systems of linear equations<U+201D>. In: Journal of computational physics 26.1 (1978), pp. 43<U+2013>65.

See Also

cgsolve

Examples

# NOT RUN {
test_A <- matrix(c(4,1,1,3), ncol = 2)
test_b <- matrix(1:2, ncol = 1)
pcgsolve(test_A, test_b, "ICC")
# }