-----------------------------------------------------------------
Contents
-----------------------------------------------------------------
(1) Introduction
(2) Requirements
(3) Installation
(4) Using the Matlab interface
(a) simple replacement for ichol
(b) general interface
(5) Control structure
(6) Information structure
-----------------------------------------------------------------
1. Introduction
-----------------------------------------------------------------
Let K be a sparse symmetric saddle-point matrix of the form
K = [ A B']
[ B -C ]
where A is (n-m) x (n-m) symmetric positive definite, B is rectangular
m x (n-m) and of full rank (m < n), and C is m x m symmetric positive
semi-definite. HSL_MI30 computes a signed incomplete Cholesky factorization.
Given an n x n sparse symmetric matrix A, HSL_MI30 computes a signed
incomplete Cholesky factorization LDL' that may be used as a preconditioner when
solving the system
K * x = b
using an iterative method.
The matrix K is optionally reordered, scaled and,
if necessary, shifted to avoid breakdown of the factorization so that the
incomplete factorization of
S Q' [ A B'] Q S + [alpha1*I 0 ]
[ B -C ] [ 0 alpha2*I ]
is computed, where Q is a permutation matrix, S is a diagonal scaling matrix
and alpha1 and alpha2 are non-negative shifts.
The incomplete factorization is based on a matrix decomposition of the form
(L+R) D (L+R)' + E,
where the incomplete factor L is a lower triangular matrix with
positive diagonal entries and D is a diagonal matrix with positive values corresponding
to the entries of the (1,1) block of K and negative values corresponding
to the entries of the (2,2) block of K. The matrix R is a strictly lower
triangular amtrix and E is an error matrix. The matrix R has entries that are
smaller in absolute value that those of L.
It is used to stabilize the factorization process but is then discarded.
The user can control the maximum number of entries within each column of
L and R and the dropping of small entries from L and R.
The matrix K is optionally reordered, scaled and, if necessary, shifted to
avoid breakdown of the factorization.
Further details are given in the Fortran user documentation and in the report:
J.A. Scott and M. Tuma. (2014).
On signed incomplete Cholesky factorization preconditioners for
saddle-point systems. Technical Report RAL-P-2014-003.
-----------------------------------------------------------------
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.
-----------------------------------------------------------------
3. Installation
-----------------------------------------------------------------
Instructions for installing the Matlab Interface for HSL_MI30 are located in
INSTALL.
-----------------------------------------------------------------
4. Using the Matlab interface
-----------------------------------------------------------------
Two interfaces are offered:
(a) hsl_mi30_sichol : this allows no preordering of K.
(b) general interface.
-----------------------------------------------------------------
4(a). Using the Matlab interface : sichol
-----------------------------------------------------------------
- If not already in the search paths, add the directory where you installed the
interface to the search paths, e.g.
>> addpath('mi30_matlab')
>> javaaddpath('mi30_matlab')
where mi30_matlab is the directory.
[ You can add these paths permanently (see 'help pathtool')]
- To compute a signed incomplete Cholesky factorization
>> [L, D, S] = hsl_mi30_sichol(K,m)
Options may be specified using a control structure (see Section 5)
>> [L, D, S] = hsl_mi30_sichol(K,m,control)
-----------------------------------------------------------------
4(a)(i). Using the Matlab interface: sichol example
-----------------------------------------------------------------
Below is an example of using the sichol interface
in conjunction with an iterative solver:
>> B = sparse(gallery('uniformdata',10,0));
>> K = [speye(10) B; B' sparse(10,10)];
>> b = ones(length(K),1);
>> control.lsize = 5;
>> control.rsize = 5;
>> [L, D, S] = hsl_mi30_sichol(K,10,control);
>> M1 = inv(S)*L*D;
>> M2 = L'*inv(S);
>> x = gmres(K,b,20,1e-08,20,M1,M2);
-----------------------------------------------------------------
4(b). Using the Matlab interface : the general interface
-----------------------------------------------------------------
- If not already in the search paths, add the directory where you installed the
interface to the search paths, e.g.
>> addpath('mi30_matlab')
>> javaaddpath('mi30_matlab')
where mi30_matlab is the directory.
[ You can add these paths permanently (see 'help pathtool')]
- A signed incomplete factorization of K is computed using hsl_mi30_precond:
>> pc = hsl_mi30_precond(K,m)
Optionally, a control structure (see Section 5) may be used:
>> pc = hsl_mi30_precond(K,m,control)
hsl_mi30_precond returns a preconditioner object pc. Information about
the constructed preconditioner is available by examining properties of
pc, see Section 6 for details.
- The preconditioner is applied to a (full) vector z and a vector y returned
as follows:
>> y = pc.apply(z)
The preconditioner may be applied multiple times.
- The memory associated with the preconditioner is freed when the object
ceases to exist. To force this to occur, the clear command may be used:
>> clear pc
-----------------------------------------------------------------
4(b)(i). Using the Matlab interface: general interface example
-----------------------------------------------------------------
Below is an example of using the Matlab interface with default controls:
>> B = sparse(gallery('uniformdata',10,0));
>> K = [speye(10) B; B' sparse(10,10)];
>> pc = hsl_mi30_precond(K,10);
>> z = ones(length(K),1);
>> y = pc.apply(z);
>> clear pc;
The interface may also be used in conjunction with an iterative solver:
>> B = sparse(gallery('uniformdata',10,0));
>> K = [speye(10) B; B' sparse(10,10)];
>> b = ones(length(K),1);
>> pc = hsl_mi30_precond(K,10);
>> x = gmres(K,b,20,1e-08,20,@pc.apply);
>> clear pc;
-----------------------------------------------------------------
5. 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. Further details of the control
components are given in the Fortran documentation.
control.lsize and control.rsize : control the number of entries in L and R.
The number of fill entries in each column of L is at most lsize and the
number of entries in each column of R is at most rsize.
In general, increasing lsize and/or rsize increases the
quality of the preconditioner but at the cost of more work (slower
factorize and precondition times) and the use of more memory.
The default values are 20.
control.alpha1 : holds the initial shift alpha1. The default value is 0.
control.alpha2 : holds the initial shift alpha2. The default value is 0.
control.iorder : controls whether K is preordered before the factorization.
Options for the general interface are:
<=0 : no ordering.
1 : reverse Cuthill-McKee (RCM) ordering.
2 : approximate minimum degree (AMD) ordering.
4 : the rows are ordered by ascending degree.
5 : METIS (nested dissection) ordering.
6 : Sloan profile reduction ordering.
All other values are treated as 6. The default value is 6.
For the hsl_mi30_sichol interface, no preordering is offered (that is, all
values are treated as 0).
control.iscale : controls whether K is to be scaled before the factorization.
Options are:
<=0 : no scaling.
1 : scaling using l_2-norm of the columns of K.
2 : MC77 (1 iteration in inf norm, 3 iteration in one norm).
3 : MC64 Bipartite matching based scaling.
4 : Diagonal scaling.
All other values are treated as 1. The default value is 1.
control.lowalpha1 : controls the choice of the first value of the
shift alpha1 in the event of factorization breakdown.
The default value is 0.001.
control.lowalpha2 : controls the choice of the first value of the
shift alpha2 in the event of factorization breakdown.
The default value is 0.001.
control.maxshift : controls the maximum number of times the shift can be
decreased after a successful factorization with a positive shift.
The default value is 3.
control.shift_factor : controls how rapidly the shift is increased after
a breakdown.
The default value is 2.
control.shift_factor2 : controls how rapidly the shift is decreased after
a successful factorization with a positive shift.
The default value is 4.
control.small : any pivot whose modulus is less than small is treated as 0.0
and, if such a pivot is encountered, the factorization breaks down, the
shift is increased and the factorization restarted.
The default value is 1e-20.
control.tau1 and control.tau2 : control the dropping of small entries from L and
R. Dropping entries leads to sparser factors but may reduce the quality
of the preconditioner.
The default values are 0.001 and 0.0001, respectively
-----------------------------------------------------------------
6. Information available to user
-----------------------------------------------------------------
The following properties of the computed preconditioner may be queried by the
user:
pc.alpha1 - Holds global shift alpha1 used.
pc.alpha2 - Holds global shift alpha1 used.
pc.nzl - Number of entries in the computed incomplete factor L.
pc.nzr - Number of entries in R.
pc.nshift - Number of non-zero shifts used.
pc.nrestart - Number of times the factorization was restarted.