Tuesday, 18 March 2014

ILU0 implementation.



I have implemented the most basic function ilu0 that do the incomplete LU decomposition with 0-fill. I have drawn a table with the execution times of that algorithm using Matlab, my version and using the code that Kai implemented last year using ITSOL. The function right now can be used with real and complex numbers and has the milu='row' option implemented.

The table shows for a NxN sparse matrix with a number of non-zero elements NNZ, the time  of execution (tic - toc was used).

N-NNZ ILU0-mine ILU0-Matlab ILU0-ITSOL
50 - 683 0.000055 s 0.000065 s 0.00085 s
400 - 72435 0.027 s 0.024 s 0.04 s
2000 - 1805571 3.35 s 3.25 s 4.88 s
5000 - 6482839 14.2 s 14.5 s 22.75 s

It can be seen that the implementation using ITSOL is the slowest. Maybe just because the overhead of copying and translating data back and forth between Octave and ITSOL. Between my version and Matlab's there is almost no difference.


Here you can download the code (ilu0.cc). It does not have any test nor documentation written yet.
 

Project goals

My intention is upgrading  some functions related with sparse matrices  so they become compliant with Matlab  and implement others that are not present in Octave right now. This is the list plus some comments about how I expect to do things.

  1. Implement  sprand/sprandn with the 4rth parameter.(link) (Patch accepted)
  2. sprandsym (implement arguments needed to be compliant with Matlab)
  3. lsqr 
  4. minres
  5. ilu (complete last year GSOC project)
  6. ichol (complete last year GSOC project)
 1: I have already submitted a patch but maybe some modifications are needed to improve performance. Anyway I have mailed Rik (the last one that modified the source code of that functions) and has told me that he will give it a look in a couple of weeks.

2: I have not yet investigated enough to give a road map for implementing that function.

3 & 4: For implementing that functions I found those links:

In this website there are codes written in several programming languages that implement the algorithms (the authors are the ones that written the paper that Matlab gives as reference in their documentation for both functions). I have emailed professor Michael Saunders about adapting them into Octave versions and he answered me that I am welcome to do while I respect the license (CPL or BSD licenses). They are meant to be very unrestrictive but I need to get informed about the compatibility with GPL3.

5: Here comes the big one. That function has a big chunk of options and the last year was almost implemented by Kai Torben as his GSOC project. He interfaced Octave with ITSOL/ZITSOL libraries but in the end there were some issues with that approach:

  • ILUTP algorithm did not work 
  • He had to patch the library to get things work!
  • modified versions of algorithms ("milu" option) were not implemented in the libraries
  • That "ugly" scenario lead to finally not being able to include ITSOL as a dependency with Octave. Bottom line, the integration with the development repository could not be achieved.
My approach: I will write from scratch all the functions needed as oct-files (ILUTP, ILU0, ILUC and ILUT) for real and complex numbers implementing the algorithms described by Yousef Saad in his book "Iterative methods for sparse linear systems Ed. 2". I can use some of the code that Kai wrote, mostly the tests and the m-file "ilu.m" that glue together all the functions.

6: That function need less work since it is almost all implemented. Just some complex implementations are missing and the modified version of the algorithms too. Since Kai implemented those functions from scratch there are no dependency problems. That is nice.

Update 1: Kai has commented me that there are some license issues on the FORTRAN code he used (see here). So ichol related functions need to be implemented from scratch.


Update 2: Following Kai's recommendations I will focus on ilu/ichol functions for the GSoC period and just in case there is some time left go ahead with the other functions.

Sunday, 16 March 2014

Introducing myself

My name is Eduardo (edu159), and that blog has the purpose of tracking the state of my project with Octave during the GSOC2014 program  (if I become selected) .

You can visit my public profile at the Octave wiki here:


Feedback would be welcome. Feel free to comment :).

Eduardo