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