8
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ ?*`*ILU* will reveal which internal functions do this. You might want to experiment with them and figure out their calling syntax. $\endgroup$ – Szabolcs Jun 14 '12 at 8:21
  • $\begingroup$ There was some discussion on this in MathGroup, but nothing too concrete showed up. $\endgroup$ – J. M. will be back soon Jun 14 '12 at 8:21
  • $\begingroup$ @Szabolcs: it seems to suffice just giving the SparseArray[] object to SparseArray`SparseMatrixILU[] ; 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. $\endgroup$ – J. M. will be back soon Jun 14 '12 at 8:43
  • $\begingroup$ @J.M. Can you post that as an answer? I'm not familiar enough with ILU. $\endgroup$ – Szabolcs Jun 14 '12 at 8:48
  • $\begingroup$ @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 SparseArray`SparseMatrixILU[] 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... $\endgroup$ – J. M. will be back soon Jun 14 '12 at 8:53
7
+100
$\begingroup$

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 = SparseArray`SparseMatrixILU[s]

And then:

SparseArray`SparseMatrixApplyILU[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[SparseArray`SparseMatrixILU]
{SparseArray`FillIn -> Automatic, Method -> "ILUT", 
 SparseArray`PermutationTolerance -> Automatic, Tolerance -> Automatic}
SparseArray`SparseMatrixILU[s, Method -> "ILUTP"]

Or,

res2 = SparseArray`SparseMatrixILU[LowerTriangularize[s]]

And as a preconditioner:

pf = With[{pm = res2},
  Function[arg,
   SparseArray`SparseMatrixApplyILU[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.

$\endgroup$
  • 2
    $\begingroup$ Finally… ;) For reference, is there a demo showing how the original matrix can be reconstituted from the triangular factors? $\endgroup$ – J. M. will be back soon Jul 4 '16 at 13:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.