# How to produce the ILU0 or ILUT as stand-alone procedures on sparse matrices?

Mathematica uses the ILU0 procedures automatically to precondition large sparse linear systems; e.g.

LinearSolve[mat, rhs, Method -> {Krylov, Preconditioner -> ILU0}]; // Timing


I wish to have the incomplete factorization of the large sparse matrix A in a standalone procedure.

I would be grateful if someone could give some tips or written code to provide the ILU0 or ILUT approximate inverses as standalone procedures.

-
?**ILU* will reveal which internal functions do this. You might want to experiment with them and figure out their calling syntax. – Szabolcs Jun 14 '12 at 8:21
There was some discussion on this in MathGroup, but nothing too concrete showed up. – J. M. Jun 14 '12 at 8:21
@Szabolcs: it seems to suffice just giving the SparseArray[] object to SparseArraySparseMatrixILU[] ; one can then set the Method option to either "ILU0" or "ILUT" as needed. One can then use LowerTriangularize[] and UpperTriangularize[] to extract the needed factors from the compressed SparseArray[] representation. – J. M. Jun 14 '12 at 8:43
@J.M. Can you post that as an answer? I'm not familiar enough with ILU. – Szabolcs Jun 14 '12 at 8:48
@Szabolcs: Maybe after some more experimentation; it seems I got lucky with my initial examples, and extracting the factors isn't as simple as I thought it was. But the output of SparseArraySparseMatrixILU[]  is similar to the output of LUDecomposition[] as expected: the compressed matrix where the $\mathbf L$ and $\mathbf U$ factors are packed, and a permutation matrix represented as an integer permutation... – J. M. Jun 14 '12 at 8:53

Here is a way to do it:

dim = 5;
s = SparseArray[{{i_, i_} -> -2., {i_, j_} /; Abs[i - j] == 1 ->
1.}, {dim, dim}, 0.];
s[[1, All]] = s[[-1, All]] = 0.;
s[[1, 1]] = s[[-1, -1]] = 1.;
f = ConstantArray[0., {dim}];
f[[1]] = 0.; f[[-1]] = 1.;

LinearSolve[s, f]
{0., 0.25, 0.5, 0.75, 1.}


Now, we can use:

res = SparseArraySparseMatrixILU[s]


And then:

SparseArraySparseMatrixApplyILU[res, f]
{0., 0.24999999999999997, 0.49999999999999994, \
0.7499999999999999, 1.}

res[[1]] // Normal
{{1., 0., 0., 0., 0.}, {1., -2., 1., 0.,
0.}, {0., -0.5, -1.5, 1., 0.}, {0.,
0., -0.6666666666666666, -1.3333333333333335, 1.}, {0., 0.,
0., 0., 1.}}


Also,

Options[SparseArraySparseMatrixILU]
{SparseArrayFillIn -> Automatic, Method -> "ILUT",
SparseArrayPermutationTolerance -> Automatic, Tolerance -> Automatic}
SparseArraySparseMatrixILU[s, Method -> "ILUTP"]


Or,

res2 = SparseArraySparseMatrixILU[LowerTriangularize[s]]


And as a preconditioner:

pf = With[{pm = res2},
Function[arg,
SparseArraySparseMatrixApplyILU[pm, arg]
]]

LinearSolve[s, f, Method -> {"Krylov", "Preconditioner" -> pf}]
{0., 0.25000000000000006, 0.5000000000000001, 0.7500000000000001, \
1.}
`

If you think this should be documented, extended etc. then you should write a comment about that to the support at wolfram only if enough people request this type of functionality will be done. Maybe.

-
Finally… ;) For reference, is there a demo showing how the original matrix can be reconstituted from the triangular factors? – J. M. Jul 4 at 13:47