*     F12FDF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          LICOMM, NIN, NOUT
      PARAMETER        (LICOMM=134,NIN=5,NOUT=6)
      INTEGER          MAXN, MAXNCV, LDV
      PARAMETER        (MAXN=256,MAXNCV=30,LDV=MAXN)
      INTEGER          LCOMM
      PARAMETER        (LCOMM=3*MAXN+MAXNCV*MAXNCV+8*MAXNCV+60)
      INTEGER          IMON, IPOINT
      PARAMETER        (IMON=0,IPOINT=1)
      DOUBLE PRECISION FOUR, ONE, SIX, TWO, ZERO
      PARAMETER        (FOUR=4.0D+0,ONE=1.0D+0,SIX=6.0D+0,TWO=2.0D+0,
     +                 ZERO=0.0D+0)
*     .. Local Scalars ..
      DOUBLE PRECISION H, R1, R2, SIGMA
      INTEGER          IFAIL, INFO, IREVCM, J, N, NCONV, NCV, NEV,
     +                 NITER, NSHIFT
*     .. Local Arrays ..
      DOUBLE PRECISION AD(MAXN), ADL(MAXN), ADU(MAXN), ADU2(MAXN),
     +                 COMM(LCOMM), D(MAXNCV,2), MX(MAXN), RESID(MAXN),
     +                 V(LDV,MAXNCV), X(MAXN)
      INTEGER          ICOMM(LICOMM), IPIV(MAXN)
*     .. External Functions ..
      DOUBLE PRECISION DNRM2
      EXTERNAL         DNRM2
*     .. External Subroutines ..
      EXTERNAL         DCOPY, DGTTRF, DGTTRS, F12FAF, F12FBF, F12FCF,
     +                 F12FDF, F12FEF, MV
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12FDF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) N, NEV, NCV
      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 F12FAF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*     We are solving a generalized problem
         CALL F12FDF('GENERALIZED',ICOMM,COMM,IFAIL)
*     Indicate that we are using the shift and invert mode.
         CALL F12FDF('SHIFTED INVERSE',ICOMM,COMM,IFAIL)
         IF (IPOINT.EQ.1) THEN
*        Use pointers to Workspace in calculating matrix vector products
*        rather than interfacing through the array X
            CALL F12FDF('POINTERS=YES',ICOMM,COMM,IFAIL)
         END IF
*
         H = ONE/DBLE(N+1)
         R1 = (FOUR/SIX)*H
         R2 = (ONE/SIX)*H
         SIGMA = ZERO
         DO 20 J = 1, N
            AD(J) = TWO/H - SIGMA*R1
            ADL(J) = -ONE/H - SIGMA*R2
   20    CONTINUE
         CALL DCOPY(N,ADL,1,ADU,1)
         CALL DGTTRF(N,ADL,AD,ADU,ADU2,IPIV,INFO)
*
         IREVCM = 0
         IFAIL = -1
   40    CONTINUE
         CALL F12FBF(IREVCM,RESID,V,LDV,X,MX,NSHIFT,COMM,ICOMM,IFAIL)
         IF (IREVCM.NE.5) THEN
            IF (IREVCM.EQ.-1) THEN
*           Perform  y <--- OP*x = inv[A-SIGMA*M]*M*x.
               IF (IPOINT.EQ.0) THEN
                  CALL MV(N,X,MX)
                  CALL DCOPY(N,MX,1,X,1)
                  CALL DGTTRS('N',N,1,ADL,AD,ADU,ADU2,IPIV,X,N,INFO)
               ELSE
                  CALL MV(N,COMM(ICOMM(1)),COMM(ICOMM(2)))
                  CALL DGTTRS('N',N,1,ADL,AD,ADU,ADU2,IPIV,COMM(ICOMM(2)
     +                        ),N,INFO)
               END IF
            ELSE IF (IREVCM.EQ.1) THEN
*           Perform y <-- OP*x = inv[A-sigma*M]*M*x;
*           M*x has been saved in COMM(ICOMM(3)) or MX.
               IF (IPOINT.EQ.0) THEN
                  CALL DCOPY(N,MX,1,X,1)
                  CALL DGTTRS('N',N,1,ADL,AD,ADU,ADU2,IPIV,X,N,INFO)
               ELSE
                  CALL DCOPY(N,COMM(ICOMM(3)),1,COMM(ICOMM(2)),1)
                  CALL DGTTRS('N',N,1,ADL,AD,ADU,ADU2,IPIV,COMM(ICOMM(2)
     +                        ),N,INFO)
               END IF
            ELSE IF (IREVCM.EQ.2) THEN
*           Perform  y <--- M*x.
               IF (IPOINT.EQ.0) THEN
                  CALL MV(N,X,MX)
               ELSE
                  CALL MV(N,COMM(ICOMM(1)),COMM(ICOMM(2)))
               END IF
            ELSE IF (IREVCM.EQ.4 .AND. IMON.NE.0) THEN
*           Output monitoring information
               CALL F12FEF(NITER,NCONV,D,D(1,2),ICOMM,COMM)
               WRITE (6,99998) NITER, NCONV, DNRM2(NEV,D(1,2),1)
            END IF
            GO TO 40
         END IF
         IF (IFAIL.EQ.0) THEN
*        Post-Process using F12FCF to compute eigenvalues/values.
            CALL F12FCF(NCONV,D,V,LDV,SIGMA,RESID,V,LDV,COMM,ICOMM,
     +                  IFAIL)
            WRITE (NOUT,99996) NCONV, SIGMA
            DO 60 J = 1, NCONV
               WRITE (NOUT,99995) J, D(J,1)
   60       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 F12FBF Returned with IFAIL = ',I6)
99996 FORMAT (1X,/' The ',I4,' Ritz values of closest to ',F8.4,' are:',
     +       /)
99995 FORMAT (1X,I8,5X,F12.4)
      END
*
      SUBROUTINE MV(N,V,W)
*     .. Parameters ..
      DOUBLE PRECISION ONE, FOUR, SIX
      PARAMETER     (ONE=1.0D+0,FOUR=4.0D+0,SIX=6.0D+0)
*     .. Scalar Arguments ..
      INTEGER       N
*     .. Array Arguments ..
      DOUBLE PRECISION V(N), W(N)
*     .. Local Scalars ..
      DOUBLE PRECISION H
      INTEGER       J
*     .. External Subroutines ..
      EXTERNAL      DSCAL
*     .. Intrinsic Functions ..
      INTRINSIC     DBLE
*     .. Executable Statements ..
      H = ONE/(DBLE(N+1)*SIX)
      W(1) = FOUR*V(1) + V(2)
      DO 20 J = 2, N - 1
         W(J) = V(J-1) + FOUR*V(J) + V(J+1)
   20 CONTINUE
      J = N
      W(J) = V(J-1) + FOUR*V(J)
      CALL DSCAL(N,H,W,1)
      RETURN
      END
