subroutine FGMRES(A, IA, JA, RHS, N, nzmax, COMPUTED_SOLUTION)
implicit none
integer N
integer SIZE
parameter (SIZE=128)
integer, intent(in) :: nzmax
!---------------------------------------------------------------------------
! Define arrays for the upper triangle of the coefficient matrix
! Compressed sparse row storage is used for sparse representation
!---------------------------------------------------------------------------
integer IA(N+1)
integer JA(nzmax)
real*8 A(nzmax)
!---------------------------------------------------------------------------
! Allocate storage for the ?par parameters and the solution/rhs/residual vectors
!---------------------------------------------------------------------------
integer IPAR(SIZE)
real*8 DPAR(SIZE), TMP(N*(2*150+1)+(150*(150+9))/2+1)
real*8 RHS(N) !, B(N)
real*8 COMPUTED_SOLUTION(N)
! real*8 RESIDUAL(N)
!---------------------------------------------------------------------------
! Some additional variables to use with the RCI (P)FGMRES solver
!---------------------------------------------------------------------------
integer ITERCOUNT, EXPECTED_ITERCOUNT
parameter (EXPECTED_ITERCOUNT=233)
integer RCI_REQUEST, I
real*8 DVAR
!---------------------------------------------------------------------------
! An external BLAS function is taken from MKL BLAS to use
! with the RCI (P)FGMRES solver
! Save the right-hand side in vector B for future use
!---------------------------------------------------------------------------
! call DCOPY(N, RHS, 1, B, 1)
!---------------------------------------------------------------------------
! Initialize the initial guess
!---------------------------------------------------------------------------
do I=1,N
COMPUTED_SOLUTION(I)=0.
END DO
!---------------------------------------------------------------------------
! Initialize the solver
!---------------------------------------------------------------------------
CALL DFGMRES_INIT(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP)
IF (RCI_REQUEST.NE.0) GOTO 999
!---------------------------------------------------------------------------
! Set the desired parameters:
! do the restart after 2 iterations: IPAR(15)=2
! LOGICAL parameters:
! do not do the stopping test for the maximal number of iterations: IPAR(8)=0
! do the Preconditioned iterations of FGMRES method: IPAR(11)=1
! DOUBLE PRECISION parameters
! set the relative tolerance to 1.0D-3 instead of default value 1.0D-6: DPAR(1)=1.0D-3
!---------------------------------------------------------------------------
IPAR(9)=1
IPAR(10)=0
IPAR(12)=1
DPAR(1)=1.0D-6
!---------------------------------------------------------------------------
! 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.NE.0) GOTO 999
!---------------------------------------------------------------------------
! Compute the solution by RCI (P)FGMRES solver with 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)))
GOTO 1
! END IF
!!---------------------------------------------------------------------------
!! If RCI_request=2, then do the user-defined stopping test
!! The residual stopping test for the computed solution is performed here
!!---------------------------------------------------------------------------
!! NOTE: from this point vector B(N) is no longer containing the right-hand
!! side of the problem! It contains the current FGMRES approximation to the
!! solution. If you need to keep the right-hand side, save it in some other
!! vector before the call to DFGMRES routine. Here we saved it in vector
!! RHS(N). The vector B is used instead of RHS to preserve the original
!! right-hand side of the problem and guarantee the proper restart of FGMRES
!! method. Vector B will be altered when computing the residual stopping
!! criterion!
!!---------------------------------------------------------------------------
! IF (RCI_REQUEST.EQ.2) THEN
!! Request to the DFGMRES_GET routine to put the solution into B(N) via IPAR(13)
! IPAR(13)=1
!! Get the current FGMRES solution in the vector B(N)
! CALL DFGMRES_GET(N, COMPUTED_SOLUTION, B, RCI_REQUEST, IPAR, DPAR, TMP, ITERCOUNT)
!! Compute the current true residual via MKL (Sparse) BLAS routines
! CALL MKL_DCSRGEMV('N', N, A, IA, JA, B, RESIDUAL)
! CALL DAXPY(N, -1.0D0, RHS, 1, RESIDUAL, 1)
! DVAR=DNRM2(N, RESIDUAL, 1)
! IF (DVAR.LT.1.0E-3) THEN
! GOTO 3
! ELSE
! GOTO 1
! END IF
!!---------------------------------------------------------------------------
!! If RCI_REQUEST=anything else, then DFGMRES subroutine failed
!! to compute the solution vector: COMPUTED_SOLUTION(N)
!!---------------------------------------------------------------------------
ELSE
! call data
WRITE(*,*) 'This example FAILED as the solver has returned the ERROR code',RCI_REQUEST
pause
GOTO 999
END IF
!---------------------------------------------------------------------------
! 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)
IF (ITERCOUNT<=EXPECTED_ITERCOUNT) THEN
! WRITE(*,*) 'This example has successfully PASSED through all steps of computation!'
! print*,'number of iteration',ITERCOUNT
RETURN
ELSE
WRITE(*,*) 'This example have FAILED as the number of iterations differs from the expected number of iterations'
pause
STOP 15
END IF
999 WRITE(*,*) 'This example FAILED as the solver has returned the ERROR code',RCI_REQUEST
CALL MKL_FREE_BUFFERS
STOP 16
END SUBROUTINE FGMRES详细解释下这段代码,中文注释
最新发布