*     F12AEF 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)
      DOUBLE PRECISION ZERO
      PARAMETER        (ZERO=0.0D+0)
*     .. Local Scalars ..
      COMPLEX *16      C1, C2, C3
      DOUBLE PRECISION DENI, DENR, NUMI, NUMR, SIGMAI, SIGMAR
      INTEGER          IFAIL, IFAIL1, INFO, IREVCM, J, N, NCONV, NCV,
     +                 NEV, NITER, NSHIFT
      LOGICAL          FIRST
*     .. Local Arrays ..
      COMPLEX *16      CDD(MAXN), CDL(MAXN), CDU(MAXN), CDU2(MAXN),
     +                 CTEMP(MAXN)
      DOUBLE PRECISION AX(MAXN), COMM(LCOMM), D(MAXNCV,3), MX(MAXN),
     +                 RESID(MAXN), V(LDV,MAXNCV), X(MAXN)
      INTEGER          ICOMM(LICOMM), IPIV(MAXN)
*     .. External Functions ..
      DOUBLE PRECISION DDOT, DNRM2, F06BNF
      EXTERNAL         DDOT, DNRM2, F06BNF
*     .. External Subroutines ..
      EXTERNAL         AV, DCOPY, F12AAF, F12ABF, F12ACF, F12ADF,
     +                 F12AEF, MV, ZGTTRF, ZGTTRS
*     .. Intrinsic Functions ..
*
      INTRINSIC        DBLE, DCMPLX
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12AEF Example Program Results'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) N, NEV, NCV, SIGMAR, SIGMAI
      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('SHIFTED REAL',ICOMM,COMM,IFAIL)
*        Set problem type
         CALL F12ADF('GENERALIZED',ICOMM,COMM,IFAIL)
*        Solve A*x = lambda*B*x in shift-invert mode.
*        The shift, sigma, is a complex number (sigmar, sigmai).
*        OP = Real_Part{inv[A-(SIGMAR,SIGMAI)*M]*M and  B = M.
         C1 = DCMPLX(-2.0D+0-SIGMAR,-SIGMAI)
         C2 = DCMPLX(2.0D+0-4.0D+0*SIGMAR,-4.0D+0*SIGMAI)
         C3 = DCMPLX(3.0D+0-SIGMAR,-SIGMAI)
*
         DO 20 J = 1, N - 1
            CDL(J) = C1
            CDD(J) = C2
            CDU(J) = C3
   20    CONTINUE
         CDD(N) = C2
*
         CALL ZGTTRF(N,CDL,CDD,CDU,CDU2,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
               CALL MV(N,X)
               DO 60 J = 1, N
                  CTEMP(J) = DCMPLX(X(J))
   60          CONTINUE
               CALL ZGTTRS('N',N,1,CDL,CDD,CDU,CDU2,IPIV,CTEMP,N,INFO)
               DO 80 J = 1, N
                  X(J) = DBLE(CTEMP(J))
   80          CONTINUE
            ELSE IF (IREVCM.EQ.1) THEN
*              Perform  x <--- OP*x = inv[A-SIGMA*M]*M*x,
*              M*X stored in MX.
               DO 100 J = 1, N
                  CTEMP(J) = DCMPLX(MX(J))
  100          CONTINUE
               CALL ZGTTRS('N',N,1,CDL,CDD,CDU,CDU2,IPIV,CTEMP,N,INFO)
               DO 120 J = 1, N
                  X(J) = DBLE(CTEMP(J))
  120          CONTINUE
            ELSE IF (IREVCM.EQ.2) THEN
*              Perform  y <--- M*x
               CALL MV(N,X)
            ELSE IF (IREVCM.EQ.4) 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)
            FIRST = .TRUE.
            DO 140 J = 1, NCONV
*              Use Rayleigh Quotient to recover eigenvalues of the original
*              problem.
               IF (D(J,2).EQ.ZERO) THEN
*                 Ritz value is real.
                  CALL AV(N,V(1,J),AX)
                  NUMR = DDOT(N,V(1,J),1,AX,1)
                  CALL DCOPY(N,V(1,J),1,MX,1)
                  CALL MV(N,MX)
                  DENR = DDOT(N,V(1,J),1,MX,1)
                  D(J,1) = NUMR/DENR
               ELSE IF (FIRST) THEN
*                 Ritz value is complex.
*                 Compute x'(Ax)
                  CALL AV(N,V(1,J),AX)
                  NUMR = DDOT(N,V(1,J),1,AX,1)
                  NUMI = DDOT(N,V(1,J+1),1,AX,1)
                  CALL AV(N,V(1,J+1),AX)
                  NUMR = NUMR + DDOT(N,V(1,J+1),1,AX,1)
                  NUMI = -NUMI + DDOT(N,V(1,J),1,AX,1)
*                 Compute x'(Mx)
                  CALL DCOPY(N,V(1,J),1,MX,1)
                  CALL MV(N,MX)
                  DENR = DDOT(N,V(1,J),1,MX,1)
                  DENI = DDOT(N,V(1,J+1),1,MX,1)
                  CALL DCOPY(N,V(1,J+1),1,MX,1)
                  CALL MV(N,MX)
                  DENR = DENR + DDOT(N,V(1,J+1),1,MX,1)
                  DENI = -DENI + DDOT(N,V(1,J),1,MX,1)
*              d=x'(Ax)/x'(Mx)
                  D(J,1) = (NUMR*DENR+NUMI*DENI)/F06BNF(DENR,DENI)
                  D(J,2) = (NUMI*DENR-NUMR*DENI)/F06BNF(DENR,DENI)
                  FIRST = .FALSE.
               ELSE
*              Second of complex conjugate pair.
                  D(J,1) = D(J-1,1)
                  D(J,2) = -D(J-1,2)
                  FIRST = .TRUE.
               END IF
  140       CONTINUE
*        Print computed eigenvalues.
            WRITE (NOUT,99996) NCONV, SIGMAR, SIGMAI
            DO 160 J = 1, NCONV
               WRITE (NOUT,99995) J, D(J,1), D(J,2)
  160       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 =',E12.4)
99997 FORMAT (1X,' NAG Routine F12ABF Returned with IFAIL = ',I6)
99996 FORMAT (1X,/' The ',I4,' generalized Ritz values closest to ','(',
     +       F8.4,', ',F8.4,')',' are:',/)
99995 FORMAT (1X,I8,5X,'(',F7.4,',',F7.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 FOUR
      PARAMETER     (FOUR=4.0D+0)
*     .. Scalar Arguments ..
      INTEGER       N
*     .. Array Arguments ..
      DOUBLE PRECISION V(N)
*     .. Local Scalars ..
      DOUBLE PRECISION VM1, VV
      INTEGER       J
*     .. Executable Statements ..
      VM1 = V(1)
      V(1) = FOUR*V(1) + V(2)
      DO 20 J = 2, N - 1
         VV = V(J)
         V(J) = VM1 + FOUR*VV + V(J+1)
         VM1 = VV
   20 CONTINUE
      V(N) = VM1 + FOUR*V(N)
      RETURN
      END
*
      SUBROUTINE AV(N,V,W)
*     .. Parameters ..
      DOUBLE PRECISION THREE, TWO
      PARAMETER     (THREE=3.0D+0,TWO=2.0D+0)
*     .. Scalar Arguments ..
      INTEGER       N
*     .. Array Arguments ..
      DOUBLE PRECISION V(N), W(N)
*     .. Local Scalars ..
      INTEGER       J
*     .. Executable Statements ..
      W(1) = TWO*V(1) + THREE*V(2)
      DO 20 J = 2, N - 1
         W(J) = -TWO*V(J-1) + TWO*V(J) + THREE*V(J+1)
   20 CONTINUE
      W(N) = -TWO*V(N-1) + TWO*V(N)
      RETURN
      END
