*     F12FGF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          LICOMM, NIN, NOUT
      PARAMETER        (LICOMM=140,NIN=5,NOUT=6)
      INTEGER          MAXBDW, MAXN, MAXNCV, LDAB, LDV
      PARAMETER        (MAXBDW=50,MAXN=1000,MAXNCV=50,LDAB=MAXBDW,
     +                 LDV=MAXN)
      INTEGER          LCOMM
      PARAMETER        (LCOMM=60)
      DOUBLE PRECISION ONE, ZERO
      PARAMETER        (ONE=1.0D+0,ZERO=0.0D+0)
*     .. Local Scalars ..
      DOUBLE PRECISION H2, SIGMA
      INTEGER          I, IDIAG, IFAIL, IFAIL1, ISUB, ISUP, J, KL, KU,
     +                 LO, N, NCONV, NCV, NEV, NX, OUTCHN
*     .. Local Arrays ..
      DOUBLE PRECISION AB(LDAB,MAXN), AX(MAXN), COMM(LCOMM),
     +                 D(MAXNCV,2), MB(1), RESID(MAXN), V(LDV,MAXNCV)
      INTEGER          ICOMM(LICOMM)
*     .. External Functions ..
      DOUBLE PRECISION DNRM2
      EXTERNAL         DNRM2
*     .. External Subroutines ..
      EXTERNAL         DAXPY, DGBMV, F06QHF, F12FFF, F12FGF, X04ABF,
     +                 X04CAF
*     .. Intrinsic Functions ..
      INTRINSIC        DABS
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12FGF Example Program Results'
      OUTCHN = NOUT
      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
*        Initialize communication arrays.
         CALL F12FFF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
*
*           Construct the matrix A in banded form and store in AB.
*           Zero out AB.
            IFAIL = 0
            CALL F06QHF('G',LDAB,N,ZERO,ZERO,AB,LDAB)
*           KU, KL are number of superdiagonals and subdiagonals within
*           the band of matrices A and M.
            KL = NX
            KU = NX
*           Main diagonal of A.
            H2 = ONE/((NX+1)*(NX+1))
            IDIAG = KL + KU + 1
            DO 20 J = 1, N
               AB(IDIAG,J) = 4.0D+0/H2
   20       CONTINUE
*           First subdiagonal and superdiagonal of A.
            ISUP = KL + KU
            ISUB = KL + KU + 2
            DO 60 I = 1, NX
               LO = (I-1)*NX
               DO 40 J = LO + 1, LO + NX - 1
                  AB(ISUP,J+1) = -ONE/H2
                  AB(ISUB,J) = -ONE/H2
   40          CONTINUE
   60       CONTINUE
*           KL-th subdiagonal and KU-th super-diagonal.
            ISUP = KL + 1
            ISUB = 2*KL + KU + 1
            DO 100 I = 1, NX - 1
               LO = (I-1)*NX
               DO 80 J = LO + 1, LO + NX
                  AB(ISUP,NX+J) = -ONE/H2
                  AB(ISUB,J) = -ONE/H2
   80          CONTINUE
  100       CONTINUE
*
*           Find eigenvalues of largest magnitude and the corresponding
*           eigenvectors.
            IFAIL = -1
            CALL F12FGF(KL,KU,AB,LDAB,MB,1,SIGMA,NCONV,D,V,LDV,RESID,V,
     +                  LDV,COMM,ICOMM,IFAIL)
            IF (IFAIL.EQ.0) THEN
*              Compute the residual norm  ||A*x - lambda*x||.
               DO 120 J = 1, NCONV
                  CALL DGBMV('NoTranspose',N,N,KL,KU,ONE,AB(KL+1,1),
     +                       LDAB,V(1,J),1,ZERO,AX,1)
                  CALL DAXPY(N,-D(J,1),V(1,J),1,AX,1)
                  D(J,2) = DNRM2(N,AX,1)
                  D(J,2) = D(J,2)/DABS(D(J,1))
  120          CONTINUE
               WRITE (NOUT,*)
               CALL X04ABF(1,OUTCHN)
               CALL X04CAF('G','N',NCONV,2,D,MAXNCV,
     +                     ' Ritz values and residuals',IFAIL1)
            END IF
         ELSE
            WRITE (NOUT,99998) IFAIL
         END IF
      END IF
*
99999 FORMAT (1X,A,I5)
99998 FORMAT (1X,' ** F12FFF returned with IFAIL = ',I5)
      END
