*     F12AQF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          IMON, LICOMM, NERR, NIN, NOUT
      PARAMETER        (IMON=0,LICOMM=140,NERR=6,NIN=5,NOUT=6)
      INTEGER          MAXN, MAXNCV, LDV
      PARAMETER        (MAXN=256,MAXNCV=30,LDV=MAXN)
      INTEGER          LCOMM
      PARAMETER        (LCOMM=3*MAXN+3*MAXNCV*MAXNCV+5*MAXNCV+60)
      COMPLEX *16      ONE
      PARAMETER        (ONE=(1.0D+0,0.0D+0))
*     .. Local Scalars ..
      COMPLEX *16      H, SIGMA
      INTEGER          IFAIL, IFAIL1, INFO, IREVCM, J, N, NCONV, NCV,
     +                 NEV, NITER, NSHIFT, NX
*     .. Local Arrays ..
      COMPLEX *16      COMM(LCOMM), D(MAXNCV,2), DD(MAXN), DL(MAXN),
     +                 DU(MAXN), DU2(MAXN), MX(MAXN), RESID(MAXN),
     +                 V(LDV,MAXNCV), X(MAXN)
      INTEGER          ICOMM(LICOMM), IPIV(MAXN)
*     .. External Functions ..
      DOUBLE PRECISION DZNRM2
      EXTERNAL         DZNRM2
*     .. External Subroutines ..
      EXTERNAL         AV, F12ANF, F12APF, F12AQF, F12ARF, F12ASF, MV,
     +                 ZGTTRF, ZGTTRS
*     .. 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
      IF (N.LT.1 .OR. N.GT.MAXN) THEN
         WRITE (NOUT,99999) 'N is out of range: N = ', N
      ELSE IF (NCV.GT.MAXNCV) THEN
         WRITE (NOUT,99999) 'NCV is out of range: NCV = ', NCV
      ELSE
         IFAIL = 1
         CALL F12ANF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
*           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=KIND(H))
*
            DO 20 J = 1, N - 1
               DL(J) = H
               DD(J) = (4.0D+0,0.0D+0)*H
               DU(J) = H
   20       CONTINUE
            DD(N) = (4.0D+0,0.0D+0)*H
*
            CALL ZGTTRF(N,DL,DD,DU,DU2,IPIV,INFO)
            IF (INFO.NE.0) THEN
               WRITE (NERR,99998) INFO
               GO TO 80
            END IF
*
            IREVCM = 0
            IFAIL = -1
   40       CONTINUE
            CALL F12APF(IREVCM,RESID,V,LDV,X,MX,NSHIFT,COMM,ICOMM,IFAIL)
            IF (IREVCM.NE.5) THEN
               IF (IREVCM.EQ.-1 .OR. IREVCM.EQ.1) THEN
*                 Perform  y <--- OP*x = inv[M]*A*x      |
                  CALL AV(NX,COMM(ICOMM(1)),COMM(ICOMM(2)))
                  CALL ZGTTRS('N',N,1,DL,DD,DU,DU2,IPIV,COMM(ICOMM(2)),
     +                        N,INFO)
                  IF (INFO.NE.0) THEN
                     WRITE (NERR,99997) INFO
                     GO TO 80
                  END IF
               ELSE IF (IREVCM.EQ.2) THEN
*                 Perform  y <--- M*x
                  CALL MV(NX,COMM(ICOMM(1)),COMM(ICOMM(2)))
               ELSE IF (IREVCM.EQ.4 .AND. IMON.NE.0) THEN
*                 Output monitoring information
                  CALL F12ASF(NITER,NCONV,D,D(1,2),ICOMM,COMM)
                  WRITE (6,99996) NITER, NCONV, DZNRM2(NEV,D(1,2),1)
               END IF
               GO TO 40
            END IF
            IF (IFAIL.EQ.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,99994) NCONV
               DO 60 J = 1, NCONV
                  WRITE (NOUT,99993) J, D(J,1)
   60          CONTINUE
            END IF
         ELSE
            WRITE (NOUT,99995) IFAIL
         END IF
   80    CONTINUE
      END IF
*
99999 FORMAT (1X,A,I5)
99998 FORMAT (1X,'** Error status returned by ZGTTRF, INFO =',I12)
99997 FORMAT (1X,'** Error status returned by ZGTTRS, INFO =',I12)
99996 FORMAT (1X,'Iteration',1X,I3,', No. converged =',1X,I3,', norm o',
     +       'f estimates =',E16.8)
99995 FORMAT (1X,' ** F12ANF returned with IFAIL = ',I5)
99994 FORMAT (1X,/' The ',I4,' Ritz values of largest magnitude are:',/)
99993 FORMAT (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )')
      END
*
      SUBROUTINE AV(NX,V,W)
*     .. Parameters ..
      COMPLEX *16      ONE, TWO, RHO
      PARAMETER        (ONE=(1.0D+0,0.0D+0),TWO=(2.0D+0,0.0D+0),
     +                 RHO=(1.0D+1,0.0D+0))
*     .. Scalar Arguments ..
      INTEGER          NX
*     .. Array Arguments ..
      COMPLEX *16      V(NX*NX), W(NX*NX)
*     .. Local Scalars ..
      COMPLEX *16      DD, DL, DU, H, S
      INTEGER          J, N
*     .. Intrinsic Functions ..
      INTRINSIC        CMPLX
*     .. Executable Statements ..
      N = NX*NX
      H = ONE/CMPLX(N+1,KIND=KIND(DD))
      S = RHO/TWO
      DD = TWO/H
      DL = -ONE/H - S
      DU = -ONE/H + S
      W(1) = DD*V(1) + DU*V(2)
      DO 20 J = 2, N - 1
         W(J) = DL*V(J-1) + DD*V(J) + DU*V(J+1)
   20 CONTINUE
      W(N) = DL*V(N-1) + DD*V(N)
      RETURN
      END
*
      SUBROUTINE MV(NX,V,W)
*     .. Parameters ..
      COMPLEX *16      ONE, FOUR
      PARAMETER        (ONE=(1.0D+0,0.0D+0),FOUR=(4.0D+0,00D+0))
*     .. Scalar Arguments ..
      INTEGER          NX
*     .. Array Arguments ..
      COMPLEX *16      V(NX*NX), W(NX*NX)
*     .. Local Scalars ..
      COMPLEX *16      H
      INTEGER          J, N
*     .. External Subroutines ..
      EXTERNAL         ZSCAL
*     .. Intrinsic Functions ..
      INTRINSIC        CMPLX
*     .. Executable Statements ..
      N = NX*NX
      W(1) = FOUR*V(1) + ONE*V(2)
      DO 20 J = 2, N - 1
         W(J) = ONE*V(J-1) + FOUR*V(J) + ONE*V(J+1)
   20 CONTINUE
      W(N) = ONE*V(N-1) + FOUR*V(N)
      H = ONE/CMPLX(N+1,KIND=KIND(H))
      CALL ZSCAL(N,H,W,1)
      RETURN
      END