Sparse Incomplete Cholesky

Question

I'm looking for an efficient, multicore, library to do incomplete cholesky (possibly modified). Many ILU code exists, but I can't find much about IC except in PETSC or Pastix. Could some of you drop me any library name ?

Thank you ! Tom

Answer 1

First of all, if you have an incomplete LU factorization

$A \approx LU$,

you can write the upper triangular factor $U$ as $U = DR$ where $D$ is diagonal and $R$ is unit right-triangular. If the matrix $A$ is symmetric positive definite, then $R = L^\top$. You can then easily modify the incomplete LU factorization to

$A = LDL^\top$,

or to an honest incomplete Cholesky factorization. The matrices you do find should be transposes of each other to within machine precision. While this doesn't answer your question of finding a library that specifically does parallel incomplete Cholesky, it will, with minor modifications, get you what you need.

The Euclid library is pretty popular for parallel ILU; PETSc interfaces to it. However, near as I can tell, it only does ILU and not IC, hence my digression above. I was going to recommend that you look at Hypre, but upon looking through their user manual, they tell you to just use Euclid instead.

The other one that comes to mind is pARMS (parallel algebraic recursive multilevel solver). pARMS isn't exactly for pure ILU decompositions; it works in a multi-level framework not unlike algebraic multigrid. (There's a good explanation in Yousef Saad's book, which is also a good reference for how parallel ILU works.) Nonetheless, it's parallel and you may be able to set the number of levels to 0 to recover the usual ILU factorization.

Finally, Trilinos has a parallel ILU preconditioner. However, it only computes ILU factorizations locally to each processor and uses some overlap to guarantee that the method is scalable. From what I gather, this isn't exactly what you're looking for.

Answer 2

Trilinos provides an incomplete Cholesky preconditioner in two packages, AztecOO and Ifpack. AztecOO provides a pattern-based incomplete factorization using the concept of level-fill. Ifpack provides both level-fill and threshold-based drop tolerance approaches.

These implementations are intended as subdomain solvers as part of an additive Schwarz algebraic domain decomposition preconditioner, or as a local smoother for an algebraic multigrid preconditioner. They can be executed in parallel on a multicore processor, but only if you enable the use of MPI. In this case, each MPI process is assigned a portion of the matrix equations on which an incomplete Cholesky preconditioner will be used. If you compile Trilinos with MPI support disabled, you can still execute the incomplete Cholesky preconditioner, but only on a single core.

Details for each are here:

http://trilinos.org/oldsite/packages/aztecoo/AztecOOUserGuide.pdf (starting about p 16)

http://trilinos.org/oldsite/packages/ifpack/IfpackUserGuide.pdf (starting about p 15)

Using MPI to obtain parallelism on a multicore processor is quite effective, but it does typically require a fairly substantial refactoring of your code if you are running a sequential or OpenMP based code at this time.

Obtaining good threaded parallel performance, e.g., with OpenMP, from incomplete Cholesky on a typical sparse matrix is challenging. The factorization phase can obtain reasonable speedups, say 5x on 8 cores (based on my own experience) for large enough problems. But the solve phase (which you will call for each iteration) is difficult to improve substantially for general sparse matrices. I have seldom seen more than a factor of 2 improvement, no matter how many cores are used.

Complete sparse Cholesky algorithms have a rich graph theory framework that enables organization of the factorization and solve into multifrontal (task) and supernodal (data) parallelism, and can benefit from using optimized dense BLAS. Also, most of the computational time for complete sparse Cholesky is spent in the factorization (again easier to parallelize) and the solve is typically only called once (since the factorization is complete). So threaded complete Cholesky is typically quite effective.

Incomplete Cholesky does not have these same favorable properties, which is why threaded parallelism typically doesn't translate over to the incomplete factorization situation.

If you are already using PETSc, I am fairly certain it has a similar capability.