Cholesky decomposition of large sparse matrices in Java

Question

I want to do Cholesky decomposition of large sparse matrices in Java. Currently I'm using the Parallel Colt library class SparseDoubleCholeskyDecomposition but it's much slower than using my code I wrote in C for dense matrices which I use in java with JNI.

For example for 5570x5570 matrices with a non-zero density of 0.25% with SparseDoubleCholeskyDecomposition takes 26.6 seconds to factor and my own code for the same matrix using dense storage only takes 1.12 seconds. However, if I set the density to 0.025% then the colt library only takes 0.13 seconds.

I have also written my own sparse matrix Cholesky decomposition in C using compressed row storage and OpenMP and it's quite a bit faster than SparseDoubleCholeskyDecomposition as well but still slower than my algorithm using dense storage. I could probably optimize it further but in any case the Cholesky factors are dense so it neither solves the speed nor the storage problem.

But I want times in the milliseconds and I eventually need to scale to over 10000x10000 matrices so even my own dense code for dense matrices will become too slow and use too much memory.

I have reason to believe the Cholesky decomposition of sparse matrices can be done much faster and use less memory. Maybe I need a better Java BLAS library? Is there a library for Java which does Cholesky decomposition for sparse matrices efficiently? Maybe I'm not using Parallel Colt optimally?

import cern.colt.matrix.tdouble.algo.decomposition.SparseDoubleCholeskyDecomposition;
import cern.colt.matrix.tdouble.impl.*;
import java.util.Random;

public class Cholesky {
    static SparseRCDoubleMatrix2D create_sparse(double[] a, int n, double r) {
        SparseRCDoubleMatrix2D result = new SparseRCDoubleMatrix2D(n, n);
        Random rand = new Random();
        for(int i=0; i<n; i++) {
            for(int j=i; j<n; j++) {
                double c = rand.nextDouble();
                double element = c<r ? rand.nextDouble() : 0;
                a[i*n+j] = element;
                a[j*n+i] = element;
                if (element != 0) {
                    result.setQuick(i, j, element);
                    result.setQuick(j, i, element);
                }
            }
        }
        for(int i=0; i<n; i++) {
            a[i * n + i] += n;
            result.setQuick(i,i, result.getQuick(i,i) + n);
        }
        return result;
    }

    public static void main(String[] args) {
        int n = 5570;
        //int n = 2048;
        double[] a = new double[n*n];
        SparseRCDoubleMatrix2D sparseMatrix = create_sparse(a, n, 0.0025);

        long startTime, endTime, duration;

        startTime = System.nanoTime();
        SparseDoubleCholeskyDecomposition sparseDoubleCholeskyDecomposition = new SparseDoubleCholeskyDecomposition(sparseMatrix, 0);
        endTime = System.nanoTime();
        duration = (endTime - startTime);
        System.out.printf("colt time construct %.2f s\n", 1E-9*duration);
    }
}

Impeller