MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: A function to do incomplete LU decomposition with a drop tolerance?

  • To: mathgroup at
  • Subject: [mg123313] Re: A function to do incomplete LU decomposition with a drop tolerance?
  • From: "Nasser M. Abbasi" <nma at>
  • Date: Thu, 1 Dec 2011 05:51:24 -0500 (EST)
  • Delivered-to:

On 11/30/2011 5:51 AM, Oliver Ruebenkoenig wrote:
> On Wed, 30 Nov 2011, Nasser M. Abbasi wrote:
>> Version 8.0.4
>> I need to use ILU as a preconditioner for a conjugate gradient solver
>> I am writing in Mathematica. But looking around, I am not able to
>> find one in Mathematica. Only reference I found is here
>> Where it say in the above the following:
>> "For sparse arrays, LinearSolve uses UMFPACK multifrontal direct solver
>> methods and with Method uses Krylov iterative methods preconditioned
>> by an incomplete LU factorization"
>> But how does a user access this ILU from Mathematica?
>> I do have the A matrix (sparse).

> Here are some ideas:
> AbsoluteTiming[
>    solutionGMRES =
>      LinearSolve[stiffness, load,
>       Method ->  {"Krylov", Method ->  "GMRES", "Preconditioner" ->  None,
>         Tolerance ->  tol, MaxIterations ->  100}];]

<SNIP more nice examples>

Thanks Oliver for the examples, they are useful.

I have 2 follow up questions:

Q1. Are these documented somewhere? looking in the documentation center, not
able to find more information. Only mention of "Krylov" is one line
in the help for LinearSolve.

For example, how does one know what ILU0,ILU1,ILUT,ILUTP etc.. and
learn more about the API?  I am sure there is a document somewhere
that goes into details, but I am not finding it now.

> What you will need however, is something different. If you want to go for
> iterative solves you'll need bandwidth reordering.
> Needs["GraphUtilities`"]
> {r, c} = MinimumBandwidthOrdering[mat["PatternArray"]];
> MatrixPlot[ mat[[r, c]] ]

I am not sure I understand this last part. How does the above fit
into the examples you shown before? How/when do I use that? When I tried:

{row, col} = MinimumBandwidthOrdering[stiffnessMatrix["PatternArray"]];
sol = LinearSolve[stiffnessMatrix[[row, col]], loadVector,
   Method -> {"Krylov", Method -> "GMRES", "Preconditioner" -> None,
     Tolerance -> 0.001, MaxIterations -> 100}];

Then the solution did not look right, compared to just writing

sol = LinearSolve[stiffnessMatrix, loadVector,
   Method -> {"Krylov", Method -> "GMRES", "Preconditioner" -> None,
     Tolerance -> 0.001, MaxIterations -> 100}];

But what I really was looking for is to do incomplete LU decomposition
on my A matrix as a stand alone operation, as I wrote the preconditioned
conjugate gradient solver in Mathematica, described here:

And you see in the above algorithm, what is called the 'M' matrix, is the
preconditioner matrix. This is what I wanted to generate in Mathematica.

One choice for a preconditioner matrix is to use ILU on A, then save
result in what is called M above, and then use it in the iterations
as decribed above.

I implemented the above conjuagte gradient myself and not using LinearSolver.
(This is for learning, and will be in a demo later.) Actually I did not
know that one can use "GMRES" in Mathematica before now.

So my question is, it would be nice if Mathematica had a direct API
to allow one to do this ILU decomposition on a matrix as a stand alone

I can try to implement ILU myself, the algorithm is described here

And in other Fortran code on the net.

But thought to ask first here before I do that.

>> Also, I see still that the Incomplete Cholesky decomposition
>> does not have a drop tolerance parameter. I hope in version 9
>> it will.
> Hm, while I agree that this would be useful - there are bigger fish to
> fry. This a personal opinion and I do not know if this will be in V9.
> Oliver

Ok, but just you know, your market competitor in this area, does has
this option on their chol() implementation and has ilu() for many years
now :)

thank you,


  • Prev by Date: Re: RE: Re: How to simplify ArcSin formula
  • Next by Date: Re: FindShortestTour Function- Error
  • Previous by thread: Re: gives different result compared to 1/Diagonal[Normal@A] when A is sparse
  • Next by thread: Re: A function to do incomplete LU decomposition with a drop tolerance?