-----------------------------------------------------------------
1. Introduction
-----------------------------------------------------------------
HSL_MA97 solves the system
A * X = B
where A is a specified n x n symmetric or Hermitian matrix, B is an n x m
specified right-hand side and X is an n x m unknown. If m=1, this simplifies
to A * x = b.
A direct method is used, performing the factorization PAP' = LDL' where P is
a permutation matrix.
If the A is positive definite, the factorization PAP' = LL' may be used
by setting control.pos_def=true.
OpenMP is used to support parallel computation on multicore. For best
performance an OpenMP supporting compiler should be used.
-----------------------------------------------------------------
2. Requirements
-----------------------------------------------------------------
These instructions are for linux-based systems (both 32-bit and 64-bit machines
are supported).
Requirements:
- Matlab 2008a or more recent version
- The gfortran or g95 compiler
- The MATLAB environment variable must point to your system matlab directory
(see INSTALL for further details)
We note that gfortran compiler is NOT SUPPORTED by Mathworks on Linux prior
to MATLAB 2011a, however g95 does not support OpenMP and has considerably
poorer code optimization.
This code requires OpenMP 3.0, which is not supported under gfortran 4.3.
Therefore you must use the serial version under Linux until MATLAB supports a
more modern Fortran compiler (gfortran 4.4 or newer). Please direct complaints
on this matter to Mathworks.
-----------------------------------------------------------------
3. Installation
-----------------------------------------------------------------
Instructions for installing the Matlab Interface for HSL_MA97 are located in
INSTALL.
-----------------------------------------------------------------
4. Using the Matlab interface
-----------------------------------------------------------------
There are two approaches for using HSL_MA97 under Matlab.
(a) Just use hsl_ma97_backslash in place of matlab's own X = A \ B.
(b) Use the "expert" interface that allows a factorization to be preserved and
multiple solves to be performed using it.
Approach (a) is offered as a subset of the functions implement (b).
-----------------------------------------------------------------
4(a). As a replacement backslash routine
-----------------------------------------------------------------
- If not already in the search paths, add the directory where you installed the
interface to the search paths, e.g.
>> addpath('ma97_matlab')
>> javaaddpath('ma97_matlab')
where ma97_matlab is the directory.
[ You can add these paths permanently (see 'help pathtool')]
- To solve the equation AX = B for X:
>> X = ma97_backslash(A, B)
This will only work when A is symmetric. For Hermitian A, or to perform
the faster LL' factorization of a positive-definite matrix, see the expert
interface.
-----------------------------------------------------------------
4(b). The "expert" interface
-----------------------------------------------------------------
The expert interface has the concept of handles. These are integers that refer
to factorizations held in memory. The factorization will continue to take up
memory until the 'destroy' call is used on that particular handle. The use
of these handles allows for storing multiple matrix factors simultaneously
without the inefficiency of translating internal data formats to MATLAB arrays.
Complex matrices are assumed to be symmetric unless they are specified as
Hermitian by setting the value control.hermitian=true.
Real symmetric and complex Hermitian matrices are assumed to be indefinite
unless they are specified as being positive definite by setting the value
control.pos_def=true.
- If not already in the search paths, add the directory where you installed the
interface to the search paths, e.g.
>> addpath('ma97_matlab')
>> javaaddpath('ma97_matlab')
where ma97_matlab is the directory.
[ You can add these paths permanently (see 'help pathtool') ]
- hsl_ma97_factor may be used to perform the factorization of a symmetric
matrix A. This is equivalent to calling the Fortran routines ma97_analyse
and ma97_factor.
>> handle = hsl_ma97_factor(A)
Optionally control and info structures may be used
>> [handle, info] = hsl_ma97_factor(A, control)
If the user wishes to provide a permutation P (rather than allow ma97 to find
its own), it should be specified as the fourth argument:
>> handle = hsl_ma97_factor(A, [], P)
where P is a vector as returned by symamd(A).
- hsl_ma97_solve may be used to perform the solution of AX=B where A has
already been factorized by a previous call to hsl_ma97_factor or
hsl_ma97_backslash. This is equivalent to calling the Fortran routine
ma97_solve or ma97_solve_fredholm.
>> X = hsl_ma97_solve(handle, B)
Optionally control and info structures may be used
>> [X, info] = hsl_ma97_solve(handle, B, control)
If the matrix is singular, a Fredholm alternative may be calculated by
>> [X, info, consist] = hsl_ma97_solve(handle, B, control)
where control and info are again optional. If B is n x k, then alternative
will be a 1 x k Logical matrix. If right-hand side i is consistent,
consist(i) will be true and X(:,i) will provide a solution. If the right-hand
side i is inconsistent, consist(i) will be false, X(:,i) will provide a
solution which is consistent for a non-singular submatrix of $A$, and
X(:,i+size(B,2)) will provide a vector that satisfies the Fredholm
alternative: A*x=0 and x'*b!=0.
- hsl_ma97_backslash may be used to perform a combined factorization and
solve. This is equivalent to calling the Fortran routines ma97_analyse
and ma97_factor_solve.
>> X = hsl_ma97_backslash(A, b)
Optionally control and info structures may be used
>> [X, info] = hsl_ma97_backslash(A, b, control)
It is possible to preserve the factorization in memory by supplying an
output variable to store the handle
>> [X, info, handle] = hsl_ma97_backslash(A, b)
As in the 'factor' call a permutation P may be supplied by the user
>> X = hsl_ma97_backslash(A, b, [], P)
- hsl_ma97_multiply may be used to form any of the products LX, L^TX, L^{-1}X,
or L^{-T}X where L is understood to include scaling and permutations.
>> Y = hsl_ma97_multiply(handle, transpose, inverse, X)
will form the appropriate product. If transpose is true, then the transpose is
used, and if inverse is used, the inverse is used.
In the case L^{-1}X where X is a sparse n x 1 matrix, an efficient method is
used, but in all other cases sparse vectors are converted to be dense before
proceeding.
Depending on the arguments, this method is equivalent to one of the Fortran
routines ma97_solve, ma97_sparse_fwd_solve or ma97_lmultiply.
- hsl_ma97_destroy may be used to free the memory and resources associated with
a previous factorization.
>> hsl_ma97_destroy(handle)
frees all memory associated with handle: the factorization may not be reused
thereafter.
-----------------------------------------------------------------
5. Control and information structures
-----------------------------------------------------------------
A limited subset of the Fortran control and information parameters may be
used through the expert MATLAB interface.
-----------------------------------------------------------------
5(a). The control structure
-----------------------------------------------------------------
The argument control may have any of the following components set. If
unrecognised components are present a warning is issued. If a component
is not present its default value is used.
control.hermitian - True or false. Determines if a complex matrix is
treated as Hermitian (true) or symmetric (false).
Default is false.
control.nemin - Maximum number of columns in candidates for
supernode amalgamation. Default is 8.
control.num_threads - Number of threads on which to run. Default is the
maximum available.
control.ordering - If an explicit P is not supplied, determines the
ordering algorithm employed:
1 : Approximate minimum Degree
2 : Minimum Degree
3 : METIS (will need to be compiled against
the METIS library)
4 : MA47 ordering
5 : Heuristic choice between AMD and METIS for
parallel computation.
6 : Heuristic choice between AMD and METIS for
serial computation.
7 : Matching based ordering using AMD.
8 : Matching based ordering using METIS.
Default is 5.
Note that for 7 and 8, MC64 scaling is calculated as
a side-effect and can therefore be used for free.
control.pos_def - True or false. Determines if a matrix is treated
as positive-definite (true) or indefinite (false).
In the complex case only Hermitian matrices can
be positive definite.
Default is false.
control.scaling - Determines if scaling is to be used with values:
<=0 : No scaling
1 : MC64 Bipartite matching based scaling
2 : MC77 (1 itr in inf norm, 3 itr in one norm)
>=4 : MC30
Default is 0.
control.solve_blas3 - True if we should use BLAS3 for solve, false for BLAS2.
In multiple rhs case BLAS3 is always used.
Default is false.
control.solve_mf - True or false. Determines whether to use a parallel
multifrontal-style solve or not. If false, a serial
supernodal-style solve is used. Typically the former is
faster on large problems, and the latter on small ones.
The default is false.
control.small - Pivots of modulus less than this are treated as zero.
Default is 1e-20.
control.u - Initial relative pivot tolerance threshold. Default
is 0.01.
A fuller description for many of the components is available in the Fortran
documentation.
-----------------------------------------------------------------
5(b). The info structure
-----------------------------------------------------------------
On return from a routine, the output argument info is a structure that will
have one or more of the following components set:
info.matrix_rank - Number of non-zero pivots.
info.num_delay - Number of delayed pivots.
info.num_factor - Number of entries in the factors (after supernode
amalgamation).
info.num_flops - Number of floating point operations to form factors
(after supernode amalgamation).
info.num_neg - Number of negative pivots in factors.
info.num_two - Number of 2x2 pivots used in factorization.
info.order - Ordering used. One of 'AMD', 'MD', 'MC47', 'MeTiS' or
'user'.
info.analyse_time - Wall clock time for Fortran ma97_analyse call
info.factor_time - Wall clock time for Fortran ma97_factor call
info.factor_solve_time - Wall clock time for Fortran ma97_factor_solve call
info.solve_time - Wall clock time for Fortran ma97_solve call
A fuller description for many of the components is available in the Fortran
documentation.
-----------------------------------------------------------------
6. Example
-----------------------------------------------------------------
The following MATLAB session shows an example of using hsl_ma97 to solve a
system using the expert 'backslash' action.
>> A = sparse ([1 1 1 2 2 3 3 3 4 4], [2 3 4 1 3 1 2 3 1 4], [1.1-i 2.2-i 3.3-i, 1.1+i 4.4-i, 2.2+i 4.4+i 5.5, 3.3+i 6.6])
A =
(2,1) 1.1000 + 1.0000i
(3,1) 2.2000 + 1.0000i
(4,1) 3.3000 + 1.0000i
(1,2) 1.1000 - 1.0000i
(3,2) 4.4000 + 1.0000i
(1,3) 2.2000 - 1.0000i
(2,3) 4.4000 - 1.0000i
(3,3) 5.5000
(1,4) 3.3000 - 1.0000i
(4,4) 6.6000
>> x = rand(size(A,1), 2)
x =
0.8092 0.3258
0.7486 0.5464
0.1202 0.3989
0.5250 0.4151
>> b = A*x;
>> control.hermitian = true;
>> [soln, info] = hsl_ma97_backslash(A, b, control)
soln =
0.8092 0.3258
0.7486 - 0.0000i 0.5464 + 0.0000i
0.1202 - 0.0000i 0.3989
0.5250 0.4151 + 0.0000i
info =
matrix_rank: 4
num_delay: 0
num_factor: 10
num_flops: 30
num_neg: 2
num_two: 1
order: 'AMD'
analyse_time: 1.0000e-4
factor_solve_time: 1.0000e-4