PROGRAM f12apfe

!      F12APF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : dznrm2, f12anf, f12apf, f12aqf, f12arf, f12asf, &
                               nag_wp, zgttrf, zgttrs
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       COMPLEX (KIND=nag_wp), PARAMETER :: one = (1.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, nerr = 6, nin = 5, nout = 6
!      .. Local Scalars ..
       COMPLEX (KIND=nag_wp)           :: h, h2, rho, s, s1, s2, s3, sigma
       INTEGER                         :: ifail, info, irevcm, j, lcomm, ldv,  &
                                          licomm, 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, ALLOCATABLE            :: icomm(:), ipiv(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          cmplx, int, max
!      .. Executable Statements ..
       WRITE (nout,*) 'F12APF Example Program Results'
       WRITE (nout,*)

!      Skip heading in data file
       READ (nin,*)

       READ (nin,*) nx, nev, ncv
       n = nx*nx

!      Initialize communication arrays.
!      Query the required sizes of the communication arrays.

       licomm = -1
       lcomm = -1
       ALLOCATE (icomm(max(1,licomm)),comm(max(1,lcomm)))

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

       licomm = icomm(1)
       lcomm = int(comm(1))
       DEALLOCATE (icomm,comm)
       ALLOCATE (icomm(max(1,licomm)),comm(max(1,lcomm)))

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

!      Set the mode.

       ifail = 0
       CALL f12arf('SHIFTED INVERSE',icomm,comm,ifail)

       sigma = cmplx(0,kind=nag_wp)
       rho = (10.0_nag_wp,0.0_nag_wp)
       h = one/cmplx(n+1,kind=nag_wp)
       h2 = h*h
       s = rho/two

       s1 = -one/h2 - s/h
       s2 = two/h2 - sigma
       s3 = -one/h2 + s/h

       ALLOCATE (dl(n-1),dd(n),du(n-1),du2(n-2),ipiv(n))
       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 40
       END IF

       ldv = n
       ALLOCATE (resid(n),v(ldv,ncv),x(n),mx(n),d(ncv,2))

       irevcm = 0
       ifail = -1

20     CONTINUE

       CALL f12apf(irevcm,resid,v,ldv,x,mx,nshift,comm,icomm,ifail)

       IF (irevcm/=5) THEN

          IF (irevcm==-1 .OR. irevcm==1) THEN

!            Perform  x <--- OP*x = inv[A-SIGMA*I]*x

!            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
                GO TO 40
             END IF

          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

          GO TO 20

       END IF

       IF (ifail==0) THEN

!         Post-Process using F12AQF to compute eigenvalues/vectors.

          ifail = 0
          CALL f12aqf(nconv,d,v,ldv,sigma,resid,v,ldv,comm,icomm,ifail)

          WRITE (nout,99996) nconv

          DO j = 1, nconv
             WRITE (nout,99995) j, d(j,1)
          END DO

       END IF

40     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 smallest magnitude are:'/)
99995  FORMAT (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )')
    END PROGRAM f12apfe