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

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

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

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

AbsoluteTiming[
  solutionILUTP =
    LinearSolve[stiffness, load,
     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
>
>
>
>