`pcg`

Solves a symmetric positive definite linear system with preconditioned conjugate gradients.

Syntax

x = pcg(A, b)

Solves A*x = b with default tolerance and iteration limit.

x = pcg(A, b, tol, maxit)

Sets the relative residual tolerance and the maximum number of iterations.

x = pcg(A, b, tol, maxit, M1, M2, x0)

Adds a preconditioner M = M1*M2 and an initial guess x0.

[x, flag, relres, iter, resvec] = pcg(___)

Also returns the convergence flag, final relative residual, iteration count and residual history.

Description

pcg is the standard method when A is symmetric and positive definite. It needs one product with A per iteration, stores only a few vectors, and in exact arithmetic finds the answer in at most size(A,1) steps. In practice the preconditioner decides how fast it gets there.

If A is not symmetric positive definite, the method can break down or give a wrong answer; use gmres or bicgstab instead.

When called with fewer than two outputs it prints a one-line message about whether it converged. With a flag output it stays silent, so check flag yourself. A good preconditioner such as an incomplete factorization (ilu, ichol) can cut the iteration count a lot.

Input arguments

A: Square coefficient matrix (usually sparse), or a function handle that returns A*x.
b: Right-hand side column vector with size(A,1) entries.
tol: Relative residual tolerance; the default is 1e-6.
maxit: Maximum number of iterations.
M1, M2: Preconditioner factors, as matrices or function handles; the solver works with M1*M2.
x0: Initial guess. The default is a vector of zeros.

Output arguments

x: The computed solution, a column vector.
flag: 0 if it converged within maxit iterations, 1 if it ran out of iterations, 2 if the preconditioner was ill-conditioned, 3 if it stagnated, and other values for a breakdown.
relres: Relative residual norm(b - A*x)/norm(b) at the returned iterate.
iter: Iteration at which x was computed.
resvec: History of residual norms.

Example

Tridiagonal system, plain and preconditioned

A = sparse(diag(4*ones(5,1)) + diag(ones(4,1),1) + diag(ones(4,1),-1));
b = (1:5)';
[x, flag, relres, iter] = pcg(A, b, 1e-10, 100);
fprintf('flag = %d, iterations = %d\n', flag, iter)
disp(x')
L = ichol(A);
[x2, flag2] = pcg(A, b, 1e-10, 100, L, L');
fprintf('preconditioned flag = %d, same answer: %d\n', flag2, norm(x - x2) < 1e-8)

Output

flag = 0, iterations = 5
    0.1679    0.3282    0.5192    0.5949    1.1013
preconditioned flag = 0, same answer: 1

pcg