*     F12AGF 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, LDMB, LDV
      PARAMETER        (MAXBDW=50,MAXN=1000,MAXNCV=50,LDAB=MAXBDW,
     +                 LDMB=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 H, RHO, SIGMAI, SIGMAR
      INTEGER          I, IDIAG, IFAIL, IFAIL1, ISUB, ISUP, J, KL, KU,
     +                 LO, N, NCONV, NCV, NEV, NX
      LOGICAL          FIRST
*     .. Local Arrays ..
      DOUBLE PRECISION AB(LDAB,MAXN), AX(MAXN), COMM(LCOMM),
     +                 D(MAXNCV,3), MB(LDAB,MAXN), MX(MAXN),
     +                 RESID(MAXN), V(LDV,MAXNCV)
      INTEGER          ICOMM(LICOMM)
*     .. External Functions ..
      DOUBLE PRECISION DNRM2, F06BNF
      EXTERNAL         DNRM2, F06BNF
*     .. External Subroutines ..
      EXTERNAL         DAXPY, DGBMV, F06QHF, F12ADF, F12AFF, F12AGF,
     +                 X04ABF, X04CAF
*     .. Intrinsic Functions ..
*
      INTRINSIC        DABS, DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12AGF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) NX, NEV, NCV, 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 = 0
*     Initialize communication arrays.
         CALL F12AFF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*     Set the mode.
         CALL F12ADF('SHIFTED IMAGINARY',ICOMM,COMM,IFAIL)
*     Set problem type
         CALL F12ADF('GENERALIZED',ICOMM,COMM,IFAIL)
*
*     Construct the matrix A in banded form and store in AB.
*     Zero out AB and MB.
         CALL F06QHF('G',LDAB,N,ZERO,ZERO,AB,LDAB)
         CALL F06QHF('G',LDAB,N,ZERO,ZERO,MB,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.
         IDIAG = KL + KU + 1
         DO 20 J = 1, N
            AB(IDIAG,J) = 4.0D+0
            MB(IDIAG,J) = 4.0D+0
   20    CONTINUE
*     First subdiagonal and superdiagonal of A.
         ISUP = KL + KU
         ISUB = KL + KU + 2
         RHO = 100.0D0
         H = ONE/DBLE(NX+1)
         DO 60 I = 1, NX
            LO = (I-1)*NX
            DO 40 J = LO + 1, LO + NX - 1
               AB(ISUB,J+1) = -ONE + 0.5D0*H*RHO
               AB(ISUP,J) = -ONE - 0.5D0*H*RHO
   40       CONTINUE
   60    CONTINUE
*
         DO 80 J = 1, N - 1
            MB(ISUB,J+1) = ONE
            MB(ISUP,J) = ONE
   80    CONTINUE
*     KL-th subdiagonal and KU-th super-diagonal.
         ISUP = KL + 1
         ISUB = 2*KL + KU + 1
         DO 120 I = 1, NX - 1
            LO = (I-1)*NX
            DO 100 J = LO + 1, LO + NX
               AB(ISUP,NX+J) = -ONE
               AB(ISUB,J) = -ONE
  100       CONTINUE
  120    CONTINUE
*
*     Find eigenvalues closest in value to SIGMA and corresponding
*     eigenvectors.
         IFAIL = -1
         CALL F12AGF(KL,KU,AB,LDAB,MB,LDMB,SIGMAR,SIGMAI,NCONV,D,D(1,2),
     +               V,LDV,RESID,V,LDV,COMM,ICOMM,IFAIL)
         IF (IFAIL.EQ.0) THEN
*        Compute the residual norm  ||A*x - lambda*x||.
            FIRST = .TRUE.
            DO 140 J = 1, NCONV
               IF (D(J,2).EQ.ZERO) THEN
                  CALL DGBMV('NoTranspose',N,N,KL,KU,ONE,AB(KL+1,1),
     +                       LDAB,V(1,J),1,ZERO,AX,1)
                  CALL DGBMV('NoTranspose',N,N,KL,KU,ONE,MB(KL+1,1),
     +                       LDMB,V(1,J),1,ZERO,MX,1)
                  CALL DAXPY(N,-D(J,1),MX,1,AX,1)
                  D(J,3) = DNRM2(N,AX,1)
                  D(J,3) = D(J,3)/DABS(D(J,1))
               ELSE IF (FIRST) THEN
                  CALL DGBMV('NOTRANSPOSE',N,N,KL,KU,ONE,AB(KL+1,1),
     +                       LDAB,V(1,J),1,ZERO,AX,1)
                  CALL DGBMV('NOTRANSPOSE',N,N,KL,KU,ONE,MB(KL+1,1),
     +                       LDAB,V(1,J),1,ZERO,MX,1)
                  CALL DAXPY(N,-D(J,1),MX,1,AX,1)
                  CALL DGBMV('NOTRANSPOSE',N,N,KL,KU,ONE,MB(KL+1,1),
     +                       LDMB,V(1,J+1),1,ZERO,MX,1)
                  CALL DAXPY(N,D(J,2),MX,1,AX,1)
                  D(J,3) = DNRM2(N,AX,1)
                  CALL DGBMV('NOTRANSPOSE',N,N,KL,KU,ONE,AB(KL+1,1),
     +                       LDAB,V(1,J+1),1,ZERO,AX,1)
                  CALL DGBMV('NOTRANSPOSE',N,N,KL,KU,ONE,MB(KL+1,1),
     +                       LDMB,V(1,J+1),1,ZERO,MX,1)
                  CALL DAXPY(N,-D(J,1),MX,1,AX,1)
                  CALL DGBMV('NOTRANSPOSE',N,N,KL,KU,ONE,MB(KL+1,1),
     +                       LDMB,V(1,J),1,ZERO,MX,1)
                  CALL DAXPY(N,-D(J,2),MX,1,AX,1)
                  D(J,3) = F06BNF(D(J,3),DNRM2(N,AX,1))
                  D(J,3) = D(J,3)/F06BNF(D(J,1),D(J,2))
                  D(J+1,3) = D(J,3)
                  FIRST = .FALSE.
               ELSE
                  FIRST = .TRUE.
               END IF
  140       CONTINUE
            WRITE (NOUT,*)
            CALL X04ABF(1,NOUT)
            CALL X04CAF('G','N',NCONV,2,D,MAXNCV,
     +                  ' Ritz values closest to sigma',IFAIL1)
         ELSE
            WRITE (NOUT,99998) IFAIL
         END IF
      END IF
*
99999 FORMAT (1X,A,I5)
99998 FORMAT (1X,' NAG Routine F12AGF Returned with IFAIL = ',I6)
      END