Incomplete Factorizations Computed Via Fixed-Point Iterations
-------------------------------------------------------------

Three factorizations are available: 

  ILU
  IC (upper triangular incomplete Cholesky factor)
  Incomplete U^T D U

For each factorization, a Matlab reference code is available (using
synchronous updates):

  parilu_ref.m
  paric_ref.m
  pariutdu_ref.m

Mex versions of the factorizations are provided (both synchronous and
asynchronous updates) and can be called via these Matlab driver functions:

  parilu.m
  paric.m
  pariutdu.m

The mex codes are implemented in the following files:

  parilu_sync_mex.c
  parilu_async_mex.c
  paric_sync_mex.c
  paric_async_mex.c
  pariutdu_sync_mex.c
  pariutdu_async_mex.c

Threshold-based incomplete factorizations are also provided (both
synchronous and asynchronous updates).  These are implemented based on the
above fixed-point incomplete factorizations for a given sparsity pattern.
Matlab driver functions are available.  They call the corresponding
mex codes.

  parilut.m
  parict.m
  pariutdut.m


Test Programs
-------------

Change "numsweeps" and "numthreads" to check that expected behavior 
is observed:

  driv_test.m

Please supply your own test matrices in these drivers and edit as necessary:

  driv_paric.m
  driv_parict.m



Initial Guesses
---------------

The computed factors will have the same sparsity pattern as the initial
guesses for those factors.  The initial guesses for the incomplete
factors have the following requirements:

  ILU:   L has a unit diagonal; U has a full diagonal (no zeros on diagonal)
  IC:    U has a full diagonal
  IUTDU: U has a unit diagonal; D is diagonal with no zeros on diagonal

It is recommended to use initial guess based on the 
Symmetric Gauss-Seidel approximation:

A \approx L_A D_A^{-1} U_A = (L_A D_A^{-1}) D_A (D_A^{-1} U_A)

where

L_A = tril(A)
D_A = diag(A)
U_A = triu(A)

thus the initial guesses are:

LU:    L = L_A D_A^{-1},   U = U_A
LL^T:  L = L_A D_A^{-1/2}
LDU^T: L = L_A D_A^{-1},   D = D_A,   U = D_A^{-1} U_A

If the matrix A has been scaled such that its diagonal is all ones,
then the standard initial guess is to use the triangular parts of A
as the initial guess.  To be clear, however, it is not necessary
to scale the matrix if the above recommended initial guesses are used.

The sparsity pattern of the initial guess does not have to be related
to the pattern of matrix A.  To specify explicit zeros in the initial
guess, enter very small values (e.g., machine epsilon).


Special Note on Using Incomplete U^T D U
----------------------------------------

The IUTDU version may be useful for SPD matrices to allow the pivots to
become negative during intermediate sweeps.  It may also be useful for
symmetric but not positive definite problems.

The IUTDU version does not require the matrix to be scaled before a
reasonable initial guess can be constructed, although scaling may still
be needed to aid convergence.  The following constructs an initial guess
U and D for the factorization of a matrix A,

        U_a = triu(A);
        D_a = diag(diag(A));
        U = D_a\U_a; % has unit diagonal
        D = D_a;

which is based on the approximation A \approx U_a^T inv(D_a) U_a.
Since the matrix does not need to be scaled, it may be useful for
solving sequences of related systems where one factorization is used as
the initial guess for another factorization.  In this case, it would be
inappropriate to scale the matrices in different ways.

Note that it has been observed for many problems that the IUTDU version 
does not produce as good preconditioners compared to the IC version 
for the same number of sweeps (they are not mathematically the same
for SPD matrices, even in the synchronous case).


Parallelization
---------------

Thread scheduling uses OpenMP dynamic scheduling with block size 4096.
This is hard-coded, but other choices may give better multithreaded 
performance.


Feedback and Bug Reports
------------------------

Please report bugs to:

Edmond Chow
[email protected]