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

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

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

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
>>
>> http://reference.wolfram.com/mathematica/tutorial/SomeNotesOnInternalImplementation.html
>>
>> 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.

http://reference.wolfram.com/mathematica/ref/LinearSolve.html

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}];

Q2)
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:

http://en.wikipedia.org/wiki/Conjugate_gradient_method#The_preconditioned_conjugate_gradient_method

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
function.

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

http://www.cfd-online.com/Wiki/Incomplete_LU_factorization_-_ILU

And in other Fortran code on the net. http://www.netlib.org/linalg/

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,

--Nasser