-----------------------------------------------------------------
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
-----------------------------------------------------------------
HSL_MI35 computes an incomplete Cholesky (IC) factorization preconditioner that may
be used in the solution of weighted sparse linear least squares problems.
Given an m x n (m >= n) sparse matrix A and an m x m diagonal matrix of weights W (which
may be the identity matrix), HSL_MI35 computes an incomplete factorization of the normal matrix
C = A' W*W A
The IC factorization is based on a matrix decomposition of the form
(L+R) (L+R)' + E,
where the incomplete Cholesky factor L is a lower triangular matrix with
positive diagonal entries. The matrix R is a strictly lower triangular 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 C may be optionally reordered, scaled and,
if necessary, shifted to avoid breakdown of the IC factorization. The
computed lower triangular matrix L is such that LL' is an incomplete factorization of
SQ' C QS + alpha I,
where Q is a permutation matrix, S is a diagonal scaling matrix and alpha is a non-negative shift.
The incomplete factorization may be used for preconditioning when solving the n x n
normal equations
A' W*W Ax = A' W*W b.
-----------------------------------------------------------------
2. Requirements
-----------------------------------------------------------------
These instructions are for linux-based systems (both 32-bit and 64-bit machines
are supported).
Requirements:
- Matlab 2011a or more recent version
- The gfortran compiler
-----------------------------------------------------------------
3. Installation
-----------------------------------------------------------------
Instructions for installing the Matlab Interface for HSL_MI35 are located in
INSTALL.
-----------------------------------------------------------------
4. Using the Matlab interface
-----------------------------------------------------------------
Two interfaces are offered:
(a) hsl_mi35_ls_ichol may be used as a replacement for applying the Matlab ichol to the
normal matrix C
(b) general interface that offers greater flexibility and more options than (a).
In particular, A may optionally be preordered (this is not offered by (a)).
Note that if (b) is used with default settings, A is preordered (see
control.iorder in Section 6) and, in this case, (a) and (b) will return
different results; the results will be the same if control.iorder is set to 0
(with all other components of control set to their defaults).
----------------------------------------------------------------------------------------
4(a). Using the Matlab interface : the replacement for ichol applied to normal equations
----------------------------------------------------------------------------------------
- If not already in the search paths, add the directory where you installed the
interface to the search paths, e.g.
>> addpath('mi35_matlab')
>> javaaddpath('mi35_matlab')
where mi35_matlab is the directory.
[ You can add these paths permanently (see 'help pathtool')]
- To compute an incomplete Cholesky factorization
>> L = hsl_mi35_ls_ichol(A)
Options may be specified using a control structure (see Section 5)
>> L = hsl_mi35_ls_ichol(A,control)
-----------------------------------------------------------------
4(a)(i). Using the Matlab interface: replacement for ichol: example
-----------------------------------------------------------------
Below is an example of using the interface that is a replacement for ichol
(with default controls) in conjunction with an iterative solver:
>> A = sprand(2000,1000,0.01);
>> b = ones(size(A,1),1);
>> L = hsl_mi35_ls_ichol(A);
>> x = lsqr(A,b,1e-8,100,L');
-----------------------------------------------------------------
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('mi35_matlab')
>> javaaddpath('mi35_matlab')
where mi35_matlab is the directory.
[ You can add these paths permanently (see 'help pathtool')]
- An incomplete factorization of A'A is computed using hsl_mi35_ls_precond:
>> pc = hsl_mi35_ls_precond(A)
Optionally, a control structure (see Section 5) may be used:
>> pc = hsl_mi35_ls_precond(A,control)
hsl_mi35_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 incomplete factors L and L' can be applied to a (full) vector z and
a vector y is returned as follows:
- to solve the system L y = z:
>> y = pc.solve('N', z);
- to solve the system L' y = z
>> y = pc.solve('T', z);
- 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:
>> A = sprand(2000,1000,0.01);
>> pc = hsl_mi35_ls_precond(A);
>> z = ones(size(A,2),1);
>> y = pc.apply(z);
>> clear pc;
or
>> A = sprand(2000,1000,0.01);
>> pc = hsl_mi35_ls_precond(A);
>> z = ones(size(A,2),1);
>> y = pc.solve('T', pc.solve('N', z));
>> clear pc;
The interface may also be used in conjunction with an iterative solver:
(You need to copy the following script into a matlab script file and run
it since Matlab Command Window does not allow defining a function)
A = sprand(2000,1000,0.01);
b = ones(size(A,1),1);
pc = hsl_mi35_ls_precond(A);
prec =@(x,t) splitprec(x,t,pc);
x = lsqr(A,b,1e-8,100,prec);
clear pc;
function y = splitprec(x,t,pc)
if(strcmp(t,'notransp'))
y = pc.solve('T',x); % we want to apply L' -- the preconditioned system is A L\^{-T} x = b
else
y = pc.solve('N',x);
end
end
-----------------------------------------------------------------
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.alpha : holds the initial shift alpha. The default value is 0.
control.iorder : controls whether A 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 (will need to be compiled
against the METIS library).
6 : Sloan profile reduction ordering.
All other values are treated as 6. The default value is 6.
For the ichol replacement interface, no preordering is offered (that is, all
values are treated as 0).
control.iscale : controls whether A is to be scaled before the factorization.
Options are:
<=0 : no scaling.
1 : scaling using l_2-norm of the columns of A.
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.lowalpha : controls the choice of the first non-zero
shift 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.alpha - Holds global shift 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.