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.
Posted by at No comments:

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 semideļ¬nite 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.
Posted by at No comments:

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.
Posted by at No comments:
Subscribe to: Posts (Atom)