FGMRES

[Million_P_]

All,

I just started learning to use the mkl and fgmres in general and I'm having troubles implementing/understanding the fgmres. I've looked at the example that came in the library and I've been able to run it successfully but I do have a couple of questions. I was trying to do a simple example just to confirm if I understand what's going on. I've reproduced part of my codes below. Looking at the scheme in the manual, the fgmres does not use the coefficient matrix at all in all the calls except when the rci_request is not equal to zero. In the example provided in the manual, the expected solution was used to initialize the right hand side. so at what point does the code actually use the original matrix? In my excerpt below, the rci_request was zero, when I called the "dfgmres".

character(1)::transa
integer ::n,m,lval,ndiag,openstatus,lvall,ndiagg,i,j
parameter (ndiag=5,lval=9,m=9,n=9)
integer :: idiag(ndiag)
double precision ::val(lval,ndiag),rhs(lval), y(m),computed_solution(n)
!-----------------------------------------------------------
! declarations used in the fgmres
!-----------------------------------------------------------
integer :: itercount
integer :: rci_request, rci_request2
integer :: size
parameter (size=128)
integer :: ipar(size)
double precision :: dpar(size),tmp(n*(2*n+1)+(n*(n+9))/2+1)
!-----------------------------------------------------------
idiag= (/-3,-1,0,1,3/)

!read array from file

open(unit=14, file="filetest",status="old",iostat=openstatus)
if(openstatus/=0) stop "cannot read data from file1"
read(14,*), ((val(lvall,ndiagg),lvall=1,lval),ndiagg=1,ndiag)    !This was successfully read  
 

open(unit=17, file="vector",status="old",iostat=openstatus)
if(openstatus/=0) stop "cannot read data from file2"
read(17, *), (rhs(i),i=1,lval)                                                !successfully read
!------------------------------------------------------------
!Initializing the initial guess
!------------------------------------------------------------
    do i=1,m
        computed_solution(i)=1.0
    end do
!--------------------------------------------------------------
!Initialize the solver
!--------------------------------------------------------------
call dfgmres_init(n,computed_solution,rhs,rci_request,ipar,dpar,tmp)
    !if (rci_request /= 0) goto 999
!-----------------------------------------------------------------
!   Setting the desired parameters
! do the restart after 2 iterations
!   Logical parameters
!   do not do the stopping test for the maximal number of iterations
!   do the preconditioned iterations of fgmres method
!   double precision parameters
!   set the relative tolerance to 1.0d-3 instead of default 1.0d-6
!-------------------------------------------------------------------
    
            ipar(15)=2
            ipar(8)=0
            ipar(11)=1
            dpar(1)=1.0d-3
!-------------------------------------------------------------------           
! Check the correctness and consistency of the newly set parameters
!-------------------------------------------------------------------
        call dfgmres_check(n,computed_solution,rhs,rci_request,ipar,dpar,tmp)
        !if(rci_request /=0) goto 999
!----------------------------------------------------------------------
!compute the solution by rci (p)fgmres solver with preconditioning
!reverse communication starts here
!-----------------------------------------------------------------------
call dfgmres(n,computed_solution, rhs, rci_request, ipar, dpar, tmp)
!-----------------------------------------------------------------------
    ipar(13)=0
call dfgmres_get(n,computed_solution,rhs,rci_request2,ipar,dpar,tmp,itercount)
!----------------------------------------------------------------------------
 

[Chao_Y_Intel]

Hello,

FGMRES is some control functions, it does have the matrix computation.  the actually computation (with matrix data ) is by some BLAS functions. 

For example, in the MKL example, the mkl_dcsrgemv() functions will use the matrix data: 

if (RCI_request == 1){

mkl_dcsrgemv (&cvar, &ivar, A, ia, ja, &tmp[ipar[21] - 1],
            &tmp[ipar[22] - 1]);
..
}

if (RCI_request == 2){
     ….
      dfgmres_get (&ivar, computed_solution, b, &RCI_request, ipar, dpar, tmp,
           &itercount);
      /* Compute the current true residual via MKL (Sparse) BLAS routines */
      mkl_dcsrgemv (&cvar, &ivar, A, ia, ja, b, residual);
  ...
 } 

Thanks,
Chao