GNU Octave Manual Version 3 by John W. Eaton, David Bateman, Søren Hauberg Paperback (6"x9"), 568 pages ISBN 095461206X RRP £24.95 ($39.95) 
20.3 Iterative Techniques applied to sparse matrices
The left division \
and right division /
operators,
discussed in the previous section, use direct solvers to resolve a
linear equation of the form x = A \ b
or
x = b / A
. Octave equally includes a number of
functions to solve sparse linear equations using iterative techniques.
 Function File: x = pcg (a, b, tol, maxit, m1, m2, x0, ...)
 Function File: [x, flag, relres, iter, resvec, eigest] = pcg (...)

Solves the linear system of equations
a * x = b
by means of the Preconditioned Conjugate Gradient iterative method. The input arguments are
a can be either a square (preferably sparse) matrix or a
function handle, inline function or string containing the name
of a function which computes
a * x
. In principle a should be symmetric and positive definite; ifpcg
finds a to not be positive definite, you will get a warning message and the flag output parameter will be set.  b is the right hand side vector.

tol is the required relative tolerance for the residual error,
b  a * x
. The iteration stops ifnorm (b  a * x) <= tol * norm (b  a * x0)
. If tol is empty or is omitted, the function setstol = 1e6
by default. 
maxit is the maximum allowable number of iterations; if
[]
is supplied formaxit
, orpcg
has less arguments, a default value equal to 20 is used. 
m = m1 * m2 is the (left) preconditioning matrix, so that the iteration is
(theoretically) equivalent to solving by
pcg
P * x = m \ b
, withP = m \ a
. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method. Instead of matrices m1 and m2, the user may pass two functions which return the results of applying the inverse of m1 and m2 to a vector (usually this is the preferred way of using the preconditioner). If[]
is supplied for m1, or m1 is omitted, no preconditioning is applied. If m2 is omitted, m = m1 will be used as preconditioner.  x0 is the initial guess. If x0 is empty or omitted, the function sets x0 to a zero vector by default.
The arguments which follow x0 are treated as parameters, and passed in a proper way to any of the functions (a or m) which are passed to
pcg
. See the examples below for further details. The output arguments are
x is the computed approximation to the solution of
a * x = b
. 
flag reports on the convergence.
flag = 0
means the solution converged and the tolerance criterion given by tol is satisfied.flag = 1
means that the maxit limit for the iteration count was reached.flag = 3
reports that the (preconditioned) matrix was found not positive definite.  relres is the ratio of the final residual to its initial value, measured in the Euclidean norm.
 iter is the actual number of iterations performed.

resvec describes the convergence history of the method.
resvec (i,1)
is the Euclidean norm of the residual, andresvec (i,2)
is the preconditioned residual norm, after the (i1)th iteration,i = 1, 2, ..., iter+1
. The preconditioned residual norm is defined asnorm (r) ^ 2 = r' * (m \ r)
wherer = b  a * x
, see also the description of m. If eigest is not required, onlyresvec (:,1)
is returned. 
eigest returns the estimate for the smallest
eigest (1)
and largesteigest (2)
eigenvalues of the preconditioned matrixP = m \ a
. In particular, if no preconditioning is used, the estimates for the extreme eigenvalues of a are returned.eigest (1)
is an overestimate andeigest (2)
is an underestimate, so thateigest (2) / eigest (1)
is a lower bound forcond (P, 2)
, which nevertheless in the limit should theoretically be equal to the actual value of the condition number. The method which computes eigest works only for symmetric positive definite a and m, and the user is responsible for verifying this assumption.
Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A)
N = 10; A = spdiag ([1:N]); b = rand (N, 1); [L, U, P, Q] = luinc (A,1.e3);
Example 1: Simplest use of
pcg
x = pcg(A,b)
Example 2:
pcg
with a function which computesa * x
function y = applyA (x) y = [1:N]'.*x; endfunction x = pcg ("applyA", b)
Example 3:
pcg
with a preconditioner: l * ux=pcg(A,b,1.e6,500,L*U);
Example 4:
pcg
with a preconditioner: l * u. Faster than Example 3 since lower and upper triangular matrices are easier to invertx=pcg(A,b,1.e6,500,L,U);
Example 5: Preconditioned iteration, with full diagnostics. The preconditioner (quite strange, because even the original matrix a is trivial) is defined as a function
function y = applyM(x) K = floor (length (x)  2); y = x; y(1:K) = x(1:K)./[1:K]'; endfunction [x, flag, relres, iter, resvec, eigest] = ... pcg (A, b, [], [], "applyM"); semilogy (1:iter+1, resvec);
Example 6: Finally, a preconditioner which depends on a parameter k.
function y = applyM (x, varargin) K = varargin{1}; y = x; y(1:K) = x(1:K)./[1:K]'; endfunction [x, flag, relres, iter, resvec, eigest] = ... pcg (A, b, [], [], "applyM", [], [], 3)
References
[1] C.T.Kelley, “Iterative methods for linear and nonlinear equations”, SIAM, 1995 (the base PCG algorithm) [2] Y.Saad, “Iterative methods for sparse linear systems”, PWS 1996 (condition number estimate from PCG). A revised version of this book is available online at http://wwwusers.cs.umn.edu/~saad/books.html.
See also sparse, pcr

a can be either a square (preferably sparse) matrix or a
function handle, inline function or string containing the name
of a function which computes
 Function File: x = pcr (a, b, tol, maxit, m, x0, ...)
 Function File: [x, flag, relres, iter, resvec] = pcr (...)

Solves the linear system of equations
a * x = b
by means of the Preconditioned Conjugate Residuals iterative method. The input arguments are
a can be either a square (preferably sparse) matrix or a
function handle, inline function or string containing the name
of a function which computes
a * x
. In principle a should be symmetric and nonsingular; ifpcr
finds a to be numerically singular, you will get a warning message and the flag output parameter will be set.  b is the right hand side vector.

tol is the required relative tolerance for the residual error,
b  a * x
. The iteration stops ifnorm (b  a * x) <= tol * norm (b  a * x0)
. If tol is empty or is omitted, the function setstol = 1e6
by default. 
maxit is the maximum allowable number of iterations; if
[]
is supplied formaxit
, orpcr
has less arguments, a default value equal to 20 is used. 
m is the (left) preconditioning matrix, so that the iteration is
(theoretically) equivalent to solving by
pcr
P * x = m \ b
, withP = m \ a
. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method. Instead of matrix m, the user may pass a function which returns the results of applying the inverse of m to a vector (usually this is the preferred way of using the preconditioner). If[]
is supplied for m, or m is omitted, no preconditioning is applied.  x0 is the initial guess. If x0 is empty or omitted, the function sets x0 to a zero vector by default.
The arguments which follow x0 are treated as parameters, and passed in a proper way to any of the functions (a or m) which are passed to
pcr
. See the examples below for further details. The output arguments are
x is the computed approximation to the solution of
a * x = b
. 
flag reports on the convergence.
flag = 0
means the solution converged and the tolerance criterion given by tol is satisfied.flag = 1
means that the maxit limit for the iteration count was reached.flag = 3
reports tpcr
breakdown, see [1] for details.  relres is the ratio of the final residual to its initial value, measured in the Euclidean norm.
 iter is the actual number of iterations performed.

resvec describes the convergence history of the method,
so that
resvec (i)
contains the Euclidean norms of the residual after the (i1)th iteration,i = 1,2, ..., iter+1
.
Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A)
N = 10; A = diag([1:N]); A = sparse(A); b = rand(N,1);
Example 1: Simplest use of
pcr
x = pcr(A, b)
Example 2:
pcr
with a function which computesa * x
.function y = applyA(x) y = [1:10]'.*x; endfunction x = pcr('applyA',b)
Example 3: Preconditioned iteration, with full diagnostics. The preconditioner (quite strange, because even the original matrix a is trivial) is defined as a function
function y = applyM(x) K = floor(length(x)2); y = x; y(1:K) = x(1:K)./[1:K]'; endfunction [x, flag, relres, iter, resvec] = ... pcr(A, b, [], [], 'applyM') semilogy([1:iter+1], resvec);
Example 4: Finally, a preconditioner which depends on a parameter k.
function y = applyM(x, varargin) K = varargin{1}; y = x; y(1:K) = x(1:K)./[1:K]'; endfunction [x, flag, relres, iter, resvec] = ... pcr(A, b, [], [], 'applyM', [], 3)
References
[1] W. Hackbusch, “Iterative Solution of Large Sparse Systems of Equations”, section 9.5.4; Springer, 1994
See also sparse, pcg

a can be either a square (preferably sparse) matrix or a
function handle, inline function or string containing the name
of a function which computes
The speed with which an iterative solver converges to a solution can be
accelerated with the use of a preconditioning matrix M. In this
case the linear equation M^1 * x = M^1 *
A \ b
is solved instead. Typical preconditioning matrices
are partial factorizations of the original matrix.
 Loadable Function: [l, u, p, q] = luinc (a, '0')
 Loadable Function: [l, u, p, q] = luinc (a, droptol)
 Loadable Function: [l, u, p, q] = luinc (a, opts)
 Produce the incomplete LU factorization of the sparse matrix a.
Two types of incomplete factorization are possible, and the type
is determined by the second argument to luinc.
Called with a second argument of '0', the zerolevel incomplete LU factorization is produced. This creates a factorization of a where the position of the nonzero arguments correspond to the same positions as in the matrix a.
Alternatively, the fillin of the incomplete LU factorization can be controlled through the variable droptol or the structure opts. The umfpack multifrontal factorization code by Tim A. Davis is used for the incomplete LU factorization.^{(21)}
droptol determines the values below which the values in the LU factorization are dropped and replaced by zero. It must be a positive scalar, and any values in the factorization whose absolute value are less than this value are dropped, expect if leaving them increase the sparsity of the matrix. Setting droptol to zero results in a complete LU factorization which is the default.
opts is a structure containing one or more of the fields
droptol

The drop tolerance as above. If opts only contains
droptol
then this is equivalent to using the variable droptol. milu

A logical variable flagging whether to use the modified incomplete LU
factorization. In the case that
milu
is true, the dropped values are subtracted from the diagonal of the matrix U of the factorization. The default isfalse
. udiag

A logical variable that flags whether zero elements on the diagonal of U
should be replaced with droptol to attempt to avoid singular
factors. The default is
false
. thresh
 Defines the pivot threshold in the interval [0,1]. Values outside that range are ignored.
All other fields in opts are ignored. The outputs from luinc are the same as for lu.
See also sparse, lu, cholinc
ISBN 095461206X  GNU Octave Manual Version 3  See the print edition 