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

• To: mathgroup at smc.vnet.net
• Subject: [mg123302] Re: A function to do incomplete LU decomposition with a drop tolerance?
• From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
• Date: Wed, 30 Nov 2011 07:50:15 -0500 (EST)
• Delivered-to: l-mathgroup@mail-archive0.wolfram.com

```
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 =
Method -> {"Krylov", Method -> "GMRES", "Preconditioner" -> None,
Tolerance -> tol, MaxIterations -> 100}];]

AbsoluteTiming[
solutionILU0 =
Method -> {"Krylov", Method -> "GMRES",
"Preconditioner" -> "ILU0", Tolerance -> tol,
MaxIterations -> 100}];]

AbsoluteTiming[
solutionILUT =
Method -> {"Krylov", Method -> "GMRES",
"Preconditioner" -> "ILUT", Tolerance -> tol,
MaxIterations -> 100}];]

Timing[solutionILUT2 =
Method -> {"Krylov", Method -> "GMRES",
"Preconditioner" -> {"ILUT", "Tolerance" -> 10^-8,
"FillIn" -> 8}, Tolerance -> tol, MaxIterations -> 100}];]

AbsoluteTiming[
solutionILUTP =
Method -> {"Krylov", Method -> "GMRES",
"Preconditioner" -> "ILUTP", Tolerance -> tol,
MaxIterations -> 100}];]

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

where mat is a SparseArray (* you do not need the PatternArry, but if you
are tight on memory - which I assume you are (otherwise use a direct
solver!) - can save some memory for coupled problems. If you have this
kind of problem you'd have to send some more detail *)

Alternative:

p3 = SparseArray`ApproximateMinimumDegree[mat];
p4 = SparseArray`NestedDissection[mat];

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

>
> thanks,
> --Nasser
>
>
>
>

```

• Prev by Date: Re: gives different result compared to 1/Diagonal[Normal@A] when A is sparse
• Previous by thread: A function to do incomplete LU decomposition with a drop tolerance?
• Next by thread: Problems with a first order differential equation