*     F12ANF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          LICOMM, NIN, NOUT, IMON, IPOINT
      PARAMETER        (LICOMM=140,NIN=5,NOUT=6,IMON=0,IPOINT=0)
      INTEGER          MAXN, MAXNCV, LDV
      PARAMETER        (MAXN=256,MAXNCV=30,LDV=MAXN)
      INTEGER          LCOMM
      PARAMETER        (LCOMM=3*MAXN+3*MAXNCV*MAXNCV+5*MAXNCV+60)
*     .. Local Scalars ..
      COMPLEX *16      SIGMA
      INTEGER          I, IFAIL, IFAIL1, IREVCM, N, NCONV, NCV, NEV,
     +                 NITER, NSHIFT, NX
*     .. Local Arrays ..
      COMPLEX *16      AX(MAXN), COMM(LCOMM), D(MAXNCV,2), MX(MAXN),
     +                 RESID(MAXN), V(LDV,MAXNCV), X(MAXN)
      INTEGER          ICOMM(LICOMM)
*     .. External Functions ..
      DOUBLE PRECISION DZNRM2
      EXTERNAL         DZNRM2
*     .. External Subroutines ..
      EXTERNAL         AV, F12ANF, F12APF, F12AQF, F12ARF, F12ASF, ZCOPY
*     .. Intrinsic Functions ..
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12ANF 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 = 0
         CALL F12ANF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
         IF (IPOINT.NE.0) THEN
*        Use pointers to Workspace in calculating matrix vector products
*        rather than interfacing through the array X
            CALL F12ARF('POINTERS=YES',ICOMM,COMM,IFAIL)
         END IF
         IREVCM = 0
         IFAIL = -1
   20    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 matrix vector multiplication y <--- Op*x
               IF (IPOINT.EQ.0) THEN
                  CALL AV(NX,X,AX)
                  CALL ZCOPY(N,AX,1,X,1)
               ELSE
                  CALL AV(NX,COMM(ICOMM(1)),COMM(ICOMM(2)))
               END IF
            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,99998) NITER, NCONV, DZNRM2(NEV,D(1,2),1)
            END IF
            GO TO 20
         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,99996) NCONV
            DO 40 I = 1, NCONV
               WRITE (NOUT,99995) I, D(I,1)
   40       CONTINUE
         ELSE
            WRITE (NOUT,99997) IFAIL
         END IF
      END IF
      STOP
*
99999 FORMAT (1X,A,I5)
99998 FORMAT (1X,'Iteration',1X,I3,', No. converged =',1X,I3,', norm o',
     +       'f estimates =',E16.8)
99997 FORMAT (1X,' NAG Routine F12APF Returned with IFAIL = ',I6)
99996 FORMAT (1X,/' The ',I4,' Ritz values of largest magnitude are:',/)
99995 FORMAT (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )')
      END
*
      SUBROUTINE AV(NX,V,W)
*     .. Scalar Arguments ..
      INTEGER       NX
*     .. Array Arguments ..
      COMPLEX *16   V(NX*NX), W(NX*NX)
*     .. Local Scalars ..
      COMPLEX *16   H2
      INTEGER       J, LO
*     .. External Subroutines ..
      EXTERNAL      TV, ZAXPY
*     .. Intrinsic Functions ..
      INTRINSIC     DCMPLX
*     .. Executable Statements ..
      H2 = DCMPLX(-(NX+1)*(NX+1))
*
      CALL TV(NX,V(1),W(1))
      CALL ZAXPY(NX,H2,V(NX+1),1,W(1),1)
*
      DO 20 J = 2, NX - 1
         LO = (J-1)*NX
         CALL TV(NX,V(LO+1),W(LO+1))
         CALL ZAXPY(NX,H2,V(LO-NX+1),1,W(LO+1),1)
         CALL ZAXPY(NX,H2,V(LO+NX+1),1,W(LO+1),1)
   20 CONTINUE
*
      LO = (NX-1)*NX
      CALL TV(NX,V(LO+1),W(LO+1))
      CALL ZAXPY(NX,H2,V(LO-NX+1),1,W(LO+1),1)
*
      RETURN
      END
*
      SUBROUTINE TV(NX,X,Y)
*     Compute the matrix vector multiplication y<---T*x where T is a nx
*     by nx tridiagonal matrix.
*     .. Parameters ..
      COMPLEX *16   RHO
      PARAMETER     (RHO=(1.0D+2,0.0D+0))
*     .. Scalar Arguments ..
      INTEGER       NX
*     .. Array Arguments ..
      COMPLEX *16   X(NX), Y(NX)
*     .. Local Scalars ..
      COMPLEX *16   DD, DL, DU, H, H2
      INTEGER       J
*     .. Intrinsic Functions ..
*
      INTRINSIC     DCMPLX
*     .. Executable Statements ..
      H = DCMPLX(NX+1)
      H2 = H*H
      DD = (4.0D+0,0.0D+0)*H2
      DL = -H2 - (5.0D-1,0.0D+0)*RHO*H
      DU = -H2 + (5.0D-1,0.0D+0)*RHO*H
*
      Y(1) = DD*X(1) + DU*X(2)
      DO 20 J = 2, NX - 1
         Y(J) = DL*X(J-1) + DD*X(J) + DU*X(J+1)
   20 CONTINUE
      Y(NX) = DL*X(NX-1) + DD*X(NX)
      RETURN
      END
