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.