topic Re: Question about using fgmres in Intel® oneAPI Math Kernel Library

Question about using fgmres

Zhanghong_T_ — Mon, 11 Jun 2007 09:30:53 GMT

Hi all,

I am trying to use the fgmres function provided by Intel MKL to solve A*x=b. Since the function provided external matrix-vector operation interface and the input matrix A is not inputted into the functions related to fgmres,I can use it to solve a symmetric matrix and then theinput matrix is lower part.I have the following code:

! FGMRES.f90 
!
! FUNCTIONS:
! FGMRES - Entry point of console application.
!

!****************************************************************************
!
! PROGRAM: FGMRES
!
! PURPOSE: Entry point for the console application.
!
!****************************************************************************

 program FGMRES

 implicit none

 ! Variables
 integer::i,j,n,nr,nz 
 integer,allocatable::irow(:),jcol(:)
 real*8,allocatable::x(:,:),b(:,:),s(:)
 nr=1
 open(1,file='s.txt')
 read(1,*)nz
 allocate(irow(nz),jcol(nz),s(nz))
 do i=1,nz
 read(1,*)irow(i),jcol(i),s(i)
 enddo
 close(1)
 open(1,file='r.txt')
 read(1,*)n
 allocate(x(n,nr),b(n,nr))
 do i=1,n
 read(1,*)b(i,1)
 x(i,1)=0.d0
 enddo
 close(1)
 call gmressolverd(n,nz,nr,irow,jcol,s,x,b)
 ! Body of FGMRES

 end program FGMRES
 
 subroutine gmressolverd(n,nz,nr,irow,jcol,s,x,b)
 implicit none
 integer::ipar(128)
 real*8::dpar(128)
 integer::i,j,n,nr,nz,irow(1),jcol(1),RCI_request,ntemp,iter
 real*8::s(1),x(n,1),b(n,1)
 real*8,allocatable::tmp(:)
 character(6)::mattype
 j=0
 do i=1,nz
 if(irow(i)>jcol(i))then
 if(j==0.or.j==2)j=j+1
 elseif(irow(i) if(j==0.or.j==1)j=j+2
 endif
 if(j==3)exit
 enddo
 if(j==1)then
 mattype='SUNF00'
 elseif(j==2)then
 mattype='SLNF00'
 elseif(j==3)then
 mattype='G00F00'
 endif
 ntemp=150 ! max iterative times
 ntemp=((2*ntemp+1)*n+ntemp*ntemp+9)/2+1
 ntemp=ntemp+ntemp
 allocate(tmp(ntemp))
 call dfgmres_init(n, x, b, RCI_request, ipar, dpar, tmp)
 if(RCI_request/=0)goto 999
 ipar(9)=1
 ipar(10)=0
 ipar(12)=1
 dpar(1)=1.d-20
&nbs
p; call dfgmres_check(n, x(:,1), b(:,1), RCI_request, ipar, dpar, tmp)
 if(RCI_request/=0)goto 999
 j=0
 do
 j=j+1
 write(*,*)'iter=',j
 call dfgmres(n,x(:,1),b(:,1),RCI_request,ipar,dpar,tmp)
 if(RCI_request==0)then
 exit
 elseif(RCI_request==1)then
 call mkl_dcoomv('n',n,n,1.d0,mattype,s,irow,jcol,nz,tmp(ipar(22)),0.d0,tmp(ipar(23)))
 cycle
 else
 goto 999
 endif
 enddo
 call dfgmres_get(n,x(:,1),b(:,1),RCI_request,ipar,dpar,tmp,iter)
 return
999 & 
 write(*,*)'the solver has returned the error code:',RCI_request
 call dfgmres_get(n,x,b,RCI_request,ipar,dpar,tmp,iter)
 end subroutine

However, the result is not correct. Could anyone point out the errors of my code?

Thanks,

Zhanghong Tang

Re: Question about using fgmres

Intel_C_Intel — Mon, 11 Jun 2007 20:54:19 GMT

Zhanghong,

Did you look at the FGMRES examples that comewith MKL? They may be of help to you. At least that is the intention for developing them and providing them with the library.

Bruce

Re: Question about using fgmres

Zhanghong_T_ — Tue, 12 Jun 2007 00:37:30 GMT

Hi Bruce,

Thank you very much for your kindly reply. That code is just from the example, except that I change the matrix-vector product function from 'MKL_DCSRGEMV' to 'mkl_dcoomv', since the input matrix is coordinate format. Could you please help me to locate the problem? I have already attached an input matrix. The 'fgmres' function is a black box for me so I have no idea to check it.

Thanks,

Zhanghong Tang

Re: Question about using fgmres

Zhanghong_T_ — Tue, 12 Jun 2007 01:17:00 GMT

Hi Bruce,

Thank you very much! I have found the problem and the code can run well now.

Thanks,

Zhanghong Tang

Re: Question about using fgmres

Intel_C_Intel — Mon, 18 Jun 2007 15:25:00 GMT

Glad to hear that things are running correctly for you now. Wish I could say that I had been a major help here, but in honesty I can't .

Bruce

Re: Question about using fgmres

Zhanghong_T_ — Mon, 25 Jun 2007 13:05:48 GMT

Hi Bruce and all,

Thanks for your kindly reply!

Another question: It seems that the dynamical allocated array'tmp' is associated with maximal iteration times given by ipar(15) (default value is 150). However, for very large martrix size, the solution can't be converged even witha preconditioner after maximal iteration times. Furthermore, we wish to decrease the size of allocated array 'tmp' to save memory (so the value of ipar(15) will be modified). Basing on these considers, we need to 'restart' the GMRES, which is the great feature of FGMRES.

Could any one tell me how to 'restart' the GMRES based on the example code?

Thanks,
Zhanghong Tang