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