`symmlq`

Solves a symmetric linear system with the symmetric LQ method.

Syntax

x = symmlq(A, b)

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

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

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

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

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

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

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

Description

symmlq is meant for symmetric A, definite or indefinite. It uses the Lanczos process, so memory stays constant, and it solves a related tridiagonal system by LQ factorization. It is a companion of minres, which minimizes the residual instead. A preconditioner must be symmetric positive definite.

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

Symmetric 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] = symmlq(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

symmlq