C
C ------------------------------------------------------------------
C
      SUBROUTINE DIPPLA(OP, IOVECT, N, XL, XR, NP, NPERM, NOP,
     *   NMVAL, VAL, NMVEC, VEC, NBLOCK, H, HV, P, Q, BOUND,
     *   D, DELTAL, DELTAR, RARITZ, EPS, IERR)
C
      INTEGER N, NP, NPERM, NOP, NMVAL, NMVEC, NBLOCK, IERR
      LOGICAL RARITZ
      DOUBLE PRECISION XL, XR, VAL(NMVAL,1), VEC(NMVEC,1),
     *   H(NP,1), HV(NP, 1), P(N,1), Q(N,1), BOUND(1), D(1),
     *   DELTAL, DELTAR, EPS
      EXTERNAL OP, IOVECT
C
C THIS SUBROUTINE POST PROCESSES THE EIGENVECTORS.  BLOCK MATRIX
C VECTOR PRODUCTS ARE USED TO MINIMIZED THE NUMBER OF CALLS TO OP.
C
      INTEGER I, J, JJ, K, KK, L, M, MOD
      DOUBLE PRECISION HMAX, HMIN, TEMP, DDOT, DNRM2, DABS,
     *   DMAX1, DMIN1, DSQRT, DZERO(1)
      EXTERNAL DAXPY, DCOPY, DDOT, DLAGER, DLAEIG
      DZERO(1) = 0.0D0
C
C IF RARITZ IS .TRUE.  A FINAL RAYLEIGH-RITZ PROCEDURE IS
C APPLIED TO THE EIGENVECTORS.
C
      IF (.NOT.RARITZ) GO TO 180
C
C ------------------------------------------------------------------
C
C THIS CONSTRUCTS H=Q*AQ, WHERE THE COLUMNS OF Q ARE THE
C APPROXIMATE EIGENVECTORS.
C
      CALL DCOPY(NPERM*NPERM, DZERO, 0, H, 1)
      M = MOD(NPERM,NBLOCK)
      IF (M.EQ.0) GO TO 40
      DO 10 I=1,M
         CALL DCOPY(N, VEC(1,I), 1, P(1,I), 1)
   10 CONTINUE
      CALL IOVECT(N, M, P, M, 0)
      CALL OP(N, M, P, Q)
      NOP = NOP + 1
      DO 30 I=1,M
         DO 20 J=I,NPERM
            JJ = J - I + 1
            H(JJ,I) = DDOT(N,VEC(1,J),1,Q(1,I),1)
   20    CONTINUE
   30 CONTINUE
      IF (NPERM.LT.NBLOCK) GO TO 90
   40 M = M + NBLOCK
      DO 80 I=M,NPERM,NBLOCK
         DO 50 J=1,NBLOCK
            L = I - NBLOCK + J
            CALL DCOPY(N, VEC(1,L), 1, P(1,J), 1)
   50    CONTINUE
         CALL IOVECT(N, NBLOCK, P, I, 0)
         CALL OP(N, NBLOCK, P, Q)
         NOP = NOP + 1
         DO 70 J=1,NBLOCK
            L = I - NBLOCK + J
            DO 60 K=L,NPERM
               KK = K - L + 1
               H(KK,L) = DDOT(N,VEC(1,K),1,Q(1,J),1)
   60       CONTINUE
   70    CONTINUE
   80 CONTINUE
C
C THIS COMPUTES THE SPECTRAL DECOMPOSITION OF H.
C
   90 HMIN = H(1,1)
      HMAX = H(1,1)
      CALL DLAGER(NPERM, NPERM, 1, H, HMIN, HMAX)
      CALL DLAEIG(NPERM, NPERM, 1, NPERM, H, VAL, NPERM, HV,
     1  BOUND, P, D, Q, EPS, HMIN, HMAX)
C
C THIS COMPUTES THE RITZ VECTORS--THE COLUMNS OF
C Y = QS WHERE S IS THE MATRIX OF EIGENVECTORS OF H.
C
      DO 110 I=1,NPERM
         CALL DCOPY(N, DZERO, 0, VEC(1,I), 1)
  110 CONTINUE
      M = MOD(NPERM,NBLOCK)
      IF (M.EQ.0) GO TO 140
      CALL IOVECT(N, M, P, M, 1)
      DO 130 I=1,M
         DO 120 J=1,NPERM
            CALL DAXPY(N, HV(I,J), P(1,I), 1, VEC(1,J), 1)
  120    CONTINUE
  130 CONTINUE
      IF (NPERM.LT.NBLOCK) GO TO 180
  140 M = M + NBLOCK
      DO 170 I=M,NPERM,NBLOCK
         CALL IOVECT(N, NBLOCK, P, I, 1)
         DO 160 J=1,NBLOCK
            L = I - NBLOCK + J
            DO 150 K=1,NPERM
               CALL DAXPY(N, HV(L,K), P(1,J), 1, VEC(1,K), 1)
  150       CONTINUE
  160    CONTINUE
  170 CONTINUE
C
C ------------------------------------------------------------------
C
C THIS SECTION COMPUTES THE RAYLEIGH QUOTIENTS (IN VAL(*,1))
C AND RESIDUAL NORMS (IN VAL(*,2)) OF THE EIGENVECTORS.
C
  180 M = MOD(NPERM,NBLOCK)
      IF (M.EQ.0) GO TO 210
      DO 190 I=1,M
         CALL DCOPY(N, VEC(1,I), 1, P(1,I), 1)
  190 CONTINUE
      CALL OP(N, M, P, Q)
      NOP = NOP + 1
      DO 200 I=1,M
         VAL(I,1) = DDOT(N,P(1,I),1,Q(1,I),1)
         CALL DAXPY(N, -VAL(I,1), P(1,I), 1, Q(1,I), 1)
         VAL(I,2) = DNRM2(N,Q(1,I),1)
  200 CONTINUE
      IF (NPERM.LT.NBLOCK) GO TO 250
  210 M = M + 1
      DO 240 I=M,NPERM,NBLOCK
         DO 220 J=1,NBLOCK
            L = I - 1 + J
            CALL DCOPY(N, VEC(1,L), 1, P(1,J), 1)
  220    CONTINUE
         CALL OP(N, NBLOCK, P, Q)
         NOP = NOP + 1
         DO 230 J=1,NBLOCK
            L = I - 1 + J
            VAL(L,1) = DDOT(N,P(1,J),1,Q(1,J),1)
            CALL DAXPY(N, -VAL(L,1), P(1,J), 1, Q(1,J), 1)
            VAL(L,2) = DNRM2(N,Q(1,J),1)
  230    CONTINUE
  240 CONTINUE
C
C THIS COMPUTES THE ACCURACY ESTIMATES.  IF AN EIGENVALUE'S RESIDUAL
C NORM INTERVAL LIES INSIDE THE EXCLUDED INTERVAL, THE EIGENPAIR
C IS DELETED.  IF AN EIGENVALUE IS IN THE EXCLUDED INTERVAL BUT ITS
C RESIDUAL NORM INTERVAL OVERLAPS THE BOUNDARY THEN THE EIGENVALUE
C IS SET EQUAL TO THE BOUNDARY, THE RESIDUAL NORM IS RECOMPUTED AND
C IS MADE NEGATIVE, AND 10 IS SUBTRACTED FROM IERR.  A DO LOOP WAS
C NOT USED TO ACCOMODATE THE POSSIBLE DELETION OF EIGENPAIRS.
C
  250 I = 0
  260 I = I + 1
      IF (I.GT.NPERM) RETURN
  270 IF (VAL(I,1).LE.XL .OR. VAL(I,1).GE.XR) GO TO 310
      IF (VAL(I,1).GE.DELTAR) GO TO 280
      IF (VAL(I,1)-VAL(I,2).GT.XL) GO TO 290
      TEMP = VAL(I,1) - XL
      VAL(I,1) = XL
      VAL(I,2) = -DSQRT(TEMP*TEMP+VAL(I,2)*VAL(I,2))
      IERR = IERR - 10
      GO TO 310
C
  280 IF (VAL(I,1)+VAL(I,2).LT.XR) GO TO 290
      TEMP = XR - VAL(I,1)
      VAL(I,1) = XR
      VAL(I,2) = -DSQRT(TEMP*TEMP+VAL(I,2)*VAL(I,2))
      IERR = IERR - 10
      GO TO 310
C
  290 NPERM = NPERM - 1
      DELTAL = DMIN1(DELTAL,VAL(I,1))
      DELTAR = DMAX1(DELTAR,VAL(I,1))
      IF (I.GT.NPERM) RETURN
      DO 300 J=I,NPERM
         CALL DCOPY(N, VEC(1,J+1), 1, VEC(1,J), 1)
         VAL(J,1) = VAL(J+1,1)
         VAL(J,2) = VAL(J+1,2)
  300 CONTINUE
      GO TO 270
C
  310 TEMP = DMAX1(DELTAL-VAL(I,1),VAL(I,1)-DELTAR)
      VAL(I,4) = 0.0D0
      IF (TEMP.GT.0.0D0) VAL(I,4) = DABS(VAL(I,2))/TEMP
      VAL(I,3) = VAL(I,4)*DABS(VAL(I,2))
      GO TO 260
C
      END

.