-----------------------------------------------------------------
Contents
-----------------------------------------------------------------
(1) Introduction
(2) Requirements
(3) Installation
(4) Using the Matlab interface
(5) control, info and rinfo arguments
-----------------------------------------------------------------
1. Introduction
-----------------------------------------------------------------
Given a sparse matrix A = a_{ij} of dimension n, ME57 solvesthe linear system
A x = b
ME57 offers two MATLAB procedures:
me57_factor
me57_solve
me57_factor performs an LDL' factorization of a permutation of a sparse A,
using the HSL2011 package ME57. It is equivalent to calling the Fortran
routines MAE7AD, ME57BD, and (if required) ME57ED.
me57_solve uses the factorization in order to compute the solution(s)
corresponding to the right-hand side(s) b. It is equivalent to calling
the Fortran routines ME57CD and ME57DD.
ME57 routines implement a multifrontal algorithm for the
factorization of A. The matrix L is a unit lower triangular matrix and
D is block diagonal with blocks of order 1 or 2. ME57 uses only
the diagonal and the lower triangle of A.
The upper triangle is assumed to be the transpose of the lower.
Before the factorization starts, scaling is performed using the HSL2011
package MC64.
In order to minimize the fill-in during the factorization process,
the user can choose to reorder the matrix A either by Approximate
Minimum Degree (AMD) (the default value) or, if the user wants to use
a permutation p, he/she can call the interfaces using A(p,p) and
setting the control.order parameter to a value > 1 ( see Section 5
below).
For stability, all the pivots are tested numerically. However, in
order to avoid excessive disruption of the AMD or the user's reorderings
threshold pivoting is used. Therefore, denoting the entries to be
processed at step k of the factorization by f_{ij}, a 1x1 numerical
pivot will be used if
|f_{kk} | \ge u \max_{j \ne k} |f_{kj} |,
or a 2x2 numerical pivot will be used if
(|f_{k,k}| |f_{k,k+1}| ) (\max_{j\ne k,k+1} |f_{k,j}| ) (1)
( ) ( ) <= u^-1( )
(|f_{k+1,k}| |f_{k+1,k+1}|) (\max_{j\ne k,k+1} |f_{k+1,j}|) (1)
where u is a threshold parameter given in the control.thres
parameter (the default value is 0.01), see Section 5 below.
During stage k of a multifrontal implementation
of Gaussian elimination, some entries may not be eligible to be chosen
as a pivot. In this case, the algorithm will proceed using the
eligible variables only and will choose the rejected pivot at a later
stage when the candidate entry will become eligible.
Additional details about the theory of the multifrontal method and the
numerical pivoting strategy can be found in [1] and [2].
References
[1] I.S. Duff "MA57: a new code for the solution of sparse symmetric
indefinite systems", Technical Report RAL-TR-2002-024, ACM
Trans. Math. Software 30 (2004), 118-144.
[2] I.S. Duff and S. Pralet, "Strategies for scaling and pivoting for
sparse symmetric indefinite problems", SIAM. J. Matrix Anal. & Appl. 27,
(2005), 313-340.
-----------------------------------------------------------------
2. Requirements
-----------------------------------------------------------------
These instructions should work on Linux, Mac OS or Windows systems.
They have been tested using the six most recent releases of Matlab.
Please note that HSL software requires a Fortran compiler to be installed
that is compatible with Matlab. Please see
https://uk.mathworks.com/support/requirements/supported-compilers.html
and look up the requirements for your OS and version of Matlab.
The MATLAB environment variable must point to your system
matlab directory (see INSTALL for further details)
-----------------------------------------------------------------
3. Installation
-----------------------------------------------------------------
Instructions for installing the Matlab interface to ME57 are
located in INSTALL. You can type me57_install without any
input arguments at the command line of Matlab to proceed to the
standard installation.
-----------------------------------------------------------------
4. Using the Matlab interface
-----------------------------------------------------------------
To solve Ax=b, first a factorization of A must be performed.
>> struct = me57_factor(A)
This is followed by a solve using the factors stored in struct, which must
not be modified by the user.
>> x = me57_solve(A, b, struct)
If multiple systems are to be solved using the same matrix A, it need only
be factorized once. For example,
>> struct = me57_factor(A)
>> x1 = me57_solve(A, b1, struct)
>> x2 = me57_solve(A, b2, struct)
Additional parameters may be used to control the behaviour of both functions
and statistics may also be collected:
>> [struct, info, rinfo] = me57_factor(A, control)
>> [x, info, rinfo] = me57_solve(A, b, struct, nitre)
The arguments are described below.
-----------------------------------------------------------------
4(a). The routine me57_factor
-----------------------------------------------------------------
struct = me57_factor(kind, A)
[struct, info, rinfo] = me57_factor(kind, A, control)
If A is not square or sparse an ERROR message is printed. If a
fatal error occurs during the factorization an error message is
printed. Fatal errors include the wrong number of input or output arguments
and insufficient memory.
If the numerical rank < length(A) a WARNING message is printed.
INPUT
kind integer parameter that the user must set to 1 if A complex
symmetric. For any other value, A will be interpreted as
Hermitian.
A sparse symmetric or Hermitian input matrix (only the lower
triangular part will be used). Note that if A is not symmetric or
Hermitian, me57_factor and me57_solve will expand the the lower
triangular entries of A to form a symmetric matrix.
control OPTIONAL structure whose components affect the behaviour of
ME57 as detailed in Section 5.
OUTPUT
struct structure containing the factorization produced by ME57. This
is not intended to be accessed by the user, and should be
passed unaltered to subsequent calls to me57_solve.
info OPTIONAL array containing statistics. The values will be set as
described in Section 5.
rinfo OPTIONAL array containing statistics. The values will be set as
described in Section 5.
-----------------------------------------------------------------
4(b). The routine me57_solve
-----------------------------------------------------------------
x = me57_solve(kind, A, b, struct)
[x, info, rinfo] = me57_solve(kind, A, b, struct, nitre)
INPUT
kind integer parameter that the user must set to 1 if A complex
symmetric. For any other value, A will be interpreted as
Hermitian. This must be the same as for the call to
me57_factor.
A sparse symmetric input matrix (only the lower triangular part
will be used). This must be the same as for the call to
me57_factor.
b an n x m matrix containing m right-hand sides.
struct the factors returned by a previous call to me57_factor.
nitre OPTIONAL number specifying maximum number of steps for
iterative refinement. Iterative refinement can be disabled
by setting this value to 0. If it is not present then at
most 10 steps of iterative refinement are performed.
OUTPUT
x an n x m matrix containing the solution to Ax=b.
info OPTIONAL array containing statistics. The values will be set as
described in Section 5.
rinfo OPTIONAL array containing statistics. The values will be set as
described in Section 5.
-----------------------------------------------------------------
5. control, info and rinfo arguments
-----------------------------------------------------------------
control
-------
If any of the fields of control is omitted the default value is
automatically set. Incorrect values or empty values produce an error
message.
control.order
Allows the user to choose a reordering strategy for the pivot
sequence:
control.order <= 0 ==> AMD ordering is used
control.order > 1 ==> User's reordering p is used and
the user must pass the
permuted matrix A(p,p) in the
input
Its default value is 0.
control.thres
The value of the threshold u (see ME57 Specification sheet).
Its default value is 0.01.
control.stpv
If stpv is set to a value greater than zero, static pivoting
is invoked. This means that if, at any stage of the
multifrontal elimination, the threshold criterion defined by
control.thres prevents us choosing a pivot, then we first weaken
the threshold by successively trying values one tenth of the
previous until a threshold of sqrt(thresh)*stvp is reached.
If the pivot is still not chosen, but involves non eligible
variables, then the diagonal entries are replaced by
control.stpv times the largest modulus of an entry in the
original matrix and are used as 1 by 1 pivots in the
factorization. The default value is 0.
We advise the user interested in this option to use first a
value close to sqrt(eps) for stpv. Iterative refinement must
be used but convergence is not guaranteed.
info
----
info(1) will be 0 after a successful execution;
a positive value indicates that the matrix A has
numerical pseudo rank < length(A). The solution can be
incorrect, and info(25) holds the numerical pseudo rank
value. If A is not Hermitian an error message is printed.
info(2) to info(4) are only used internally.
info(5) Forecast number of real in the factors.
info(6) Forecast number of integers in the factors.
info(7) Forecast maximum frontal matrix size.
info(8) Forecast number of nodes in the assembly tree.
info(9) Forecast of the length of fact array(real)
without numerical pivoting.
info(10) Forecast of the length of ifact array(integer)
without numerical pivoting.
info(11) Length of fact required for a successful
completion of the numerical phase allowing
data compression (without numerical pivoting).
info(12) Length of ifact required for a successful
completion of the numerical phase allowing
data compression (without numerical pivoting).
info(13) Number of data compresses.
info(14) Number of entries in factors.
info(15) Storage for real data in factors.
info(16) Storage for integer data in factors.
info(17) Minimum length of fact required for a successful
completion of the numerical phase.
info(18) Minimum length of ifact required for a successful
completion of the numerical phase.
info(19) Minimum length of fact required for a successful
completion of the numerical phase without
data compression.
info(20) Minimum length of ifact required for a successful
completion of the numerical phase without data compression.
info(21) Order of the largest frontal matrix.
info(22) Number of 2x2 pivots.
info(23) Total number of fully-summed variables passed to
the father node because of numerical pivoting.
info(24) Number of negative eigenvalues.
info(25) Rank of the factorization of the matrix.
info(26) Only used internally
info(27) Pivot step where static pivoting commences and
it is set to zero if no modification performed.
info(28) Number of data compresses on real data structures.
info(29) Number of data compresses on integer data structures.
info(30) Number of steps performed by Iterative Refinement.
info(31) Number of block pivots in factors (see info(8)).
info(32) Number of zeros in the triangle of the factors.
info(33) Number of zeros in the rectangle of the factors.
info(34) Number of zero columns in rectangle of the factors.
info(35) Number of pivots chosen in the static pivoting phase.
info(36) to info(40) are only used internally.
rinfo
-----
rinfo(1) Forecast number of floating_point additions
for the assembly processes.
rinfo(2) Forecast number of floating_point operations
to perform the elimination operations
counting a multiply-add pair as two operations
rinfo(3) Number of floating_point additions
for the assembly processes.
rinfo(4) Number of floating_point operations
to perform the elimination operations
counting multiply-add pair as two operations.
rinfo(5) Only used internally
rinfo(6) Omega1 (backward error) (see ME57 spec).
rinfo(7) Omega2 (backward error) (see ME57 spec).
rinfo(8) to rinfo(10) are only used internally
rinfo(11) Cond number for Omega1
rinfo(12) Cond number for Omega2
rinfo(13) Estimated forward error
rinfo(14) to rinfo(15) are only used internally.
rinfo(16) Minimum value of the scaling factor (MC64).
rinfo(17) Maximum value of the scaling factor (MC64).
rinfo(18) Maximum modulus of matrix entry.
rinfo(19) and rinfo(20) are only used internally.
rinfo(21) Analysis phase cputime (sec).
rinfo(22) Factorization phase cputime (sec).
rinfo(23) Solution phase cputime (sec) without iterative refinement.
rinfo(24) Solution phase cputime (sec) with iterative refinement.