This project has moved. For the latest updates, please go here.

MILU(0) preconditioner

Sep 20, 2013 at 12:12 PM
Edited Sep 20, 2013 at 12:17 PM
I wrote a milu(0) preconditioner which is a lot faster than your incomplete LU implementation and can reduce the number of iterations considerably (compared to standard ilu(0)).

Feel free to include it in your project: MILU0-TestProject.rar

Some benchmarks:
fidap004.mtx: 1601x1601 (31837 nz)
        Milu0   IncompleteLU
Init :      *          964ms
Solve:      *          312ms

fidap007.mtx: 1633x1633 (46570 nz)
        Milu0   IncompleteLU
Init :    5ms         1458ms
Solve:    2ms          322ms

fidap011.mtx: 16614x16614 (1091362 nz)
        Milu0   IncompleteLU
Init :  246ms       345174ms
Solve:   48ms        38161ms

utm1700a.mtx: 1700x1700 (21313 nz)
        Milu0   IncompleteLU
Init :    1ms          564ms
Solve:    0ms          305ms

utm5940.mtx: 5940x5940 (83842 nz)
        Milu0   IncompleteLU
Init :    6ms         7730ms
Solve:    3ms         3825ms
* = Error: Zero pivot encountered on row 24 during ILU process

The error is conforming to what Matlab does. Though your IncompleteLU works just fine for this matrix, I wonder if it does anything good regarding iteration speed!?

Sep 21, 2013 at 10:29 AM
Excellent, thanks! Yes, I'll include it into the library.

I've noticed that the existing ILU(0) and ILUTP implementations do not actually leverage the sparse storage format at all yet, which may explain at least a major part of their slowness (but of course has no effect on the number of iterations).

Sep 21, 2013 at 10:36 AM
PS: how should I refer to you for credits at least in the release notes? Just "wo80"?

Sep 23, 2013 at 8:31 AM
Edited Sep 23, 2013 at 9:06 AM
You can use my full name there (Christian Woltering).

BTW, the code is part of a major rewrite of the sparse matrix operations (with focus on performance). Most of the code comes from Saad's SPARSKIT. If you're interested, I can upload it as well...
Sep 23, 2013 at 8:02 PM
Edited Sep 23, 2013 at 8:03 PM
For completeness: the complex version works just the same. Replace all 'double' with 'Complex' and that's it.
Sep 23, 2013 at 8:04 PM

Yes, that would be great!

Btw, I've noticed in MILU(0) you had converted to the common CSR format where the number of non-zero values is stored in an additional field of the row pointers array. I actually had planned switching to this representation for a while and finally did the change; now in master. Beside of the simpler loop boundaries, this format is also directly supported by MKL in case we want to extend the native providers to sparse routines in the future.

Sep 23, 2013 at 8:05 PM
Complex: thanks - I'll make it available for all 4 types
Sep 24, 2013 at 3:09 PM
Finally the fix for CSR storage! Really like what you're doing for v3 ...

But you have introduced a bug in the MILU0 code. Arrays for MSR storage have to be (nonzeros + 1) long! So
_alu = new double[ia.Length];
_jlu = new int[ia.Length];
has to be
_alu = new double[ia[n] + 1];
_jlu = new int[ia[n] + 1];
I will upload the other stuff I have later on this week.
Sep 24, 2013 at 6:05 PM
doh - maybe I shoudn't do any refactoring at 2am. Thanks for pointing this out!