Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
1  
?*`*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 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. – 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 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... – 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 = 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.

share|improve this answer
1  
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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