*     F12FBF Example Program Text
*     Mark 21 Release. NAG Copyright 2004.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. External Subroutines ..
      EXTERNAL         EX1, EX2
*     .. Executable Statements ..
      WRITE (NOUT,*) 'F12FBF Example Program Results'
      CALL EX1
      CALL EX2
      STOP
      END
*
      SUBROUTINE EX1
*     .. Parameters ..
      INTEGER          IMON, LICOMM, NIN, NOUT
      PARAMETER        (IMON=0,LICOMM=134,NIN=5,NOUT=6)
      INTEGER          MAXN, MAXNCV, LDV
      PARAMETER        (MAXN=256,MAXNCV=30,LDV=MAXN)
      INTEGER          LCOMM
      PARAMETER        (LCOMM=3*MAXN+MAXNCV*MAXNCV+8*MAXNCV+60)
      DOUBLE PRECISION ONE, TWO, ZERO
      PARAMETER        (ONE=1.0D+0,TWO=2.0D+0,ZERO=0.0D+0)
*     .. Local Scalars ..
      DOUBLE PRECISION H2, SIGMA
      INTEGER          IFAIL, INFO, IREVCM, J, N, NCONV, NCV, NEV,
     +                 NITER, NSHIFT
*     .. Local Arrays ..
      DOUBLE PRECISION AD(MAXN), ADL(MAXN), ADU(MAXN), ADU2(MAXN),
     +                 COMM(LCOMM), D(MAXNCV,2), 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, F12FAF, F12FBF, F12FCF,
     +                 F12FDF, F12FEF
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'F12FBF Example 1'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*)
      READ (NIN,*) N, NEV, NCV
      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 F12FAF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*     Set the region of the spectrum that is required.
         CALL F12FDF('LARGEST MAGNITUDE',ICOMM,COMM,IFAIL)
*     Use the Shifted Inverse mode.
         CALL F12FDF('SHIFTED INVERSE',ICOMM,COMM,IFAIL)
*
         H2 = ONE/DBLE((N+1)*(N+1))
         SIGMA = ZERO
         DO 20 J = 1, N
            AD(J) = TWO/H2 - SIGMA
            ADL(J) = -ONE/H2
   20    CONTINUE
         CALL DCOPY(N,ADL,1,ADU,1)
         CALL DGTTRF(N,ADL,AD,ADU,ADU2,IPIV,INFO)
*
         IREVCM = 0
         IFAIL = -1
   40    CONTINUE
         CALL F12FBF(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 matrix vector multiplication y <--- inv[A-sigma*I]*x
               CALL DGTTRS('N',N,1,ADL,AD,ADU,ADU2,IPIV,X,N,INFO)
            ELSE IF (IREVCM.EQ.4 .AND. IMON.NE.0) THEN
*           Output monitoring information
               CALL F12FEF(NITER,NCONV,D,D(1,2),ICOMM,COMM)
               WRITE (6,99998) NITER, NCONV, DNRM2(NEV,D(1,2),1)
            END IF
            GO TO 40
         END IF
         IF (IFAIL.EQ.0) THEN
*        Post-Process using F12FCF to compute eigenvalues/vectors.
            CALL F12FCF(NCONV,D,V,LDV,SIGMA,RESID,V,LDV,COMM,ICOMM,
     +                  IFAIL)
            WRITE (NOUT,99996) NCONV, SIGMA
            DO 60 J = 1, NCONV
               WRITE (NOUT,99995) J, D(J,1)
   60       CONTINUE
         ELSE
            WRITE (NOUT,99997) IFAIL
         END IF
      END IF
      RETURN
*
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 F12FBF Returned with IFAIL = ',I6)
99996 FORMAT (1X,/' The ',I4,' Ritz values of closest to ',F8.4,' are:',
     +       /)
99995 FORMAT (1X,I8,5X,F12.4)
      END
*
      SUBROUTINE EX2
*     .. Parameters ..
      INTEGER          MAXM, MAXN, MAXNEV, MAXNCV, NIN, NOUT
      PARAMETER        (MAXM=500,MAXN=250,MAXNEV=10,MAXNCV=25,NIN=5,
     +                 NOUT=6)
      INTEGER          LCOMM, LDU, LDV, LICOMM
      PARAMETER        (LCOMM=3*MAXN+MAXNCV*MAXNCV+8*MAXNCV+60,
     +                  LDU=MAXM,LDV=MAXN,LICOMM=140)
*     .. Local Scalars ..
      DOUBLE PRECISION SIGMA, TEMP
      INTEGER          IFAIL, IFAIL1, IREVCM, J, M, N, NCONV, NCV, NEV,
     +                 NSHIFT
*     .. External Functions ..
      DOUBLE PRECISION DNRM2
      EXTERNAL         DNRM2
*     .. External Subroutines ..
      EXTERNAL         ATV, AV, DAXPY, DCOPY, DSCAL, F12FAF, F12FBF,
     +                 F12FCF, X04CAF
*     .. Local Arrays ..
      DOUBLE PRECISION AX(MAXM), COMM(LCOMM), D(MAXNCV,2), MX(MAXM),
     +                 RESID(MAXN), U(LDU,MAXNEV), V(LDV,MAXNCV),
     +                 X(MAXM)
      INTEGER          ICOMM(LICOMM)
*     .. Intrinsic Functions ..
      INTRINSIC        DSQRT
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'F12FBF Example 2'
      WRITE (NOUT,*)
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*) M, N, NEV, NCV
      IF (N.LT.1 .OR. N.GT.MAXN) THEN
         WRITE (NOUT,99999) 'N is out of range: N = ', N
      ELSE IF (M.LT.1 .OR. M.GT.MAXM) THEN
         WRITE (NOUT,99999) 'M is out of range: M = ', M
      ELSE IF (NCV.GT.MAXNCV) THEN
         WRITE (NOUT,99999) 'NCV is out of range: NCV = ', NCV
      ELSE
*        Initialize for problem.
         IFAIL = 0
         CALL F12FAF(N,NEV,NCV,ICOMM,LICOMM,COMM,LCOMM,IFAIL)
*        Main reverse communication loop.
         IREVCM = 0
   20    CONTINUE
         CALL F12FBF(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 the operation X <- A'AX
               CALL AV(M,N,X,AX)
               CALL ATV(M,N,AX,X)
            END IF
            GO TO 20
         END IF
         IF (IFAIL.EQ.0) THEN
*           Post-Process using F12FCF.
*           Computed singular values may be extracted.
*           Singular vectors may also be computed now if desired.
            IFAIL1 = 0
            CALL F12FCF(NCONV,D,V,LDV,SIGMA,RESID,V,LDV,COMM,ICOMM,
     *                  IFAIL1)
*
*           Singular values (squared) are returned in the first column
*           of D and the corresponding right singular vectors are
*           returned in the first NEV columns of V.
*
            DO 40 J = 1, NCONV
               D(J,1) = DSQRT(D(J,1))
*
*              Compute the left singular vectors from the formula
*                     u = Av/sigma
*              u should have norm 1 so divide by norm(Av).
*
               CALL AV(M,N,V(1,J),AX)
               CALL DCOPY(M,AX,1,U(1,J),1)
               TEMP = 1.0D0/DNRM2(M,U(1,J),1)
               CALL DSCAL(M,TEMP,U(1,J),1)
*
*              Compute the residual norm ||A*v - sigma*u|| for the NCONV
*              accurately computed singular values and vectors.
*              Store the result in 2nd column of array D.
*
               CALL DAXPY(M,-D(J,1),U(1,J),1,AX,1)
               D(J,2) = DNRM2(M,AX,1)
*
   40       CONTINUE
*
*           Print computed residuals
*
            CALL X04CAF('G','N',NCONV,2,D,MAXNCV,
     +                  ' Singular values and direct residuals',IFAIL1)
         ELSE
            WRITE (NOUT,99999) IFAIL
         END IF
      END IF
      RETURN
*
99999 FORMAT (1X,' NAG Routine F12FBF Returned with IFAIL = ',I6)
      END
*
*     Matrix vector subroutines
*
      SUBROUTINE AV(M,N,X,W)
*
*     Computes  w <- A*x.
*
*     .. Parameters ..
      DOUBLE PRECISION ONE, ZERO
      PARAMETER        (ONE=1.0D+0,ZERO=0.0D+0)
*     .. Scalar Arguments ..
      INTEGER          M, N
*     .. Array Arguments ..
      DOUBLE PRECISION W(M), X(N)
*     .. Local Scalars ..
      DOUBLE PRECISION H, K, S, T
      INTEGER          I, J
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      CONTINUE
      H = ONE/DBLE(M+1)
      K = ONE/DBLE(N+1)
      DO 20 I = 1, M
         W(I) = ZERO
   20 CONTINUE
      T = ZERO
*
      DO 80 J = 1, N
         T = T + K
         S = ZERO
         DO 40 I = 1, J
            S = S + H
            W(I) = W(I) + K*S*(T-ONE)*X(J)
   40    CONTINUE
         DO 60 I = J + 1, M
            S = S + H
            W(I) = W(I) + K*T*(S-ONE)*X(J)
   60    CONTINUE
   80 CONTINUE
*
      RETURN
      END
*
      SUBROUTINE ATV(M,N,W,Y)
*
*     Computes  y <- A'*w.
*
*     .. Parameters ..
      DOUBLE PRECISION ONE, ZERO
      PARAMETER        (ONE=1.0D+0,ZERO=0.0D+0)
*     .. Scalar Arguments ..
      INTEGER          M, N
*     .. Array Arguments ..
      DOUBLE PRECISION W(M), Y(N)
*     .. Local Scalars ..
      DOUBLE PRECISION H, K, S, T
      INTEGER          I, J
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Executable Statements ..
      CONTINUE
      H = ONE/DBLE(M+1)
      K = ONE/DBLE(N+1)
      DO 20 I = 1, N
         Y(I) = ZERO
   20 CONTINUE
      T = ZERO
*
      DO 80 J = 1, N
         T = T + K
         S = ZERO
         DO 40 I = 1, J
            S = S + H
            Y(J) = Y(J) + K*S*(T-ONE)*W(I)
   40    CONTINUE
         DO 60 I = J + 1, M
            S = S + H
            Y(J) = Y(J) + K*T*(S-ONE)*W(I)
   60    CONTINUE
   80 CONTINUE
*
      RETURN
      END