!   F12ARF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE f12arfe_mod

!      F12ARF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       COMPLEX (KIND=nag_wp), PARAMETER    :: four = (4.0_nag_wp,0.0_nag_wp)
       COMPLEX (KIND=nag_wp), PARAMETER    :: one = (1.0_nag_wp,0.0_nag_wp)
       COMPLEX (KIND=nag_wp), PARAMETER    :: six = (6.0_nag_wp,0.0_nag_wp)
       COMPLEX (KIND=nag_wp), PARAMETER    :: two = (2.0_nag_wp,0.0_nag_wp)
       INTEGER, PARAMETER                  :: imon = 0, licomm = 140,          &
                                              nerr = 6, nin = 5, nout = 6
    CONTAINS
!dec$ attributes offload : mic :: mv
       SUBROUTINE mv(nx,v,w)
!         Compute the out-of--place matrix vector multiplication Y<---M*X,
!         where M is mass matrix formed by using piecewise linear elements
!         on [0,1].

!         .. Use Statements ..
          USE nag_library, ONLY : zscal
!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: nx
!         .. Array Arguments ..
          COMPLEX (KIND=nag_wp), INTENT (IN)  :: v(nx*nx)
          COMPLEX (KIND=nag_wp), INTENT (OUT) :: w(nx*nx)
!         .. Local Scalars ..
          COMPLEX (KIND=nag_wp)               :: h
          INTEGER                             :: j, n
!         .. Intrinsic Functions ..
          INTRINSIC                              cmplx
!         .. Executable Statements ..
          n = nx*nx
          w(1) = (four*v(1)+v(2))/six
          DO j = 2, n - 1
             w(j) = (v(j-1)+four*v(j)+v(j+1))/six
          END DO
          w(n) = (v(n-1)+four*v(n))/six

          h = one/cmplx(n+1,kind=nag_wp)
!         The NAG name equivalent of zscal is f06gdf
          CALL zscal(n,h,w,1)
          RETURN
!dec$ attributes offload : mic :: mv
       END SUBROUTINE mv
    END MODULE f12arfe_mod
    PROGRAM f12arfe

!      F12ARF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : dznrm2, f12anf, f12apf, f12aqf, f12arf, f12asf, &
                               zgttrf, zgttrs
       USE f12arfe_mod, ONLY : four, imon, licomm, mv, nag_wp, nerr, nin,      &
                               nout, one, six, two
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       COMPLEX (KIND=nag_wp)               :: h, rho, s, s1, s2, s3, sigma
       INTEGER                             :: ifail, ifail1, info, irevcm, j,  &
                                              lcomm, ldv, n, nconv, ncv, nev,  &
                                              niter, nshift, nx
!      .. Local Arrays ..
       COMPLEX (KIND=nag_wp), ALLOCATABLE  :: ax(:), comm(:), d(:,:), dd(:),   &
                                              dl(:), du(:), du2(:), mx(:),     &
                                              resid(:), v(:,:), x(:)
       INTEGER                             :: icomm(licomm)
       INTEGER, ALLOCATABLE                :: ipiv(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              cmplx
!      .. Executable Statements ..
       WRITE (nout,*) 'F12ARF Example Program Results'
       WRITE (nout,*)
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) nx, nev, ncv

       n = nx*nx
       lcomm = 3*n + 3*ncv*ncv + 5*ncv + 60
       ldv = n
       ALLOCATE (comm(lcomm),ax(n),d(ncv,2),dd(n),dl(n),du(n),du2(n),mx(n), &
          resid(n),v(ldv,ncv),x(n),ipiv(n))

       ifail = 0
       CALL f12anf(n,nev,ncv,icomm,licomm,comm,lcomm,ifail)

!      Set the mode.
       ifail = 0
       CALL f12arf('SHIFTED INVERSE',icomm,comm,ifail)
!      Set problem type.
       CALL f12arf('GENERALIZED',icomm,comm,ifail)
       sigma = (500.0_nag_wp,0.0_nag_wp)
       rho = (10.0_nag_wp,0.0_nag_wp)
       h = one/cmplx(n+1,kind=nag_wp)
       s = rho/two
       s1 = -one/h - s - sigma*h/six
       s2 = two/h - four*sigma*h/six
       s3 = -one/h + s - sigma*h/six

       dl(1:n-1) = s1
       dd(1:n-1) = s2
       du(1:n-1) = s3
       dd(n) = s2

!      The NAG name equivalent of zgttrf is f07crf
       CALL zgttrf(n,dl,dd,du,du2,ipiv,info)
       IF (info/=0) THEN
          WRITE (nerr,99999) info
          GO TO 20
       END IF

       irevcm = 0
       ifail = -1
REVCM: DO
          CALL f12apf(irevcm,resid,v,ldv,x,mx,nshift,comm,icomm,ifail)
          IF (irevcm==5) THEN
             EXIT REVCM
          ELSE IF (irevcm==-1) THEN
!            Perform  x <--- OP*x = inv[A-SIGMA*M]*M*x
             CALL mv(nx,x,ax)
             x(1:n) = ax(1:n)
!            The NAG name equivalent of zgttrs is f07csf
             CALL zgttrs('N',n,1,dl,dd,du,du2,ipiv,x,n,info)
             IF (info/=0) THEN
                WRITE (nerr,99998) info
                EXIT REVCM
             END IF
          ELSE IF (irevcm==1) THEN
!            Perform  x <--- OP*x = inv[A-SIGMA*M]*M*x,
!            MX stored in COMM from location IPNTR(3)
!            The NAG name equivalent of zgttrs is f07csf
             CALL zgttrs('N',n,1,dl,dd,du,du2,ipiv,mx,n,info)
             x(1:n) = mx(1:n)
             IF (info/=0) THEN
                WRITE (nerr,99998) info
                EXIT REVCM
             END IF
          ELSE IF (irevcm==2) THEN
!            Perform  y <--- M*x
             CALL mv(nx,x,ax)
             x(1:n) = ax(1:n)
          ELSE IF (irevcm==4 .AND. imon/=0) THEN
!            Output monitoring information
             CALL f12asf(niter,nconv,d,d(1,2),icomm,comm)
!            The NAG name equivalent of dznrm2 is f06jjf
             WRITE (6,99997) niter, nconv, dznrm2(nev,d(1,2),1)
          END IF
       END DO REVCM

       IF (ifail==0 .AND. info==0) THEN
!         Post-Process using F12AQF to compute eigenvalues/vectors.
          ifail1 = 0
          CALL f12aqf(nconv,d,v,ldv,sigma,resid,v,ldv,comm,icomm,ifail1)
          WRITE (nout,99996) nconv, sigma
          WRITE (nout,99995) (j,d(j,1),j=1,nconv)
       END IF
20     CONTINUE

99999  FORMAT (1X,'** Error status returned by ZGTTRF, INFO =',I12)
99998  FORMAT (1X,'** Error status returned by ZGTTRS, INFO =',I12)
99997  FORMAT (1X,'Iteration',1X,I3,', No. converged =',1X,I3,', norm o', &
          'f estimates =',E16.8)
99996  FORMAT (1X/' The ',I4,' generalized Ritz values closest to (',F7.3,',', &
          F7.3,') are:'/)
99995  FORMAT (1X,I8,5X,'( ',F10.4,' , ',F10.4,' )')
    END PROGRAM f12arfe
