C
C ------------------------------------------------------------------
C
      SUBROUTINE DIWLA(OP, IOVECT, N, NBAND, XL, XR, NFIG, NPERM, VAL,
     *   NMVEC, MAXVEC, VEC, NBLOCK, MAXOP, MAXJ, NOP, P1,
     *   P0, RES, TAU, OTAU, T, ALP, BET, S, P2, BOUND, ATEMP, VTEMP,
     *   D, IND, RARITZ, DELTAL, DELTAR, EPS, IERR)
C
      INTEGER N, NBAND, NFIG, NPERM, NMVEC, MAXVEC, NBLOCK, MAXOP,
     *   MAXJ, NOP, IND(1), IERR
      LOGICAL RARITZ
      DOUBLE PRECISION XL, XR, VAL(1), VEC(NMVEC,1), P0(N,1),
     *   P1(N,1), P2(N,1), RES(1), TAU(1), OTAU(1),
     *   T(NBAND,1), ALP(NBLOCK,1), BET(NBLOCK,1), BOUND(1),
     *   ATEMP(1), VTEMP(1), D(1), S(MAXJ,1), DELTAL, DELTAR,
     *   EPS
      EXTERNAL OP, IOVECT
C
C DIWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE
C ORTHOGONALIZATION.
C
C   NPERM   THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS
C     INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE
C     ALGORITHM IS ITERATED).  PERMANENT VECTORS ARE ORTHOGONAL
C     TO THE CURRENT KRYLOV SUBSPACE.
C
C   NOP   THE NUMBER OF CALLS TO OP.
C
C   P0, P1, AND P2   THE CURRENT BLOCKS OF LANCZOS VECTORS.
C
C   RES   THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS.
C
C   TAU AND OTAU   USED TO MONITOR THE NEED FOR ORTHOGONALIZATION.
C
C   T   THE BAND MATRIX.
C
C   ALP   THE CURRENT DIAGONAL BLOCK.
C
C   BET   THE CURRENT OFF DIAGONAL BLOCK.
C
C   BOUND, ATEMP, VTEMP, D  TEMPORARY STORAGE USED BY DLAEIG.
C
C   S   EIGENVECTORS OF T.
C
C   RARITZ   RETURNED AS .TRUE. IF A FINAL RAYLEIGH-RITZ PROCEDURE
C     IS TO BE DONE.
C
C   DELTAL   RETURNED AS THE EXCLUDED EIGENVALUE CLOSEST TO XR.
C     USED TO ESTIMATE THE ACCURACY OF THE COMPUTED EIGENVALUES.
C
C   DELTAR   RETURNED AS THE EXCLUDED EIGENVALUE CLOSEST TO XR.
C     USED TO ESTIMATE THE ACCURACY OF THE COMPUTED EIGENVALUES.
C
C
C INTERNAL VARIABLES USED  
C
      INTEGER I, II, I1, J, K, L, M, NG, NGOOD, NL, NSTART, NTHETA, 
     *   NUMBER, NUML, NUMR, MIN0
      LOGICAL DONE, TEST
      DOUBLE PRECISION ALPMAX, ALPMIN, AXL, ANORM, AXR, BETMAX,
     *   BETMIN, EPSRT, RNORM, TEMP, TMAX, TMIN, TOLA, TOLG, UTOL,
     *   DMAX1, DMIN1, DSQRT, DDOT, DNRM2, DZERO(1), TARR(1)
      EXTERNAL DMVPC, DORTQR, DAXPY, DCOPY, DDOT,
     *   DNRM2, DSCAL, DSWAP, DLARAN, DLAGER, DLAEIG, DVSORT
C
C   J   THE CURRENT DIMENSION OF T.  (THE DIMENSION OF THE CURRENT
C     KRYLOV SUBSPACE)
C
C   NGOOD   THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS
C     LIE IN THE CURRENT KRYLOV SUBSPACE).
C
C   NUMBER = NPERM + NGOOD
C
C   ANORM   AN EXTIMATE OF THE NORM OF THE MATRIX.
C
C   EPS   THE RELATIVE MACHINE PRECISION.
C
C   UTOL   THE USER TOLERANCE.
C
C   DZERO IS AN ARRAY OF SIZE 1 WITH VALUE ZERO.  USED FOR CONSISTENCY OF
C         VARIABLE TYPING IN CALLS TO DCOPY.
C
C   TARR  IS AN ARRAY OF LENGTH 1 USED FOR CONSISTENCE OF VARIABLE TYPING
C         IN CALLS TO DLAEIG.
C
C  INITIALIZE
C
      DZERO(1) = 0.0D0
      RNORM = 0.0D0
      IF (NPERM.NE.0) RNORM = DMAX1(-VAL(1),VAL(NPERM))
      EPSRT = DSQRT(EPS)
      NOP = 0
      NUMBER = NPERM
      RARITZ = .FALSE.
      UTOL = DMAX1(DBLE(FLOAT(N))*EPS,10.0D0**DBLE(-FLOAT(NFIG)))
      TEMP = UTOL*DMAX1(XR, -XL)
      AXL = XL + TEMP
      AXR = DMAX1(AXL,XR-TEMP)
      J = MAXJ
C
C ------------------------------------------------------------------
C
C ANY ITERATION OF THE ALGORITHM BEGINS HERE.
C
   30 DO 50 I=1,NBLOCK
         TEMP = DNRM2(N,P1(1,I),1)
         IF (TEMP.EQ.0.0D0) CALL DLARAN(N,P1(1,I))
   50 CONTINUE
      IF (NPERM.EQ.0) GO TO 70
      DO 60 I=1,NPERM
         TAU(I) = 1.0D0
         OTAU(I) = 0.0D0
   60 CONTINUE
   70 CALL DCOPY(N*NBLOCK, DZERO, 0, P0, 1)
      CALL DCOPY(NBLOCK*NBLOCK, DZERO, 0, BET, 1)
      CALL DCOPY(J*NBAND, DZERO, 0, T, 1)
      M = MAXVEC + 2
      DO 75 I = 1, M
         CALL DCOPY(J, DZERO, 0, S(1,I), 1)
   75 CONTINUE
      TEST = .TRUE.
      NGOOD = 0
      BETMAX = 0.0D0
      TMIN = 1.0D30
      TMAX = -1.0D30
      NUML = 1
      NUMR = 1
      J = 0
C
C ------------------------------------------------------------------
C
C THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP.
C
   80 J = J + NBLOCK
C
C THIS IS THE SELECTIVE ORTHOGONALIZATION.
C
      IF (NUMBER.EQ.0) GO TO 110
      DO 100 I=1,NUMBER
         IF (TAU(I).LT.EPSRT) GO TO 100
         TEST = .TRUE.
         TAU(I) = 0.0D0
         IF (OTAU(I).NE.0.0D0) OTAU(I) = 1.0D0
         DO 90 K=1,NBLOCK
            TEMP = -DDOT(N,VEC(1,I),1,P1(1,K),1)
            CALL DAXPY(N, TEMP, VEC(1,I), 1, P1(1,K), 1)
C
C THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A
C NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR.  THE ALGORITHM IS
C TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST.
C
            IF (DABS(TEMP*BET(K,K)).GT.DBLE(FLOAT(N))*EPSRT*
     *       ANORM .AND. I.GT.NPERM) GO TO 300
   90    CONTINUE
  100 CONTINUE
C
C IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET.
C
  110 IF (.NOT.TEST) GO TO 160
      TEST = .FALSE.
      CALL DORTQR(N, N, NBLOCK, P1, ALP)
      IF(J .EQ. NBLOCK) GO TO 160
      DO 130 I=1,NBLOCK
         IF(ALP(I,I) .GT. 0.0D0) GO TO 130
         M = J - 2*NBLOCK + I
         L = NBLOCK + 1
         DO 120 K=I,NBLOCK
            BET(I,K) = -BET(I,K)
            T(L,M) = -T(L,M)
            M = M + 1
            L = L - 1
  120    CONTINUE
  130 CONTINUE
C
C THIS IS THE LANCZOS STEP.
C
  160 CALL OP(N, NBLOCK, P1, P2)
      NOP = NOP + 1
      CALL IOVECT(N, NBLOCK, P1, J, 0)
C
C THIS COMPUTES P2=P2-P0*BET(TRANSPOSE)
C
      DO 180 I=1,NBLOCK
         DO 170 K=I,NBLOCK
            CALL DAXPY(N, -BET(I,K), P0(1,K), 1, P2(1,I), 1)
  170    CONTINUE
  180 CONTINUE
C
C THIS COMPUTES ALP AND P2=P2-P1*ALP.
C
      DO 200 I=1,NBLOCK
         DO 190 K=1,I
            I1 = I - K + 1
            ALP(I1,K) = DDOT(N,P1(1,I),1,P2(1,K),1)
            CALL DAXPY(N, -ALP(I1,K), P1(1,I), 1, P2(1,K), 1)
            IF (K.NE.I) CALL DAXPY(N, -ALP(I1,K), P1(1,K),
     *       1, P2(1,I), 1)
  190    CONTINUE
  200 CONTINUE
C
C  REORTHOGONALIZATION OF THE SECOND BLOCK
C
      IF(J .NE. NBLOCK) GO TO 220
      DO 215 I=1,NBLOCK
         DO 210 K=1,I
            TEMP = DDOT(N,P1(1,I),1,P2(1,K),1)
            CALL DAXPY(N, -TEMP, P1(1,I), 1, P2(1,K), 1)
            IF (K.NE.I) CALL DAXPY(N, -TEMP, P1(1,K),
     *        1, P2(1,I), 1)
            II = I - K + 1
            ALP(II,K) = ALP(II,K) + TEMP            
  210   CONTINUE
  215 CONTINUE
C
C THIS ORTHONORMALIZES THE NEXT BLOCK.
C
  220 CALL DORTQR(N, N, NBLOCK, P2, BET)
C
C THIS STORES ALP AND BET IN T.
C
      DO 250 I=1,NBLOCK
         M = J - NBLOCK + I
         DO 230 K=I,NBLOCK
            L = K - I + 1
            T(L,M) = ALP(L,I)
  230    CONTINUE
         DO 240 K=1,I
            L = NBLOCK - I + K + 1
            T(L,M) = BET(K,I)
  240    CONTINUE
  250 CONTINUE
C
C THIS SHIFTS THE LANCZOS VECTORS
C
      CALL DCOPY(NBLOCK*N, P1, 1, P0, 1)
      CALL DCOPY(NBLOCK*N, P2, 1, P1, 1)
      CALL DLAGER(J, NBAND, J-NBLOCK+1, T, TMIN, TMAX)
      ANORM = DMAX1(RNORM, TMAX, -TMIN)
      IF(NUMBER .EQ. 0) GO TO 265
C
C THIS COMPUTE THE EXTREME EIGENVALUES OF ALP
C
      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK,
     1  P2, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      ALPMIN = TARR(1)
      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, TARR,
     1  TARR, NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS,
     2  TMIN, TMAX)
      ALPMAX = TARR(1)
C
C  THIS COMPUTES ALP = BET(TRANSPOSE)*BET
C
  265 DO 275 I = 1,NBLOCK
         DO 270 K = 1,I
            L = I - K + 1
            ALP(L,K) = DDOT(NBLOCK-I+1, BET(I,I), NBLOCK, BET(K,I),
     1        NBLOCK)
  270    CONTINUE
  275 CONTINUE
      IF(NUMBER.EQ.0) GO TO 290
C
C  THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET
C
      CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, 1, 1, ALP, TARR, NBLOCK,
     1  P2, BOUND, ATEMP, D, VTEMP, EPS, 0.0D0, ANORM*ANORM)
      BETMIN = DSQRT(TARR(1))
C
C THIS UPDATES TAU AND OTAU.
C
      DO 280 I=1,NUMBER
         TEMP = (TAU(I)*DMAX1(ALPMAX-VAL(I),VAL(I)-ALPMIN) +
     *    OTAU(I)*BETMAX + EPS*ANORM)/BETMIN
         IF (I.LE.NPERM) TEMP = TEMP + RES(I)/BETMIN
         OTAU(I) = TAU(I)
         TAU(I) = TEMP
  280 CONTINUE
C
C THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET
C
  290 CALL DCOPY(NBLOCK, DZERO, 0, P2, 1)
      CALL DLAEIG(NBLOCK, NBLOCK, NBLOCK, NBLOCK, ALP,
     1  TARR, NBLOCK, P2, BOUND, ATEMP, D, VTEMP, EPS,
     2  0.0D0, ANORM*ANORM)
      BETMAX = DSQRT(TARR(1))
      IF (J.LE.2*NBLOCK) GO TO 80
C
      DONE = .TRUE.
      GO TO 310
C
C THIS INDICATES UNEXPECTED LOSS OF ORTHOGONALITY.
C
  300 J = J - NBLOCK
      IERR = -8
C
C ------------------------------------------------------------------
C
C THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO
C SEE IF FURTHER ACTION IS INDICATED.
C
  310 IF(NOP .GE. MAXOP .OR. J .GT. MAXJ-NBLOCK) TEST = .TRUE.
C
C FIRST COMPUTE ALL THE EIGENVALUES LESS THAN AXL AND ONE MORE
C
  320 NL = NUML
      M = NUMBER + NUML
      CALL DLAEIG(J, NBAND, 1, NUML, T, VAL(NUMBER+1),
     1  MAXJ, S, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      IF(VAL(M) .GT. AXL .OR. NUML .EQ. J) GO TO 340
      NUML = NUML + 1
      IF(NUML+NUMR .GT. MAXVEC+2) GO TO 330
C
      DO 325 II = 1, NUMR
         I = NUML + NUMR - II
         CALL DCOPY(J, S(1,I), 1, S(1,I+1), 1)
  325 CONTINUE
      CALL DCOPY(J, DZERO, 0, S(1,NUML), 1)
      GO TO 320
C
C FIND OUT HOW MANY THERE WHEN THERE ARE TOO MANY TO STORE
C
  330 NL = MAXVEC + 2 - NUMR
      M = NUMBER + NL
      CALL DCOPY(J, DZERO, 0, S(1,NL), 1)
      CALL DLAEIG(J, NBAND, NUML, NUML, T, VAL(M),
     1  MAXJ, S(1,NL), BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      IF(VAL(M) .GT. AXL .OR. NUML .EQ. J) GO TO 340
      NUML = NUML + 1
      CALL DCOPY(J, DZERO, 0, S(1,NL), 1)
      GO TO 330
C
C NOW COMPUTE ALL THE EIGENVALUES GREATER THAN AXR AND ONE MORE
C
  340 M = NUMBER + NL + 1
      K = J - NUMR + 1
      CALL DLAEIG(J, NBAND, K, J, T, VAL(M), MAXJ,
     1  S(1, NL+1), BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      IF(VAL(M) .LT. AXR .OR. NUMR .EQ. J) GO TO 370
      NUMR = NUMR + 1
      IF(NUML+NUMR .GT. MAXVEC+2) GO TO 360
      DO 350 II = 2, NUMR
         I = NUML + NUMR - II + 1
         CALL DCOPY(J, S(1,I), 1, S(1,I+1), 1)
  350 CONTINUE
      CALL DCOPY(J, DZERO, 0, S(1,NL+1), 1)
      GO TO 340
C
C FIND OUT HOW MANY THERE ARE WHEN THERE ARE TOO MANY TO STORE
C
  360 K = J - NUMR + 1
      CALL DLAEIG(J, NBAND, K, K, T, VAL(M), MAXJ,
     1  S(1, NL+1), BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
      IF(VAL(M) .LT. AXR .OR. NUMR .EQ. J) GO TO 370
      CALL DCOPY(J, DZERO, 0, S(1,NL+1), 1)
      NUMR = NUMR + 1
      GO TO 360
C
C THIS COMPUTES THE CORRESPONDING RESIDUAL NORMS AND
C ORTHOGONALITY COEFFICIENTS.
C
  370 NTHETA = MIN0(NUML+NUMR, MAXVEC+2)
      CALL DMVPC(NBLOCK, BET, MAXJ, J, S, NTHETA, ATEMP, VTEMP, D)
C
C  THIS CHECKS FOR TOO MANY VECTORS
C
      IF(NUML+NUMR .LE. MAXVEC+2) GO TO 380
      IERR = -100*(NUML+NUMR-MAXVEC-2)
      DONE = .FALSE.
      TEST = .TRUE.
      GO TO 440
C
C THIS CHECKS IF TERMINATION IS ALLOWED
C
  380 M = NUMBER + NL
      DELTAL = VAL(M)
      DELTAR = VAL(M+1)
      IF(J .LT. 6*NBLOCK .OR. DELTAL - ATEMP(NL) .LT. AXL .OR.
     1  DELTAR + ATEMP(NL+1) .GT. AXR) DONE = .FALSE.
      IF(DONE .AND. NTHETA .EQ. 2) RETURN
      IF(NTHETA .EQ. 2 .AND. .NOT. TEST) GO TO 80
C
C ------------------------------------------------------------------
C
C THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL.
C
C THIS CHECKS WHETHER ALL THE EIGENVALUES ARE ACCEPTABLE.
C
  440 IF (.NOT.(TEST .OR. DONE)) GO TO 470
      NSTART = 2
      M = NUMBER + NTHETA
      TOLA = UTOL*DMAX1(RNORM, VAL(M), VAL(NUMBER+1))
      DO 460 I=1,NTHETA
         IF(I .EQ. NL .OR. I .EQ. NL+1) GO TO 460
         M = NUMBER + I
         IF (DMIN1(ATEMP(I)*ATEMP(I)/DMAX1(DELTAL-VAL(M),VAL(M)
     *    -DELTAR),ATEMP(I)).GT.TOLA) GO TO 450
         IND(I) = -1
         GO TO 460
C
  450    DONE = .FALSE.
         IF (.NOT.TEST) GO TO 470
         NSTART = NSTART + 1
         IND(I) = 1
  460 CONTINUE
C
C COPY IND VALUES INTO VTEMP
C
      DO 465 I = 1,NTHETA
         IF(I .EQ. NL .OR. I .EQ. NL+1) GO TO 465
         VTEMP(I) = DBLE(FLOAT(IND(I)))
  465 CONTINUE 
      GO TO 500
C
C THIS CHECKS FOR NEW GOOD VECTORS.
C
  470 NG = 0
      TOLG = EPSRT*ANORM
      DO 490 I=1,NTHETA
         IF(I .EQ. NL .OR. I .EQ. NL+1) GO TO 490
         IF (VTEMP(I).GT.TOLG) GO TO 480
         NG = NG + 1
         VTEMP(I) = -1.0D0
         GO TO 490
  480    VTEMP(I) = 1.0D0
  490 CONTINUE
      IF (NG.LE.NGOOD) GO TO 80
      NSTART = NTHETA - NG
C
C ------------------------------------------------------------------
C
C THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS.
C IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED.
C
  500 TEST = TEST .AND. .NOT.DONE
      NGOOD = MIN0(NTHETA-NSTART,MAXVEC-NPERM)
      NSTART = NTHETA - NGOOD
      VTEMP(NL) = 1.0D0
      VTEMP(NL+1) = 1.0D0
C
C THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND
C EIGENVECTORS OF T.  THE OTHER EIGENVECTORS ARE SAVED FOR
C FORMING STARTING VECTORS, IF NECESSARY.  IT ALSO SHIFTS THE
C EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS
C PAUSE.
C
      CALL DCOPY(NTHETA, VAL(NUMBER+1), 1, VAL(NPERM+1), 1)
      IF (NSTART.EQ.NTHETA) GO TO 530
      CALL DVSORT(NTHETA, VTEMP, ATEMP, 1, VAL(NPERM+1), MAXJ, J, S)
C
C THIS ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING
C VECTORS.
C COMPUTE MINIMUM VALUE OF ATEMP TO AVOID POSSIBLE OVERFLOW
C
  530 IF(.NOT. TEST) NSTART = 0
      IF(.NOT.TEST) GO TO 580
      TEMP = ATEMP(1)
      DO 535 I = 1, NSTART
         TEMP = DMIN1(TEMP, ATEMP(I))
  535 CONTINUE
      M = NGOOD + 1
      L = NGOOD + MIN0(NSTART,NBLOCK)
      DO 540 I=M,L
         CALL DSCAL(J, TEMP/ATEMP(I), S(1,I), 1)
  540 CONTINUE
      M = (NSTART-1)/NBLOCK
      IF (M.EQ.0) GO TO 570
      L = NGOOD + NBLOCK
      DO 560 I=1,M
         DO 550 K=1,NBLOCK
            L = L + 1
            IF (L.GT.NTHETA) GO TO 570
            I1 = NGOOD + K
            CALL DAXPY(J, TEMP/ATEMP(L), S(1,L), 1, S(1,I1), 1)
  550    CONTINUE
  560 CONTINUE
  570 NSTART = MIN0(NSTART,NBLOCK)
C
C THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS.
C
  580 IF (NGOOD.EQ.0 .OR. .NOT.(TEST .OR. DONE)) GO TO 610
      DO 590 I=1,NGOOD
         M = NPERM + I
         RES(M) = ATEMP(I)
  590 CONTINUE
C
C THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE
C LANCZOS VECTORS.
C
  610 NUMBER = NPERM + NGOOD
      IF (TEST .OR. DONE) CALL DCOPY(N*NBLOCK, DZERO, 0, P1, 1)
      NSTART = MIN0(NBLOCK,NSTART)
      IF (NGOOD.EQ.0) GO TO 630
      M = NPERM + 1
      DO 620 I=M,NUMBER
         CALL DCOPY(N, DZERO, 0, VEC(1,I), 1)
  620 CONTINUE
  630 DO 680 I=NBLOCK,J,NBLOCK
         CALL IOVECT(N, NBLOCK, P2, I, 1)
         DO 670 K=1,NBLOCK
            M = I - NBLOCK + K
            IF (NSTART.EQ.0) GO TO 650
            DO 640 L=1,NSTART
               I1 = NGOOD + L
               CALL DAXPY(N, S(M,I1), P2(1,K), 1,
     *          P1(1,L), 1)
  640       CONTINUE
  650       IF (NGOOD.EQ.0) GO TO 670
            DO 660 L=1,NGOOD
               I1 = L + NPERM
               CALL DAXPY(N, S(M,L), P2(1,K), 1, VEC(1,I1), 1)
  660       CONTINUE
  670    CONTINUE
  680 CONTINUE
C
      IF (TEST .OR. DONE) GO TO 700
C
C THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE
C TAU RECURRENCE.
C
      M = NPERM + 1
      DO 690 I=M,NUMBER
         TEMP = 1.0D0/DNRM2(N,VEC(1,I),1)
         CALL DSCAL(N, TEMP, VEC(1,I), 1)
         TAU(I) = 1.0D0
         OTAU(I) = 1.0D0
  690 CONTINUE
      GO TO 80
C
C ------------------------------------------------------------------
C
C THIS SET IERR IF NECESSARY
C
  700 IF(NTHETA .NE. 2 .AND. NOP .GE. MAXOP) IERR = IERR - 2
      IF(NTHETA .EQ. 2 .AND. NOP .GE. MAXOP) IERR =  -6
      IF(NTHETA .EQ. 2 .AND. J .GT. MAXJ-NBLOCK .AND.
     1  IERR .EQ. 0) IERR = -5
C
C THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE
C PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING
C THE PERMANENT VECTORS.
C
      IF (NGOOD.EQ.0 .AND. IERR.EQ.0) GO TO 30
      IF (NGOOD.EQ.0) GO TO 800
C
C THIS ORTHONORMALIZES THE THE EIGENVECTORS
C
      CALL DORTQR(NMVEC,N,NPERM+NGOOD,VEC,S)
C
C THIS SORTS THE VALUES AND VECTORS.
C
      IF(NPERM.NE.0) CALL DVSORT(NPERM+NGOOD, VAL, RES, 0, TARR,
     *  NMVEC, N, VEC)
      NPERM = NPERM + NGOOD
C
C THIS DECIDES WHERE TO GO NEXT.
C
      IF (IERR.NE.0) GO TO 800
      IF (.NOT.DONE) GO TO 30
C
C THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED
C TO LOOK FOR UNDISCLOSED MULTIPLICITIES.
C
      M = NPERM - NBLOCK + 1
      IF (M.LE.0) RETURN
      DO 790 I=1,M
         L = I + NBLOCK - 1
         IF (VAL(L)-VAL(I).LT.TOLA) GO TO 30
  790 CONTINUE
C
C THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ
C PROCEDURE IS NEEDED.
C
  800 M = NPERM - NBLOCK
      IF (M.LE.0) RETURN
      DO 810 I=1,M
         L = I + NBLOCK
         IF (VAL(L)-VAL(I).GE.TOLA) GO TO 810
         RARITZ = .TRUE.
         RETURN
C
  810 CONTINUE
      RETURN
      END

.