symmlq

Symmetric LQ method

Syntax

x = symmlq(A,b)
symmlq(A,b,tol)
symmlq(A,b,tol,maxit)
symmlq(A,b,tol,maxit,M)
symmlq(A,b,tol,maxit,M1,M2)
symmlq(A,b,tol,maxit,M1,M2,x0)
[x,flag] = symmlq(A,b,...)
[x,flag,relres] = symmlq(A,b,...)
[x,flag,relres,iter] = symmlq(A,b,...)
[x,flag,relres,iter,resvec] = symmlq(A,b,...)
[x,flag,relres,iter,resvec,resveccg] = symmlq(A,b,...)


Description

x = symmlq(A,b) attempts to solve the system of linear equations A*x=b for x. The n-by-n coefficient matrix A must be symmetric but need not be positive definite. It should also be large and sparse. The column vector b must have length n. You can specify A as a function handle, afun, such that afun(x) returns A*x.

Parameterizing Functions explains how to provide additional parameters to the function afun, as well as the preconditioner function mfun described below, if necessary.

If symmlq converges, a message to that effect is displayed. If symmlq fails to converge after the maximum number of iterations or halts for any reason, a warning message is printed displaying the relative residual norm(b-A*x)/norm(b) and the iteration number at which the method stopped or failed.

symmlq(A,b,tol) specifies the tolerance of the method. If tol is [], then symmlq uses the default, 1e-6.

symmlq(A,b,tol,maxit) specifies the maximum number of iterations. If maxit is [], then symmlq uses the default, min(n,20).

symmlq(A,b,tol,maxit,M) and symmlq(A,b,tol,maxit,M1,M2) use the symmetric positive definite preconditioner M or M = M1*M2 and effectively solve the system inv(sqrt(M))*A*inv(sqrt(M))*y = inv(sqrt(M))*b for y and then return x = in(sqrt(M))*y. If M is [] then symmlq applies no preconditioner. M can be a function handle mfun such that mfun(x) returns M\x.

symmlq(A,b,tol,maxit,M1,M2,x0) specifies the initial guess. If x0 is [], then symmlq uses the default, an all-zero vector.

[x,flag] = symmlq(A,b,...) also returns a convergence flag.

Flag	Convergence
0	symmlq converged to the desired tolerance tol within maxit iterations.
1	symmlq iterated maxit times but did not converge.
2	Preconditioner M was ill-conditioned.
3	symmlq stagnated. (Two consecutive iterates were the same.)
4	One of the scalar quantities calculated during symmlq became too small or too large to continue computing.
5	Preconditioner M was not symmetric positive definite.

Whenever flag is not 0, the solution x returned is that with minimal norm residual computed over all the iterations. No messages are displayed if the flag output is specified.

[x,flag,relres] = symmlq(A,b,...) also returns the relative residual norm(b-A*x)/norm(b). If flag is 0, relres <= tol.

[x,flag,relres,iter] = symmlq(A,b,...) also returns the iteration number at which x was computed, where 0 <= iter <= maxit.

[x,flag,relres,iter,resvec] = symmlq(A,b,...) also returns a vector of estimates of the symmlq residual norms at each iteration, including norm(b-A*x0).

[x,flag,relres,iter,resvec,resveccg] = symmlq(A,b,...) also returns a vector of estimates of the conjugate gradients residual norms at each iteration.

Examples

Example 1

n = 100; 
on = ones(n,1); 
A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
b = sum(A,2); 
tol = 1e-10; 
maxit = 50; M1 = spdiags(4*on,0,n,n);

x = symmlq(A,b,tol,maxit,M1);
symmlq converged at iteration 49 to a solution with relative 
residual 4.3e-015

Example 2

This example replaces the matrix A in Example 1 with a handle to a matrix-vector product function afun. The example is contained in the function run_symmlq that:

• Calls symmlq with the function handle @afun as its first argument.

• Contains afun as a nested function, so that all variables in run_symmlq are available to afun.

The following shows the code for run_symmlq:

function x1 = run_symmlq
n = 100; 
on = ones(n,1); 
A = spdiags([-2*on 4*on -on],-1:1,n,n);
b = sum(A,2); 
tol = 1e-8; 
maxit = 15;
M1 = spdiags([on/(-2) on],-1:0,n,n); 
M2 = spdiags([4*on -on],0:1,n,n);
x1 = symmlq(@afun,b,tol,maxit,M1);

       function y = afun(x)
          y = 4 * x;
          y(2:n) = y(2:n) - 2 * x(1:n-1);
          y(1:n-1) = y(1:n-1) - 2 * x(2:n);
       end
end

When you enter

x1=run_symmlq;

MATLAB^® software displays the message

symmlq converged at iteration 49 to a solution with relative 
residual 4.3e-015

Example 3

Use a symmetric indefinite matrix that fails with pcg.

A = diag([20:-1:1,-1:-1:-20]);
b = sum(A,2);      % The true solution is the vector of all ones.
x = pcg(A,b);      % Errors out at the first iteration.
pcg stopped at iteration 1 without converging to the desired
tolerance 1e-006 because a scalar quantity became too small or 
too large to continue computing. 
The iterate returned (number 0) has relative residual 1

However, symmlq can handle the indefinite matrix A.

x = symmlq(A,b,1e-6,40);
symmlq converged at iteration 39 to a solution with relative 
residual 1.3e-007

References

[1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

[2] Paige, C. C. and M. A. Saunders, "Solution of Sparse Indefinite Systems of Linear Equations." SIAM J. Numer. Anal., Vol.12, 1975, pp. 617-629.

Extended Capabilities

Distributed Arrays

Partition large arrays across the combined memory of your cluster using Parallel Computing Toolbox™.

Usage notes and limitations:

• If M1 is a function, then it is applied independently to each row.

For more information, see Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox).