Sunday, September 15, 2013

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.

Monday, September 9, 2013

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



Sunday, September 1, 2013

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.

Monday, August 26, 2013

A Fortran story with Algorithm 740

Last week I was looking for an good existing implementation of the Incomplete Cholesky Factorization and finally decided for Algorithm 740 by Mark T. Jones and Paul E. Plassmann. Next to this code, I had a deeper look at many promising libraries, but found nothing I could work upon. Algorithm 740 is a Fortran implementation of the no-fill Incomplete Cholesky Factorization [IC(0) in short] for column wise stored matrices. What do you want more?

To see if this code works, I started writing an interface for Octave and tested it using two small toy problems from Quarteroni, Saleri, Gervasio: Scientific Computing with MATLAB and Octave: Third Edition, Springer 2010.

% Problem 5.1 (Hydraulic network)

A =[-0.37 0.05 0.05 0.07;
    0.05 -0.116 0 0.05;
    0.05 0 -0.116 0.05;
    0.07 0.05 0.05 -0.202];
A = -A;

% Test for positive definite
for i = 1:length (A)
    assert (det (A (1:i,1:i)) > 0)
endfor

% Prepare input for solver
D = diag (A);
A_off = sparse (tril (A - diag (D)));

[L_off, D] = JPICC_wrapper (A_off, D);
L = L_off + diag (D);

% Compute error
R = A - L * L'
norm (R,'fro')

>> norm = 0.0095646


% Problem 5.4 (Capillary networks)

A = [-1/4 1/10 1/10 0 0 0 0 0 0 0 0 0 0 0 0;
     0 -1/2 0 1/5 1/5 0 0 0 0 0 0 0 0 0 0;
     0 0 -1/2 0 0 1/5 1/5 0 0 0 0 0 0 0 0;
     0 0 0 -1 0 0 0 0.4 0.4 0 0 0 0 0 0;
     0 0 0 0 -1 0 0 0 0 0.4 0.4 0 0 0 0;
     0 0 0 0 0 -1 0 0 0 0 0 0.4 0.4 0 0;
     0 0 0 0 0 0 -1 0 0 0 0 0 0 0.4 0.4;
     0 0 0 0 0 0 0 -2 0 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 -2 0 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 -2 0 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0;
     0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0;
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2];
A = -(A + triu (A, 1)');

% Test for positive definite
for i = 1:length (A)
    assert(det(A(1:i,1:i)) > 0)
endfor

% Prepare input for solver
D = diag (A);
A_off = sparse (tril (A - diag(D)));

[L_off, D] = JPICC_wrapper (A_off, D);
L = L_off + diag(D);

% Compute error
R = A - L * L'
norm(R,'fro')

>> norm = 0.53208

The lack of a well-tested library like ITSOL makes this part of the GSoC more difficult. But MATLAB implements "only" IC(0), ICT and the modified IC option. So it will be much work left to do, but it should remain possible within the remaining weeks.

Tuesday, August 20, 2013

Little update

The repository mess

I recreated the structure of my public repository, because of non-conforming Mercurial commit messages. If you already downloaded changesets from my old repository, I recommend to strip them:

   hg strip -r MY_OLD_CHANGESET

And to get the newer ones as described here.

On my way to IC(0) and ICT

To get these two implementations I had a look at several libraries again (many of them suggested by Nir):

  • Contains a very simple implementation of IC(0) in SparseLib++/1.7/src/icpre.cc
  • Implementation just consists of coping the lower triangular part, this can be done smarter according to the publications I read so far.
  • Also looks promising. But didn't got it to work so far with the MEX interface.
  • Good start for implementing an own version, Mex-Interface with C-Code
  • A deeper look at the underlying paper (http://machinelearning.wustl.edu/mlpapers/paper_files/FineS01.pdf) revealed, that their understanding of an "incomplete Cholesky Factorization" is: "If $Q$ is positive semidefinite and singular, then it is still possible to compute an “incomplete” Cholesky factorization $GG^T$, where some columns of $G$ are zero." That is definitely not what MATLABs ichol does. 
  • Not free, but authors are requested. I started to port this function to C++ in order to test it, because it looked as one of the most promising ones.


New insights about Scilab and ILU

During the research about Incompete Cholesky factorizations, I found an interesting toolbox Spilu at Scilab. The developers of Scilab took Saads SpimplementedarseKit as basis of their seemingly well tested Fortran functions for nearly all ILUs I interfaced from ITSOL. It is definitly worth to have a look at this library, but I don't know if it was reasonable to change ITSOL to Scilabs implementation. Unfortunately Scilab doesn't seem to have any IC factorization...

News from Yousef Saad

I created another public repository with the changes Marco Atzeri (http://matzeri.altervista.org/works/itsol/) and I made to ITSOL https://code.google.com/p/itsol-clone-siko1056/. The purpose of the repository is not to maintain an own version, it is to facilitate the decision process for upstream changes by Yousef Saad. He wants to consider my changes, when Ruipeng Li is soon available again.

Sunday, August 4, 2013

Mercurial commit message cleaning

Now the midterm is over, one can already think about to make my code ready to be pushed to the official Savannah repository of GNU/Octave. So on the maintainers-mailing list there was a discussion about how to make my changesets fit for the official repository.

Again thanks to Jordi I really learned a lot, how to "manipulate" the Mercurial history. As a first step one should get Mercurial 2.7, because older versions might not support all features used in this post. The current version from the Debian repository does not suffice.

A good point to start is getting a list of all my changesets:

  hg log --user k.ohlhus

Now my first change in current development is, that I no longer want to merge the Savannah repository into mine. Instead I got to know about the rebase command from Mercurial, which sets all your nonofficial changesets on top of the history, that makes it easier to track what has changed in my repository.

The old state of my repository is very good explained in Scenario B of the rebase manual page. To fix that I had to run a mix of two commands for all my merges, e.g.:

  hg phase --draft MERGE_CHANGESET --force
  hg rebase --source MERGE_CHANGESET --dest LAST_CHANGESET_FROM_SAVANNAH
  
From now all my changes are on top. The second task of cleaning is to tidy up my commit messages, as I did not stick to the commit message guidelines of Octave. By typing

  hg histedit -r LAST_BAD_CHANGESET_OF_MINE

I can create an edit plan for my changeset commit messages and carry this plan out.


edit plan example:

mess a623d25a6b42 17121 Startet fixing the transposition problem with all ILU. S
pick 466fb4342423 17122 Interface change in itsol_util.h and further addition of
pick 5afd5d5f42d0 17123 Added MATLAB compatibility wrapper for ILU-factorization
pick 0dbc3caefb6e 17124 Integrating ZITSOL-library to Octave build system.
pick e0b86073ce45 17125 Extended ILU-Preconditioner for complex input using ZITS
pick 9436002c0a0b 17126 Updated test cases. Fixing ILUC output.
pick 9d4e24a01f42 17127 Test cases for complex case increased.
pick 1b6ac1f0991b 17128 Test cases for complex case in ilut.
pick 286aa419eb58 17129 Redirected output of ITSOL and ZITSOL from stderr to tem
pick 975f30ec1c1b 17130 Make most recent development on ILUTP public. Still not
pick 65c0db01a558 17131 Final tidy up on Matlab wrapper with tests.
pick 4244faff69d6 17132 Moved ilu.m from scripts/linear-algebra to scripts/spars
pick 31a15f423dd8 17133 removed and new build entries.

# Edit history between a623d25a6b42 and 31a15f423dd8
#
# Commands:
#  p, pick = use commit
#  e, edit = use commit, but stop for amending
#  f, fold = use commit, but fold into previous commit (combines N and N-1)
#  d, drop = remove commit from history
#  m, mess = edit message without changing commit content
#


This tedious, but necessary work has to find its way to my development repository after my journey (13th of August). Then I'll provide instructions how to get  the new changesets and how to remove the older ones. In any case I'm sorry for the inconvenience.

Tuesday, July 30, 2013

GSoC Part 2

Now follows the second part of the Google Summer of Code. Here is what I planned to do in that time.

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

TODOs: