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