
!Intel MKL RCI (P)FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method) example
!C  Example program for solving non-symmetric indefinite system of equations
!C  Simplest case: no preconditioning and no user-defined stopping tests

PROGRAM FGMRES_NO_PRECON_F
 include  'mkl_rci.fi'
 include  'mkl_blas.fi'
 INTEGER(KIND =8)  :: N 
 PARAMETER (N=5)
 INTEGER(KIND= 8) :: SIZE
 PARAMETER (SIZE=128)

!C Define arrays for the upper triangle of the coefficient matrix
!C Compressed sparse row storage is used for sparse representation

 INTEGER(KIND = 8), DIMENSION(6) :: IA = (/1,3,6,9,12,14/)
 INTEGER(KIND = 8) , DIMENSION(13) :: JA = (/1,3,1,2,4,2,3,5,3,4,5,4,5/)
REAL*8 ,DIMENSION(13)::A=(/1.0,-1.0,-1.0,1.0,-1.0,1.0,-2.0,1.0,-1.0,2.0,-1.0,-1.0,-3.0/)
! Allocate storage for the ?par parameters and the solution/rhs vectors
INTEGER(KIND =8 ), DIMENSION(128) :: IPAR 
REAL*8 ::  DPAR(SIZE), TMP(N*(2*N+1)+(N*(N+9))/2+1)
REAL*8 , DIMENSION(N):: EXPECTED_SOLUTION = (/-1.0,1.0,0.0,1.0,-1.0/)
!      DATA EXPECTED_SOLUTION /-1.0,1.0,0.0,1.0,-1.0/
REAL*8 , DIMENSION(N) :: RHS
REAL*8 , DIMENSION(N) ::COMPUTED_SOLUTION
!---------------------------------------------------------------------------
! Some additional variables to use with the RCI (P)FGMRES solver
!---------------------------------------------------------------------------
 INTEGER(KIND = 8)  :: ITERCOUNT, EXPECTED_ITERCOUNT
 PARAMETER (EXPECTED_ITERCOUNT=5)
 INTEGER(KIND = 8)  :: RCI_REQUEST, I
 REAL*8 :: DVAR
  
      PRINT *,'--------------------------------------------------'
      PRINT *,'The SIMPLEST example of usage of RCI FGMRES solver'
      PRINT *,'to solve a non-symmetric indefinite non-degenerate'
      PRINT *,'       algebraic system of linear equations'
      PRINT *,'--------------------------------------------------'

! Initialize variables and the right hand side through matrix-vector product

      CALL MKL_DCSRGEMV('N', N, A, IA, JA, EXPECTED_SOLUTION, RHS)

! Initialize the initial guess

      DO I=1,N
         COMPUTED_SOLUTION(I)=1.0
      ENDDO

! Initialize the solver

      CALL DFGMRES_INIT(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP)
      PRINT*, 'RCI_REQUEST =   DFGMRES_INIT     ;  ' , RCI_REQUEST
      IF (RCI_REQUEST.NE.0) GOTO 999

!---------------------------------------------------------------------------
! Set the desired parameters:
! LOGICAL parameters:
! do residual stopping test
! do not request for the user defined stopping test
! do the check of the norm of the next generated vector automatically
! DOUBLE PRECISION parameters
! set the relative tolerance to 1.0D-3 instead of default value 1.0D-6
!---------------------------------------------------------------------------
      IPAR(9)=1
      IPAR(10)=0
      IPAR(12)=1
      DPAR(1)=0.001
!Check the correctness and consistency of the newly set parameters
     CALL DFGMRES_CHECK(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP)
     PRINT* , ' IPAR(9)=1 =' , 1
     PRINT*, 'IPAR(10)=0 =', 0 
     PRINT*,  'IPAR(12)=1 =' ,  IPAR(12)
     PRINT*, ' IPAR(13) =',  IPAR(13)
    
     PRINT*, 'RCI_REQUEST =  DFGMRES_CHECK   ;    ' , RCI_REQUEST
     IF (RCI_REQUEST.NE.0) GOTO 999

!--------------------------------------------
!Print the info about the RCI FGMRES method
!--------------------------------------------
     ! PRINT *, ''
     ! PRINT *,'Some info about the current run of RCI FGMRES method:'
     ! PRINT *, ''
     ! IF (IPAR(8).NE.0) THEN
     !    WRITE(*,'(A,I1,A,A)') 'As IPAR(8)=',IPAR(8),', the automatic', ' test for the maximal number of iterations will be'
     !    PRINT *,'performed'
     ! ELSE
     !  WRITE(*,'(A,I1,A,A)') 'As IPAR(8)=',IPAR(8),', the automatic', ' test for the maximal number of iterations will be'
     !  PRINT *,'skipped'
     ! ENDIF
     ! PRINT *,'+++'
     ! IF (IPAR(9).NE.0) THEN
     !  WRITE(*,'(A,I1,A,A)') 'As IPAR(9)=',IPAR(9),', the automatic',' residual test will be performed'
     ! ELSE
     !  WRITE(*,'(A,I1,A,A)') 'As IPAR(9)=',IPAR(9),', the automatic', ' residual test will be skipped'
     ! ENDIF
     ! PRINT *,'+++'
     ! IF (IPAR(10).NE.0) THEN
     !  WRITE(*,'(A,I1,A,A)') 'As IPAR(10)=',IPAR(10),', the', ' user-defined stopping test will be requested via'
     !  PRINT *,'RCI_REQUEST=2'
     ! ELSE
     !  WRITE(*,'(A,I1,A,A)') 'As IPAR(10)=',IPAR(10),', the', ' user-defined stopping test will not be requested, thus,'
     !  PRINT *,'RCI_REQUEST will not take the value 2'
     ! ENDIF
     ! PRINT *,'+++'
     ! IF (IPAR(11).NE.0) THEN
     !  WRITE(*,'(A,I1,A,A)') 'As IPAR(11)=',IPAR(11),', the',' Preconditioned FGMRES iterations will be performed, thus,'
    !   WRITE(*,'(A,A)') 'the preconditioner action will be requested', ' via RCI_REQUEST=3'
    !  ELSE
    !   WRITE(*,'(A,I1,A,A)') 'As IPAR(11)=',IPAR(11),', the', ' Preconditioned FGMRES iterations will not be performed,'
    !   WRITE(*,'(A)') 'thus, RCI_REQUEST will not take the value 3'
    !  ENDIF
    !  PRINT *,'+++'
    !  IF (IPAR(12).NE.0) THEN
    !   WRITE(*,'(A,I1,A,A)') 'As IPAR(12)=',IPAR(12),', the automatic', ' test for the norm of the next generated vector is'
    !   WRITE(*,'(A,A)') 'not equal to zero up to rounding and', ' computational errors will be performed,'
    !   PRINT *,'thus, RCI_REQUEST will not take the value 4'
    !  ELSE
    !   WRITE(*,'(A,I1,A,A)') 'As IPAR(12)=',IPAR(12),', the automatic',' test for the norm of the next generated vector is'
    !   WRITE(*,'(A,A)') 'not equal to zero up to rounding and',' computational errors will be skipped,'
    !   WRITE(*,'(A,A)') 'thus, the user-defined test will be requested',' via RCI_REQUEST=4'
    !  ENDIF
    !  PRINT *,'+++'
!Compute the solution by RCI (P)FGMRES solver without preconditioning
!Reverse Communication starts here


1     CALL DFGMRES(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP)

!If RCI_REQUEST=0, then the solution was found with the required precision

 IF (RCI_REQUEST.EQ.0) GOTO 3
!--------------------------------------------------------------------------
!If RCI_REQUEST=1, then compute the vector A*TMP(IPAR(22))
!and put the result in vector TMP(IPAR(23))
!--------------------------------------------------------------------------
      IF (RCI_REQUEST.EQ.1) THEN
       CALL MKL_DCSRGEMV('N',N, A, IA, JA, TMP(IPAR(22)), TMP(IPAR(23)))
       PRINT*, 'TMP(IPAR(23)) =' , TMP(IPAR(23))
       PRINT*, 'A =' , A
       PRINT*, 'TMP(IPAR(22)) =' , TMP(IPAR(22))
       PAUSE
       GOTO 1
!--------------------------------------------------------------------------
!If RCI_REQUEST=anything else, then DFGMRES subroutine failed
!to compute the solution vector: COMPUTED_SOLUTION(N)
!--------------------------------------------------------------------------
      ELSE
       GOTO 999
      ENDIF
!--------------------------------------------------------------------------
!Reverse Communication ends here
!Get the current iteration number and the FGMRES solution (DO NOT FORGET to
!call DFGMRES_GET routine as COMPUTED_SOLUTION is still containing
!the initial guess!)
!--------------------------------------------------------------------------
3     CALL DFGMRES_GET(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP, ITERCOUNT)
!--------------------------------------------------------------------------
!Print solution vector: COMPUTED_SOLUTION(N) and
!the number of iterations: ITERCOUNT
!--------------------------------------------------------------------------
      PRINT*, 'RCI_REQUEST =   DFGMRES_GET    ;' , RCI_REQUEST     
      PRINT*, 'COMPUTED_SOLUTION =' , COMPUTED_SOLUTION
      PRINT*, 'ITERCOUNT = ', ITERCOUNT

      PRINT *, ''
      PRINT *,' The system has been solved'
      PRINT *, ''
      PRINT *,' The following solution has been obtained:'
      DO I=1,N
         WRITE(*,'(A18,I1,A2,E10.3)') 'COMPUTED_SOLUTION(',I,')=',COMPUTED_SOLUTION(I)
      ENDDO
      PRINT *, ''
      PRINT *,' The expected solution is:'
      DO I=1,N
         WRITE(*,'(A18,I1,A2,E10.3)') 'EXPECTED_SOLUTION(',I,')=', EXPECTED_SOLUTION(I)
         EXPECTED_SOLUTION(I)=EXPECTED_SOLUTION(I)-COMPUTED_SOLUTION(I)
      ENDDO
      PRINT *, ''
      PRINT *,' Number of iterations: ',ITERCOUNT
      
      DVAR = DNRM2(N,EXPECTED_SOLUTION,1);
      IF (ITERCOUNT.EQ.EXPECTED_ITERCOUNT .AND. DVAR.LE.1.0D-14) THEN
         WRITE( *,'(A,A)') 'This example has successfully PASSED', ' through all steps of computation!'
         STOP 0
      ELSE
         WRITE( *,'(A,A,A,I,A,A,A,E,A)') 'This example may have',  ' FAILED as either the number of iterations differs from the', ' expected number of iterations ',EXPECTED_ITERCOUNT,' or the',' computed solution differs much from the expected solution (', 'Euclidean norm is ',DVAR,'), or both.'
         STOP 1
      ENDIF
      
999   WRITE( *,'(A,A,I)') 'This example FAILED as the solver has', ' returned the ERROR code', RCI_REQUEST
      STOP 1
      
      END
