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.2 Linear Algebra on Sparse Matrices
Octave includes a polymorphic solver for sparse matrices, where the exact solver used to factorize the matrix depends on the properties of the sparse matrix itself. Generally, the cost of determining the matrix type is small relative to the cost of factorizing the matrix, and the type is cached once it is calculated, so that it is not redetermined each time the matrix is used.
The selection tree for how the linear equation is solved is:
 If the matrix is diagonal, solve directly and goto 8
 If the matrix is a permuted diagonal, solve directly taking into account the permutations. Goto 8
 If the matrix is square, banded and if the band density is less
than that given by
spparms ("bandden")
continue, else goto 4. If the matrix is tridiagonal and the righthand side is not sparse
continue, else goto 3b.
 If the matrix is hermitian, with a positive real diagonal, attempt Cholesky factorization using lapack xPTSV.
 If the above failed or the matrix is not hermitian with a positive real diagonal use Gaussian elimination with pivoting using lapack xGTSV, and goto 8.
 If the matrix is hermitian with a positive real diagonal, attempt Cholesky factorization using lapack xPBTRF.
 if the above failed or the matrix is not hermitian with a positive real diagonal use Gaussian elimination with pivoting using lapack xGBTRF, and goto 8.
 If the matrix is tridiagonal and the righthand side is not sparse
continue, else goto 3b.
 If the matrix is upper or lower triangular perform a sparse forward or backward substitution, and goto 8
 If the matrix is a upper triangular matrix with column permutations or lower triangular matrix with row permutations, perform a sparse forward or backward substitution, and goto 8
 If the matrix is square, hermitian with a real positive diagonal, attempt sparse Cholesky factorization using cholmod.
 If the sparse Cholesky factorization failed or the matrix is not hermitian with a real positive diagonal, and the matrix is square, factorize using umfpack.
 If the matrix is not square, or any of the previous solvers flags a singular or near singular matrix, find a minimum norm solution using cxsparse^{(19)}.
The band density is defined as the number of nonzero values in the matrix
divided by the number of nonzero values in the matrix. The banded matrix
solvers can be entirely disabled by using spparms
to set bandden
to 1 (i.e. spparms ("bandden", 1)
).
The QR solver factorizes the problem with a DulmageMendhelsohn decomposition, to
separate the problem into blocks that can be treated as overdetermined,
multiple well determined blocks, and a final overdetermined block. For
matrices with blocks of strongly connected nodes this is a big win as
LU decomposition can be used for many blocks. It also significantly
improves the chance of finding a solution to overdetermined problems
rather than just returning a vector of NaN
's.
All of the solvers above can calculate an estimate of the condition number. This can be used to detect numerical stability problems in the solution and force a minimum norm solution to be used. However, for narrow banded, triangular or diagonal matrices, the cost of calculating the condition number is significant, and can in fact exceed the cost of factoring the matrix. Therefore the condition number is not calculated in these cases, and Octave relies on simpler techniques to detect singular matrices or the underlying lapack code in the case of banded matrices.
The user can force the type of the matrix with the matrix_type
function. This overcomes the cost of discovering the type of the matrix.
However, it should be noted that identifying the type of the matrix incorrectly
will lead to unpredictable results, and so matrix_type
should be
used with care.
 Function File: [n, c] = normest (a, tol)
 Estimate the 2norm of the matrix a using a power series
analysis. This is typically used for large matrices, where the cost
of calculating the
norm (a)
is prohibitive and an approximation to the 2norm is acceptable.tol is the tolerance to which the 2norm is calculated. By default tol is 1e6. c returns the number of iterations needed for
normest
to converge.
 Function File: [est, v] = condest (a, t)
 Function File: [est, v] = condest (a, solve, solve_t, t)
 Function File: [est, v] = condest (apply, apply_t, solve, solve_t, n, t)

Estimate the 1norm condition number of a matrix A using t test vectors using a randomized 1norm estimator. If t exceeds 5, then only 5 test vectors are used.
If the matrix is not explicit, e.g. when estimating the condition number of a given an LU factorization,
condest
uses the following functions: apply

A*x
for a matrixx
of size n by t.  apply_t

A'*x
for a matrixx
of size n by t.  solve

A \ b
for a matrixb
of size n by t.  solve_t

A' \ b
for a matrixb
of size n by t.
The implicit version requires an explicit dimension n.
condest
uses a randomized algorithm to approximate the 1norms.condest
returns the 1norm condition estimate est and a vector v satisfyingnorm (A*v, 1) == norm (A, 1) * norm (v, 1) / est
. When est is large, v is an approximate null vector.References:
 Nicholas J. Higham and Françoise Tisseur, “A Block Algorithm for Matrix 1Norm Estimation, with an Application to 1Norm Pseudospectra.” SIMAX vol 21, no 4, pp 11851201.^{(20)}
See also norm, cond, onenormest
 Loadable Function: spparms ()
 Loadable Function: vals = spparms ()
 Loadable Function: [keys, vals] = spparms ()
 Loadable Function: val = spparms (key)
 Loadable Function: spparms (vals)
 Loadable Function: spparms ('defaults')
 Loadable Function: spparms ('tight')
 Loadable Function: spparms (key, val)
 Sets or displays the parameters used by the sparse solvers and factorization
functions. The first four calls above get information about the current
settings, while the others change the current settings. The parameters are
stored as pairs of keys and values, where the values are all floats and the
keys are one of the following strings:
spumoni
 Printing level of debugging information of the solvers (default 0)
ths_rel
 Included for compatibility. Not used. (default 1)
ths_abs
 Included for compatibility. Not used. (default 1)
exact_d
 Included for compatibility. Not used. (default 0)
supernd
 Included for compatibility. Not used. (default 3)
rreduce
 Included for compatibility. Not used. (default 3)
wh_frac
 Included for compatibility. Not used. (default 0.5)
autommd

Flag whether the LU/QR and the
\
and/
operators will automatically use the sparsity preserving mmd functions (default 1) autoamd

Flag whether the LU and the
\
and/
operators will automatically use the sparsity preserving amd functions (default 1) piv_tol
 The pivot tolerance of the umfpack solvers (default 0.1)
bandden
 Band density value for polymorphic solver (default 0.5)
umfpack

Flag whether the umfpack or mmd solvers are used for the LU,
\
and/
operations (default 1)
The value of individual keys can be set with
spparms (key, val)
. The default values can be restored with the special keyword 'defaults'. The special keyword 'tight' can be used to set the mmd solvers to attempt for a sparser solution at the potential cost of longer running time.
 Loadable Function: p = sprank (s)

Calculates the structural rank of a sparse matrix s. Note that only the structure of the matrix is used in this calculation based on a DulmageMendelsohn permutation to block triangular form. As such the numerical rank of the matrix s is bounded by
sprank (s) >= rank (s)
. Ignoring floating point errorssprank (s) == rank (s)
.See also dmperm
 Loadable Function: [count, h, parent, post, r] = symbfact (s, typ, mode)

Performs a symbolic factorization analysis on the sparse matrix s. Where
 s
 s is a complex or real sparse matrix.
 typ

Is the type of the factorization and can be one of
sym
 Factorize s. This is the default.
col

Factorize
s' * s
. row

Factorize
s * s'
. lo

Factorize
s'
 mode
 The default is to return the Cholesky factorization for r, and if mode is 'L', the conjugate transpose of the Cholesky factorization is returned. The conjugate transpose version is faster and uses less memory, but returns the same values for count, h, parent and post outputs.
The output variables are
 count
 The row counts of the Cholesky factorization as determined by typ.
 h
 The height of the elimination tree.
 parent
 The elimination tree itself.
 post
 A sparse boolean matrix whose structure is that of the Cholesky factorization as determined by typ.
ISBN 095461206X  GNU Octave Manual Version 3  See the print edition 