Incomplete Cholesky Factorization

[Mikhail_Matrosov]

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)))

[Alexander_K_Intel2]

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

[Mikhail_Matrosov]

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?

[Alexander_K_Intel2]

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

[Mikhail_Matrosov]

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.

[Alexander_K_Intel2]

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

[Mikhail_Matrosov]

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?

[Alexander_K_Intel2 · Accepted solution]

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

[Royi]

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 :-).

[Mikhail_Matrosov]

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

[IDZ_A_Intel]

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