*     F12ADF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          IMON, LICOMM, NIN, NOUT
      PARAMETER        (IMON=0,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)
      DOUBLE PRECISION ONE, SIX, TWO
      PARAMETER        (ONE=1.0D+0,SIX=6.0D+0,TWO=2.0D+0)
*     .. Local Scalars ..
      DOUBLE PRECISION H, RHO, S, S1, S2, S3, SIGMAI, SIGMAR
      INTEGER          IFAIL, IFAIL1, INFO, IREVCM, J, N, NCONV, NCV,
     +                 NEV, NITER, NSHIFT, NX
*     .. Local Arrays ..
      DOUBLE PRECISION COMM(LCOMM), D(MAXNCV,3), 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 DNRM2
      EXTERNAL         DNRM2
*     .. External Subroutines ..
      EXTERNAL         DCOPY, DGTTRF, DGTTRS, F12AAF, F12ABF, F12ACF,
     +                 F12ADF, F12AEF, MV
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12ADF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) NX, NEV, NCV, RHO, SIGMAR, SIGMAI
      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 F12AAF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
*           Set the mode.
            CALL F12ADF('SHIFTED REAL',ICOMM,COMM,IFAIL)
*           Set problem type
            IFAIL = 0
            CALL F12ADF('GENERALIZED',ICOMM,COMM,IFAIL)
*           Construct C = A - SIGMA*I, and factor C using DGTTRF/F07CDF.
            H = ONE/DBLE(N+1)
            S = RHO/TWO
            S1 = -ONE/H - S - SIGMAR*H/SIX
            S2 = TWO/H - 4.0D+0*SIGMAR*H/SIX
            S3 = -ONE/H + S - SIGMAR*H/SIX
            DO 20 J = 1, N - 1
               DL(J) = S1
               DD(J) = S2
               DU(J) = S3
   20       CONTINUE
            DD(N) = S2
*
            CALL DGTTRF(N,DL,DD,DU,DU2,IPIV,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) THEN
*                 Perform  x <--- OP*x = inv[A-SIGMA*M]*M*x using
*                 DGGTRS/F07CEF
                  CALL MV(N,X)
                  CALL DGTTRS('N',N,1,DL,DD,DU,DU2,IPIV,X,N,INFO)
               ELSE IF (IREVCM.EQ.1) THEN
*                 Perform  x <--- OP*x = inv[A-SIGMA*M]*M*x.
                  CALL DGTTRS('N',N,1,DL,DD,DU,DU2,IPIV,MX,N,INFO)
                  CALL DCOPY(N,MX,1,X,1)
               ELSE IF (IREVCM.EQ.2) THEN
*                 Perform  y <--- M*x
                  CALL MV(N,X)
               ELSE IF (IREVCM.EQ.4 .AND. IMON.NE.0) THEN
*                 Output monitoring information
                  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
            END IF
         ELSE
            WRITE (NOUT,99997) IFAIL
         END IF
      END IF
*
99999 FORMAT (1X,A,I5)
99998 FORMAT (1X,'Iteration',1X,I3,', No. converged =',1X,I3,', norm o',
     +       'f estimates =',E16.8)
99997 FORMAT (1X,' ** F12AAF returned with IFAIL = ',I5)
99996 FORMAT (1X,/' The ',I4,' generalized Ritz values closest to ',
     +       'unity are:',/)
99995 FORMAT (1X,I8,5X,'( ',F12.4,' , ',F12.4,' )')
      END
*
      SUBROUTINE MV(N,V)
*     Compute the in-place matrix vector multiplication X<---M*X,
*     where M is mass matrix formed by using piecewise linear elements
*     on [0,1].
*     .. Parameters ..
      DOUBLE PRECISION ONE, FOUR, SIX
      PARAMETER        (ONE=1.0D0,FOUR=4.0D+0,SIX=6.0D+0)
*     .. Scalar Arguments ..
      INTEGER          N
*     .. Array Arguments ..
      DOUBLE PRECISION V(N)
*     .. Local Scalars ..
      DOUBLE PRECISION H, VM1, VV
      INTEGER          J
*     .. External Subroutines ..
      EXTERNAL         DSCAL
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      VM1 = V(1)
      V(1) = (FOUR*V(1)+V(2))/SIX
      DO 20 J = 2, N - 1
         VV = V(J)
         V(J) = (VM1+FOUR*VV+V(J+1))/SIX
         VM1 = VV
   20 CONTINUE
      V(N) = (VM1+FOUR*V(N))/SIX
*
      H = ONE/DBLE(N+1)
      CALL DSCAL(N,H,V,1)
      RETURN
      END