MKL preconditioned conjugate gradient (PCG)

[Idefix]

Preconditioners for conjugate gradient solver of linear system A * x = b in Intel MKL

[Idefix]

Quoting - Idefix

Preconditioners for conjugate gradient solver of linear system A * x = b in Intel MKL

Hi all, let me come to my question now. I am trying to implement a precondtioned conjugate gradient solver for a system A*x=b were A is a symmetric matrix. Without precondtioning I get the correct solution in a transient flow simulation but the solution is too slow. From the MKL Reference manual I infer that preconditioners are supported by the DCG routine but it is also stated that:

'Both ILU0 and ILUT preconditioners can apply to any non-degenerate matrix. They can be used
alone or together with the Intel MKL RCI FGMRES solver (see Sparse Solver Routines). Avoid
using this preconditioners with MKL RCI CG solver because in general, they produce
non-symmetric resulting matrix even if the original matrix is symmetric.'

THis is bad news because other people have used incomplete lower upper decomposition with great success in their PCG codes. So what preconditioners are recommended by the Intel team and where can I find them? Actually, I tried ILUO in my problem despite the above warning but the solution did not converge anymore. Here is a brief extract from my code:

call dcg_init(N,xopt,B,istat,ipar,dpar,tmp)
call dcg_check(N,xopt,B,istat,ipar,dpar,tmp)
call DCSRILU0(N,AP,IA,JA,CP,ipar,dpar,ierr)
201 call dcg(N,xopt,B,istat,ipar,dpar,tmp)
if (istat.eq.0) then
goto 210
elseif (istat.eq.1) then
call MKL_DCSRSYMV('U',N,AP,IA,JA,TMP,TMP(1,2))
goto 201
elseif (istat.eq.3) then
call mkl_dcsrtrsv('L','N','U',N,CP,IA,JA,tmp(1,3),trvec)
call mkl_dcsrtrsv('U','N','N',N,CP,IA,JA,trvec,tmp(1,4))
goto 201
else
write(*,*) 'Error in MKL solver.'
stop
endif
210 call dcg_get(NumNP,xopt,B,istat,ipar,dpar,tmp,it)

If anyone has an idea where my bug is or what preconditioner I could use instead of ILU, I would appreciate it!

Idefix

[Alexander_K_Intel2]

Quoting - Idefix

Hi all, let me come to my question now. I am trying to implement a precondtioned conjugate gradient solver for a system A*x=b were A is a symmetric matrix. Without precondtioning I get the correct solution in a transient flow simulation but the solution is too slow. From the MKL Reference manual I infer that preconditioners are supported by the DCG routine but it is also stated that:

'Both ILU0 and ILUT preconditioners can apply to any non-degenerate matrix. They can be used
alone or together with the Intel MKL RCI FGMRES solver (see Sparse Solver Routines). Avoid
using this preconditioners with MKL RCI CG solver because in general, they produce
non-symmetric resulting matrix even if the original matrix is symmetric.'

THis is bad news because other people have used incomplete lower upper decomposition with great success in their PCG codes. So what preconditioners are recommended by the Intel team and where can I find them? Actually, I tried ILUO in my problem despite the above warning but the solution did not converge anymore. Here is a brief extract from my code:

call dcg_init(N,xopt,B,istat,ipar,dpar,tmp)
call dcg_check(N,xopt,B,istat,ipar,dpar,tmp)
call DCSRILU0(N,AP,IA,JA,CP,ipar,dpar,ierr)
201 call dcg(N,xopt,B,istat,ipar,dpar,tmp)
if (istat.eq.0) then
goto 210
elseif (istat.eq.1) then
call MKL_DCSRSYMV('U',N,AP,IA,JA,TMP,TMP(1,2))
goto 201
elseif (istat.eq.3) then
call mkl_dcsrtrsv('L','N','U',N,CP,IA,JA,tmp(1,3),trvec)
call mkl_dcsrtrsv('U','N','N',N,CP,IA,JA,trvec,tmp(1,4))
goto 201
else
write(*,*) 'Error in MKL solver.'
stop
endif
210 call dcg_get(NumNP,xopt,B,istat,ipar,dpar,tmp,it)

If anyone has an idea where my bug is or what preconditioner I could use instead of ILU, I would appreciate it!

Idefix

Hi Idefix,
Could you please describe some details? First of all what the value of ipar(5) and ipar(11) before calling dcg_check? The value of error with/without preconditioner after computation? Is it matrix A positive define? It can help us to understand your problem with ILU and find the way to solve it.
With best regards,
Alexander

[Idefix]

Quoting - Alexander Kalinkin (Intel)

Hi Idefix,
Could you please describe some details? First of all what the value of ipar(5) and ipar(11) before calling dcg_check? The value of error with/without preconditioner after computation? Is it matrix A positive define? It can help us to understand your problem with ILU and find the way to solve it.
With best regards,
Alexander

Hi Alexander,

thanks for your reply. I computed the eigenvalues of the matrix A in my example by the LAPACK routines DSYTRD and DSTERF. All eigenvalues are positive, hence the matrix is positive definite. The maximum eigenvalue is 46.6 and the minimum eigenvalue is 0.164.

Here are some more details:

IPAR(5)= 150
This maximum number of iterations is sufficient if I run DCG without preconditioner. I have made an additional test by setting it to 1000 if I use the preconditioner but the solution does still not converge. The routine DCG returns values of RCI_Request (istat in my example) of 1, 3, 1, 3 [] and finally -1 (maximum number of iterations exceeded). Without preconditioner, it takes between 30 to 50 iterations until convergence if the relative tolerance DPAR(1) is set to 1.D-6.

IPAR(8)=1
Routine DCG performs stopping test for number of iterations

IPAR(9)=1
Routine DCG performs residual stopping test

IPAR(10)=0
No user-defined stopping test

IPAR(11)=1
Use preconditioner.

DPAR(1)=1.0D-6
Relative tolerance. Increasing this value to 1.D-3 does not help to reach convergence if I use the preconditioner, without preconditioner it lowers the number of iterations required to reach convergence to about 20. All other values of array DPAR are at their default values.

I am looking forward to your reply,

Idefix

[Alexander_K_Intel2]

Quoting - Idefix

Hi Alexander,

thanks for your reply. I computed the eigenvalues of the matrix A in my example by the LAPACK routines DSYTRD and DSTERF. All eigenvalues are positive, hence the matrix is positive definite. The maximum eigenvalue is 46.6 and the minimum eigenvalue is 0.164.

Here are some more details:

IPAR(5)= 150
This maximum number of iterations is sufficient if I run DCG without preconditioner. I have made an additional test by setting it to 1000 if I use the preconditioner but the solution does still not converge. The routine DCG returns values of RCI_Request (istat in my example) of 1, 3, 1, 3 [] and finally -1 (maximum number of iterations exceeded). Without preconditioner, it takes between 30 to 50 iterations until convergence if the relative tolerance DPAR(1) is set to 1.D-6.

IPAR(8)=1
Routine DCG performs stopping test for number of iterations

IPAR(9)=1
Routine DCG performs residual stopping test

IPAR(10)=0
No user-defined stopping test

IPAR(11)=1
Use preconditioner.

DPAR(1)=1.0D-6
Relative tolerance. Increasing this value to 1.D-3 does not help to reach convergence if I use the preconditioner, without preconditioner it lowers the number of iterations required to reach convergence to about 20. All other values of array DPAR are at their default values.

I am looking forward to your reply,

Idefix

Hi Idefix,
All parameters look nice so the situation that you described below could be caused by non positive define of preconditioner. I think that the best way to solve this problem use not CG that could work correct only with positive define matrices but dfgmres routines. Interface of thus functions similar with interface of CG and describe in IntelMKL Manual near CG interface.
With best regards,
Alexander

[Idefix]

Quoting - Alexander Kalinkin (Intel)

Hi Idefix,
All parameters look nice so the situation that you described below could be caused by non positive define of preconditioner. I think that the best way to solve this problem use not CG that could work correct only with positive define matrices but dfgmres routines. Interface of thus functions similar with interface of CG and describe in IntelMKL Manual near CG interface.
With best regards,
Alexander

Hi Alexander,

ok, DFGMRES might probably work but it should be less efficient than PCG for my problem. We have used an unparallelized version of ILU for preconditioning CG for this problem before and it is robust and very quick. I am wondering now about something else. Could it be that DCSRILUO cannot be used at all for symmetric matrices? I mean it should return a nonsymmetric matrix composed of the upper and lower which means that I cannot use the same arrays IA and JA in the call to MKL_DCSRTRSV because the structure of the lower part is not passed to this routine. This structure is basically the same as the upper but maybe MKL_DCSRTRSV cannot extract it because my arrays IA and JA are for a symmetric matrix. Could this be the problem? And if so, could DCSRILUT be a way out because it returns modified arrays ibilut and jbilut?

Idefix

[Alexander_K_Intel2]

Quoting - Idefix

Hi Alexander,

ok, DFGMRES might probably work but it should be less efficient than PCG for my problem. We have used an unparallelized version of ILU for preconditioning CG for this problem before and it is robust and very quick. I am wondering now about something else. Could it be that DCSRILUO cannot be used at all for symmetric matrices? I mean it should return a nonsymmetric matrix composed of the upper and lower which means that I cannot use the same arrays IA and JA in the call to MKL_DCSRTRSV because the structure of the lower part is not passed to this routine. This structure is basically the same as the upper but maybe MKL_DCSRTRSV cannot extract it because my arrays IA and JA are for a symmetric matrix. Could this be the problem? And if so, could DCSRILUT be a way out because it returns modified arrays ibilut and jbilut?

Idefix

Hi Idefix,
Am I right that youstore your simmetric matrix A in upper or lowercase? May be it is nottotally correct for ILU0 preconditioner cause DCSRILU0 calculate preconditioner based on full arrays of IA and JA. If your ia and ja arrays are for upper/lower triangular matrix DCSRILU0 would create preconditioner for upper/lower triangular matrix. Will be such preconditioner sufficiently good for matrix A or not nobody can say.
With best regards,
Alexander

[Idefix]

Quoting - Alexander Kalinkin (Intel)

Hi Idefix,
Am I right that youstore your simmetric matrix A in upper or lowercase? May be it is nottotally correct for ILU0 preconditioner cause DCSRILU0 calculate preconditioner based on full arrays of IA and JA. If your ia and ja arrays are for upper/lower triangular matrix DCSRILU0 would create preconditioner for upper/lower triangular matrix. Will be such preconditioner sufficiently good for matrix A or not nobody can say.
With best regards,
Alexander

Alexander,

sure, I do only store the upper part of the symmetric matrix A, as this is in agreement with the sparse matrix storage format. From your last reply I conclude that use of DCSRILUO is impossible for CG in MKL because symmetric matrices stored in sparse matrix format cannot be handled by DCSRILUO? It looks like I will have to switch to dfgmres (as you suggested earlier) or stick to my old solver package and parallelize it with OpenMP. Is it planned to support ILU as preconditioner for CG in future releases? I think this would be great because it is so often used.

Thanks again for your help.

Idefix

[Alexander_K_Intel2]

Quoting - Idefix

Alexander,

sure, I do only store the upper part of the symmetric matrix A, as this is in agreement with the sparse matrix storage format. From your last reply I conclude that use of DCSRILUO is impossible for CG in MKL because symmetric matrices stored in sparse matrix format cannot be handled by DCSRILUO? It looks like I will have to switch to dfgmres (as you suggested earlier) or stick to my old solver package and parallelize it with OpenMP. Is it planned to support ILU as preconditioner for CG in future releases? I think this would be great because it is so often used.

Thanks again for your help.

Idefix

Idefix
You can use DCSRILU0 for CG in way you suggest butno one can say anything about condition number of preconditioned system and, as result, the number of iteration.If you want to see some new functionality in MKL you can file feature request at https://premier.intel.com
With best regards,
Alexander