That's all folks?

Definitely not! Today is the GSoC "pencils down" day and I want to summarize the state of the project. You can get everything so far from here.

Goals for the Midterm 2013-07-29 (Monday) :

• Integration of ITSOL completed (iluk.cc, ilut.cc, ilutp.cc, iluc.cc), callable from Octave
• Wrapper file for MATLAB compatibility (ilu.m)
• Small test cases for incomplete LU-factorization, showing their correctness
• Documentation of incomplete LU-factorization types
• Comprehensive test cases for all types of incomplete LU-factorization types (also  proofing compatibility to MATLAB)
• (optional) Complex version for the incomplete-factorizations

Goals for the Final Evaluation 2013-09-16 (Monday, "Suggested 'pencils down' date") :

• Implementations for IC(0) and ICT completed
• Wrapper file for MATLAB compatibility (ichol.m) finished
• Small test cases for all incomplete-factorizations
• Documentation of incomplete LU-factorization types
• Comprehensive test cases for all types of incomplete-factorization types (also  proofing compatibility to MATLAB)
• (optional) Complex version for the incomplete-factorizations

What there is left to do?

• ILUTP must get to work
• Make ITSOL a dependency for Octave (there are still some flaws in the official distribution, that make it difficult to use)
• IC routines (ichol0jp and icholt) still suffer from Fortran indexing and missing complex versions
• Implement the Modified version for all preconditioners.
• Everything has to find it's way to the official Octave development repository.

One week left

The final week is coming closer. Here is the current state of my project:

Goals for the Final Evaluation 2013-09-16 (Monday, "Suggested 'pencils down' date") :

• Implementations for IC(0) and ICT completed
• Wrapper file for MATLAB compatibility (ichol.m) finished
• Small test cases for all incomplete-factorizations (ICT missing)
• Documentation of incomplete LU-factorization types (ICT missing)
• Comprehensive test cases for all types of incomplete-factorization types (also  proofing compatibility to MATLAB) (ICT missing)
• (optional) Complex version for the incomplete-factorizations (ICT missing)

This week ICT and the remaining documentation and test cases will be on my schedule. Last week I implemented a Modified Incomplete Factorization (MIC) option for IC(0), which unfortunately is not 100% compatible to MATLAB. As far as I understood the MIC (see Saad (chapter 10.3.5) for example) you add after a completed step the dropped elements (here marked red) to the main diagonal. MATLAB seems to handle this differently. If anyone knows more, let me know.

A = [ 0.37, -0.05,  -0.05,  -0.07;

     -0.05,  0.116,  0.0,   -0.05;

     -0.05,  0.0,    0.116, -0.05;

     -0.07, -0.05,  -0.05,   0.202];

Normal Cholesky Factorization

L1 = chol (A, 'lower');

full(L1)

ans =

   0.608276253029822                   0                   0                   0

  -0.082199493652679   0.330519656364403                   0                   0

  -0.082199493652679  -0.020442828820163   0.329886850288205                   0

  -0.115079291113750  -0.179896893617439  -0.191390050272693   0.346068942669187

Incomplete Cholesky Factorization no MIC

opts.type = 'nofill';

opts.michol = 'off';

L2 = ichol (sparse (A), opts);

full(L2)

ans =

   0.608276253029822                   0                   0                   0

  -0.082199493652679   0.330519656364403                   0                   0

  -0.082199493652679                   0   0.330519656364403                   0

  -0.115079291113750  -0.179896893617439  -0.179896893617439   0.352180311900522

Incomplete Cholesky Factorization with MIC

opts.michol = 'on';

L3 = ichol (sparse (A), opts);

full(L3)

ans =

   0.608276253029822                   0                   0                   0

  -0.082199493652679   0.320135106613577                   0                   0

  -0.082199493652679                   0   0.320135106613577                   0

  -0.115079291113750  -0.185732393077497  -0.185732393077497   0.346068942669187

The way the dropped element is computed, no problem.

drop = (A(3,2) - L2(2,1) * L2(3,1)) / L2 (2,2)

drop =

  -0.020442828820163

Expected compensated diagonal value

D_2_2 = L2(2,2) + drop % * x

D_2_2 =

   0.310076827544240

Some factor, the results differ with each other

x = full(L3 (2,2) - L2(2,2)) / drop

x =

   0.507980076641026

But the main diagonal D_2_2 must be computed first. The off diagonal can be computed with it.

L_4_2 = (A(4,2) - L3(2,1) * L3(4,1)) / L3 (2,2)

L_4_2 =

  -0.185732393077497

The complex input handling works compatible to MATALB. An example:

A = [ 0.37, -0.05,         -0.05,  -0.07;

     -0.05,  0.116,         0.0,   -0.05 + 0.05i;

     -0.05,  0.0,           0.116, -0.05;

     -0.07, -0.05 - 0.05i, -0.05,   0.202];

opts.type = 'nofill';

opts.michol = 'off';

L = ichol (sparse (A), opts);

norm (A - L * ctranspose (L), 'fro') / norm (A, 'fro')

ans =  0.0195288512358234

Lost in translation with nice results

Last week I spend reimplementing the IC(0) Fortran routines by Mark T. Jones and Paul E. Plassmann (see my last blog post) in C++. This effort was done, because

• Octave is naively C++. Thus it is possible to adapt the code to more sophisticated data types (especially I think in complex sparse data types) so it will be able to handle more general input.
• I am much more familiar with C++. This makes it easier for me to extend the code with features like the Modified Incomplete Cholesky Factorization (MIC) option and diagonal shifting.

The result is accessible from my repository:

• ichol0.cc - The standard IC(0) routine
• ichol0jp.cc - IC(0) with Jones and Plassmann strategy (preserve number of non-zeros per row/column, but take the largest ones)

For testing the implementation, I took the two examples from the last post and added one example test case from http://www.mathworks.com/help/matlab/ref/ichol.html. The m-files are available from here.

Octaves output

Test Problem 1

  err_ichol0   =  0.0197360236512918

  err_ichol0jp =  0.0197548896386716

Test Problem 2

  err_ichol0   =  0.0840790133633081

  err_ichol0jp =  0.0846790338305638

Test Problem 3

  err_ichol0   =  0.0915989875674983

  err_ichol0jp =  0.0926062848306512

MATLAB (R2012b) output

Test Problem 1

  err_ichol = 0.019736023651292

Test Problem 2

  err_ichol = 0.084079013363308

Test Problem 3

  err_ichol = 0.091598987567498

Interesting was to observe, that the results of the standard IC(0) match very well with those of MATLAB. The Jones and Plassmann strategy does not differ that much from the standard IC(0), but according to their paper, it should speed up the preconditioned solution. So it will be offered as an additional options.type='jp' for example.