Hi,

How I can efficiently invert a sparse lower triangular matrix ?
My matrixes are large (size up to a million) and very sparse.

The context
I try to implement an Incomplete Cholesky Preconditionner to feed a Preconditionned Conjugate Gradient solver. The Incomplete Cholesky Preconditionner is explained here http://en.wikipedia.org/wiki/Incomplete ... torization

1) I have a matrix A, ok. Cholesky magic, I got a matrix L which is lower triangular.
2) I generate a matrix K from L and A, K is also lower triangular.
3) I invert K
4) Profit !

I miss the step 3)

The code
Here is is my attempt at doing the stuff described in Wikipedia

Code: Select all

// Go go Cholesky, give me L
Eigen::SparseLLT< Eigen::SparseMatrix<double>, Eigen::Cholmod > lLLT;
lLLT.compute(lA);

// Generate K from L and A
Eigen::SparseMatrix<double> lK(nbCells(), nbCells());
lK.reserve(lA.nonZeros());

for(int j = 0; j < lA.outerSize(); ++j)
  for(Eigen::SparseMatrix<double>::InnerIterator itA(lA, j); itA; ++itA)
    if (itA.row() <= itA.col())
      lK.insert(itA.row(), itA.col()) = lLLT.matrixL().coeff(itA.row(), itA.col());

lK.finalize();

// Compute the preconditionner from K
lInvK = ???
lPreCond = lInvK * lInvK.adjoint();


If your matrix is sparse lower triangular, and you want to compute X = A^-1 B simply do:

X = A.triangularView<Lower>.solve(B);

No need to explicitly invert it. Moreover, the inverse of a sparse matrix can be very dense, so it is really not recommended to try to compute its inverse.

Thank you, indeed, makes sense !

So, for those who came here with the same question, here it is the code I used, to compute X = M^-1 * B:

Code: Select all

X = M.triangularView<Eigen::Lower>().solve(B);


Note than you can do that too:

Code: Select all

Eigen::SparseTriangularView<Eigen::SparseMatrix<double>, Eigen::Lower> MView = M.triangularView<Eigen::Lower>();

X = MView.solve(B);