`cgs`

Solves a square linear system with the conjugate gradients squared method.

Syntax

x = cgs(A, b)

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

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

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

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

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

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

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

Description

cgs squares the residual polynomial of biconjugate gradients, which can double the convergence speed when the method behaves well but makes the residual curve irregular and sensitive to rounding. It needs two products with A per iteration and none with A'.

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

A = sparse(diag(4*ones(5,1)) + diag(ones(4,1),1) + diag(ones(4,1),-1));
b = (1:5)';
[x, flag, relres, iter] = cgs(A, b, 1e-10, 100);
fprintf('flag = %d, relres < 1e-10: %d\n', flag, relres < 1e-10)
disp(x')

Output

flag = 0, relres < 1e-10: 1
    0.1679    0.3282    0.5192    0.5949    1.1013

cgs