Re: GSoC: Incomplete Cholesky Factorizations

[Kai Torben Ohlhus]

On Tue, Oct 8, 2013 at 11:59 AM, Alexander Barth <address@hidden> wrote:

Dear Kai,

Thank you for your work on implementing these functions!

I wanted to test the modified incomplete cholesky factorization (the ichol function). But I had a problem.

I used the itsol package from https://code.google.com/p/itsol-clone-siko1056/ and zitsol from
http://www-users.cs.umn.edu/~saad/software/ITSOL/ZITSOL_1.tar.gz with the zitsol patch as described here: http://wiki.octave.org/Building_%28Z%29ITSOL. I did not apply the itsol.patch since I assumed that your hg repository already includes it (right?).

Compilation was fine and I can use the ichol function for small matrices. However with the attached 144x144 matrix, I get the following error:

octave>> load A; L = ichol(A,struct('michol','on'));
error: ichol0: zero diagonal found.

The error is generated around the line 200:
      // compensate diagonal MIC(0)
      if (michol)
    {
      diag[k] += dropsum;
      if (diag[k] <= 0.0)
        {
          std::cout << "diag1 " << diag[k] << " " << k << " " << dropsum << std::endl;
          error ("ichol0: zero diagonal found.");
          retval (0) = octave_value (k);
          return retval;
        }
    }

Before the error, I have the following values diag[k]= -1.97638, k = 125  and dropsum = -34.6556. The matrix A is well positive-defined  (min(eig(A)) =  13.804) and symmetric.

The modified incomplete cholesky factorization did work on matlab (R2013a). The variable L_matlab contains the result from matlab in case you want to know.

Do you have some ideas on what could be wrong?

Thanks again!

Alex

Dear Alex,

Your matrix and assumptions are perfectly alright, my Modified Incomplete Cholesky (MIC) factorization not yet. Without using MIC, MATLAB and my implementation both fail. But MATLABs MIC implementation still differs from mine. I'm working on that. 

%MATLAB

L = ichol(A,struct('michol','on'));

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

ans = 0.0471

If you only have this particular matrix to factorize, you can try out diagonal compensation for now:

%MATLAB and Octave

L = ichol(A,struct('michol','off','diagcomp',0.02));

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

ans =  0.046908

According to ITSOL https://code.google.com/p/itsol-clone-siko1056/, this repository already contains the patch, but ZITSOL from http://www-users.cs.umn.edu/~saad/software/ITSOL/ZITSOL_1.tar.gz has to be patched in order to work as I intended it to. If both libraries are setup on your system, it should be possible to use them from Octave.

Thank you for your interest!

Kai