MKL preconditioned conjugate gradient (PCG)
Preconditioners for conjugate gradient solver of linear system A * x = b in Intel MKL
| |
Re: MKL preconditioned conjugate gradient (PCG)
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
| |
Re: MKL preconditioned conjugate gradient (PCG)
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
| |
Re: MKL preconditioned conjugate gradient (PCG)
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
| |
Re: MKL preconditioned conjugate gradient (PCG)
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
Intel®MKL Manual near CG interface.
With best regards,
Alexander
| |
Re: MKL preconditioned conjugate gradient (PCG)
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
Intel®MKL 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
| |
Re: MKL preconditioned conjugate gradient (PCG)
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 you store your simmetric matrix A in upper or lower case? May be it is
not totally 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
| |
Re: MKL preconditioned conjugate gradient (PCG)
Hi Idefix,
Am I right that you store your simmetric matrix A
in upper or lower case? May be it is not totally 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
| |
Re: MKL preconditioned conjugate gradient (PCG)
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 but no 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
| | |
