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

    MODULE f12aqfe_mod

!      F12AQF 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.0E+0_nag_wp,0.0E+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 :: av
       SUBROUTINE av(nx,v,w)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          COMPLEX (KIND=nag_wp), PARAMETER    :: rho = (10.0_nag_wp,0.0_nag_wp)
!         .. 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)               :: dd, dl, du, h, s
          INTEGER                             :: j, n
!         .. Intrinsic Functions ..
          INTRINSIC                              cmplx
!         .. Executable Statements ..
          n = nx*nx
          h = one/cmplx(n+1,kind=nag_wp)
          s = rho/two
          dd = two/h
          dl = -one/h - s
          du = -one/h + s
          w(1) = dd*v(1) + du*v(2)
          DO j = 2, n - 1
             w(j) = dl*v(j-1) + dd*v(j) + du*v(j+1)
          END DO
          w(n) = dl*v(n-1) + dd*v(n)
          RETURN
!dec$ attributes offload : mic :: av
       END SUBROUTINE av

!dec$ attributes offload : mic :: mv
       SUBROUTINE mv(nx,v,w)

!         .. 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) + one*v(2)
          DO j = 2, n - 1
             w(j) = one*v(j-1) + four*v(j) + one*v(j+1)
          END DO
          w(n) = one*v(n-1) + four*v(n)
          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 f12aqfe_mod
    PROGRAM f12aqfe

!      F12AQF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : dznrm2, f12anf, f12apf, f12aqf, f12arf, f12asf, &
                               zgttrf, zgttrs
       USE f12aqfe_mod, ONLY : av, four, imon, licomm, mv, nag_wp, nerr, nin,  &
                               nout, one
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       COMPLEX (KIND=nag_wp)               :: h, sigma
       INTEGER                             :: ifail, ifail1, info, irevcm, j,  &
                                              lcomm, ldv, n, nconv, ncv, nev,  &
                                              niter, nshift, nx
!      .. Local Arrays ..
       COMPLEX (KIND=nag_wp), ALLOCATABLE  :: comm(:), d(:,:), dd(:), dl(:),   &
                                              du(:), du2(:), mx(:), resid(:),  &
                                              v(:,:), x(:)
       INTEGER                             :: icomm(licomm)
       INTEGER, ALLOCATABLE                :: ipiv(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              cmplx
!      .. Executable Statements ..
       WRITE (nout,*) 'F12AQF 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),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('REGULAR INVERSE',icomm,comm,ifail)
!      Set problem type.
       CALL f12arf('GENERALIZED',icomm,comm,ifail)
!      Use pointers to Workspace rather than interfacing through the array X.
       CALL f12arf('POINTERS=YES',icomm,comm,ifail)
       h = one/cmplx(n+1,kind=nag_wp)

       dl(1:n-1) = h
       dd(1:n-1) = four*h
       du(1:n-1) = h
       dd(n) = four*h

!      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 .OR. irevcm==1) THEN
!            Perform  y <--- OP*x = inv[M]*A*x      |
             CALL av(nx,comm(icomm(1)),comm(icomm(2)))
!            The NAG name equivalent of zgttrs is f07csf
             CALL zgttrs('N',n,1,dl,dd,du,du2,ipiv,comm(icomm(2)),n,info)
             IF (info/=0) THEN
                WRITE (nerr,99998) info
                EXIT REVCM
             END IF
          ELSE IF (irevcm==2) THEN
!            Perform  y <--- M*x
             CALL mv(nx,comm(icomm(1)),comm(icomm(2)))
          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
          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,' Ritz values of largest magnitude are:'/)
99995  FORMAT (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )')
    END PROGRAM f12aqfe
