`qmr`

Solves a square linear system with the quasi-minimal residual method.

Syntax

x = qmr(A, b)

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

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

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

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

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

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

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

Description

qmr is a variant of biconjugate gradients that minimizes a quantity close to the residual, which smooths the convergence curve that bicg produces. Like bicg it needs a product with A' as well as with A, and it accepts a left and a right preconditioner.

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] = qmr(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

qmr