topic Incomplete Cholesky Factorization in Intel® oneAPI Math Kernel Library

Incomplete Cholesky Factorization

Mikhail_Matrosov — Wed, 31 Aug 2011 10:21:59 GMT

Hello,

I'm using MKL RCI CG solver to solve large sparse SLE with symmetric and positive definite matrix. I would like to use an incomplete Cholesky factorization as a preconditioner. How can I compute it with MKL?

There are routines for generating ILU0 and ILUT preconditioners described in "Preconditioners based on Incomplete LU Factorization Technique" section. However, it said that:

Avoid using these preconditioners with MKL RCI CG solver because in general, they produce non-symmetric resulting matrix even if the original matrix is symmetric. Usually, an inverse of the preconditioner is required in this case. To do this the Intel MKL triangular solver routine mkl_dcsrtrsv must be applied twice: for the low triangular part of the preconditioner, and then for its upper triangular part.

So does these routines can be used to obtain incomplete Cholesky factorization? I.e. to compute L and U such that simply L = U^{T}? And what the phrase "inverse of the preconditioner is required in this case ... solver routine mkl_dcsrtrsv must be applied twice" does exactly means? There is an example of using this preconditioner on general matrix for MKL RCI FGMRES solver and it already utilizes two calls:

call mkl_dcsrtrsv('L','N','U', n, bilut, ibilut, jbilut, tmp(ipar(22)),trvec)

call mkl_dcsrtrsv('U','N','N', n, bilut, ibilut, jbilut, trvec, tmp(ipar(23)))

Incomplete Cholesky Factorization

Alexander_K_Intel2 — Wed, 31 Aug 2011 11:12:48 GMT

Hi Mikhail,

RCI solver compute lower triangular matrix L and upper triangular matrix U that has sparse format and their multiplication is near initial matrix A. But in general case matrix L in not equal U^T. To apply preconditioner on your iterative solver you need to calulate z =(LU)^-1*r, where r is vector.

z=((LU)^-1)*r could be realized in two step:

y = (L^-1)*r;

z =(U^-1)*y using 2 call of mkl_dcsrtrsv.

With best regards,

Alexander Kalinkin

Incomplete Cholesky Factorization

Mikhail_Matrosov — Wed, 31 Aug 2011 12:19:06 GMT

Ok, several more questions:

Do I need to provide my symmetric matrix in full (non-symmetric) format for dcsrilu0 routine? How do the L and U matrices stored in the only output parameter bilu0 of this routine?

Incomplete Cholesky Factorization

Alexander_K_Intel2 — Wed, 31 Aug 2011 12:39:40 GMT

Hi,

You need to provide full matrix format if you want to obtain incomplete Cholesky decomposition for this matrix. If you provide only lower or upper part of initial matrix as input parameter for dcsrilu0 routine you'd gotten incomplete Cholesky decomposition only for this part of matrix. About second question: L and U matrices stored in one matrix described by arrays ia, ja andbilu0.

With best regards,

Alexander Kalinkin

Incomplete Cholesky Factorization

Mikhail_Matrosov — Wed, 31 Aug 2011 14:34:19 GMT

I did all as you said, but it doesn't seem to work. At least, convergance rate is much worse than for simple Jacobi preconditioner, so I assume something went wrong.

Applying (LU)^-1 as a preconditioner is just what is done for RCI FGMRES for non-symmetric matrices. It is also done in MKL example dcsrilu0_exampl1.c for the same purpose. So, are you really sure, I should use the same scenario when working with symmetric matrix and PCI CG?.. I'm not very expirienced in numerical methods, so maybe I'm missing something obvious.

Incomplete Cholesky Factorization

Alexander_K_Intel2 — Wed, 31 Aug 2011 15:03:37 GMT

Mikhail,

Incomplete Cholesky preconditioner could be better than Jacobi in condition number sense in several cases but could be worse than it in other. There is not any theoretical results on quality of such preconditioner in general case. Moreover, incomplete Cholesky preconditioner is unsymmetrical so you can't use it with CG.

With best regards,

Alexander Kalinkin

Incomplete Cholesky Factorization

Mikhail_Matrosov — Wed, 31 Aug 2011 15:27:24 GMT

incomplete Cholesky preconditioner is unsymmetrical so you can't use it with CG

Thats's not true. Incomplete Cholesky factorization is given by A = L * L^T, so it is symmetrical by design, in distinction from incomplete LU factorization. Yes, I cannot use LU factorization, that's from where my questions did arise.

So, again: how can I compute incomplete Cholesky factorization with MKL?

Incomplete Cholesky Factorization

Alexander_K_Intel2 — Wed, 31 Aug 2011 18:46:53 GMT

Ok, I've undersood you. Current version of MKL could calculate only incomplete LU decomposition but not incomplete L*L^T decomposition.

With best regards,

Alexander Kalinkin

Incomplete Cholesky Factorization

Mikhail_Matrosov — Fri, 02 Sep 2011 16:25:13 GMT

That's sad. Is there a chance it will be introduced in the upcoming releases?

Incomplete Cholesky Factorization

IDZ_A_Intel — Sat, 03 Sep 2011 03:28:03 GMT

Hi Mikhail,

If you provide data about your company or university I'd include into Feature Request on Incomplete Cholesky Factorization. You could do it using private answer.

With best regards,

Alexander Kalinkin

Re: Incomplete Cholesky Factorization

Royi — Fri, 10 Jul 2020 07:29:03 GMT

Could you please fix that?

The Incomplete Cholesky is very important.

There are modern algorithms to parallelize it and optimize it. We expect the fastest factorization from the Intel guys :-).