Question about using fgmres

[Zhanghong_T_]

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

[Intel_C_Intel]

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

[Zhanghong_T_]

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

[Zhanghong_T_]

Hi Bruce,

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

Thanks,

Zhanghong Tang

[Intel_C_Intel]

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

[Zhanghong_T_]

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