----------------------------------------------------------------- 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.