*     F12ACF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          LICOMM, NIN, NOUT
      PARAMETER        (LICOMM=140,NIN=5,NOUT=6)
      INTEGER          MAXN, MAXNCV, LDV
      PARAMETER        (MAXN=256,MAXNCV=30,LDV=MAXN)
      INTEGER          LCOMM
      PARAMETER        (LCOMM=3*MAXN+3*MAXNCV*MAXNCV+6*MAXNCV+60)
      INTEGER          IMON
      PARAMETER        (IMON=0)
      DOUBLE PRECISION ONE
      PARAMETER        (ONE=1.0D+0)
*     .. Local Scalars ..
      DOUBLE PRECISION H, RHO, SIGMAI, SIGMAR
      INTEGER          IFAIL, IFAIL1, INFO, IREVCM, J, N, NCONV, NCV,
     +                 NEV, NITER, NSHIFT, NX
*     .. Local Arrays ..
      DOUBLE PRECISION COMM(LCOMM), D(MAXNCV,3), MD(MAXN), ME(MAXN-1),
     +                 MX(MAXN), RESID(MAXN), V(LDV,MAXNCV), X(MAXN)
      INTEGER          ICOMM(LICOMM)
*     .. External Functions ..
      DOUBLE PRECISION DNRM2
      EXTERNAL         DNRM2
*     .. External Subroutines ..
      EXTERNAL         AV, DPTTRF, DPTTRS, F12AAF, F12ABF, F12ACF,
     +                 F12ADF, F12AEF, MV
*     .. Intrinsic Functions ..
*
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12ACF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) NX, NEV, NCV, RHO
      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 F12AAF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*        Set the mode.
         CALL F12ADF('REGULAR INVERSE',ICOMM,COMM,IFAIL)
*        Set problem type.
         CALL F12ADF('GENERALIZED',ICOMM,COMM,IFAIL)
*        Use pointers to Workspace in calculating matrix vector products
*        rather than interfacing through the array X
         CALL F12ADF('POINTERS=YES',ICOMM,COMM,IFAIL)
*
*        Construct M, and factorize using DPTTRF/F07JDF.
         H = ONE/DBLE(N+1)
         DO 20 J = 1, N - 1
            MD(J) = 4.0D+0*H
            ME(J) = H
   20    CONTINUE
         MD(N) = 4.0D+0*H
*
         CALL DPTTRF(N,MD,ME,INFO)
*
         IREVCM = 0
         IFAIL = -1
   40    CONTINUE
         CALL F12ABF(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 using DPTTRS/F07JEF.
               CALL AV(NX,RHO,COMM(ICOMM(1)),COMM(ICOMM(2)))
               CALL DPTTRS(N,1,MD,ME,COMM(ICOMM(2)),N,INFO)
            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 if required.
               CALL F12AEF(NITER,NCONV,D,D(1,2),D(1,3),ICOMM,COMM)
               WRITE (6,99998) NITER, NCONV, DNRM2(NEV,D(1,3),1)
            END IF
            GO TO 40
         END IF
         IF (IFAIL.EQ.0) THEN
*           Post-Process using F12ACF to compute eigenvalues/vectors.
            IFAIL1 = 0
            CALL F12ACF(NCONV,D,D(1,2),V,LDV,SIGMAR,SIGMAI,RESID,V,LDV,
     +                  COMM,ICOMM,IFAIL1)
*           Print computed eigenvalues.
            WRITE (NOUT,99996) NCONV
            DO 60 J = 1, NCONV
               WRITE (NOUT,99995) J, D(J,1), D(J,2)
   60       CONTINUE
         ELSE
            WRITE (NOUT,99998) 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 F12ABF Returned with IFAIL = ',I6)
99996 FORMAT (1X,/' The ',I4,' generalized Ritz values of largest ',
     +       'magnitude are:',/)
99995 FORMAT (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )')
      END
*
      SUBROUTINE AV(NX,RHO,V,W)
*     .. Parameters ..
      DOUBLE PRECISION ONE, TWO
      PARAMETER     (ONE=1.0D+0,TWO=2.0D+0)
*     .. Scalar Arguments ..
      DOUBLE PRECISION RHO
      INTEGER       NX
*     .. Array Arguments ..
      DOUBLE PRECISION V(NX*NX), W(NX*NX)
*     .. Local Scalars ..
      DOUBLE PRECISION DD, DL, DU, H, S
      INTEGER       J, N
*     .. Intrinsic Functions ..
      INTRINSIC     DBLE
*     .. Executable Statements ..
      N = NX*NX
      H = ONE/DBLE(N+1)
      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 ..
      DOUBLE PRECISION ONE, FOUR
      PARAMETER     (ONE=1.0D+0,FOUR=4.0D+0)
*     .. Scalar Arguments ..
      INTEGER       NX
*     .. Array Arguments ..
      DOUBLE PRECISION V(NX*NX), W(NX*NX)
*     .. Local Scalars ..
      DOUBLE PRECISION H
      INTEGER       J, N
*     .. External Subroutines ..
      EXTERNAL      DSCAL
*     .. Intrinsic Functions ..
*
      INTRINSIC     DBLE
*     .. 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/DBLE(N+1)
      CALL DSCAL(N,H,W,1)
      RETURN
      END