C
C***********************************************************************
C
      SUBROUTINE DLABCM(N, NBAND, NL, NR, A, EIGVAL,
     1  LDE, EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP)
C
C  THIS SUBROUTINE ORGANIZES THE CALCULATION OF THE EIGENVALUES
C  FOR THE BNDEIG PACKAGE.  EIGENVALUES ARE COMPUTED BY 
C  A MODIFIED RAYLEIGH QUOTIENT ITERATION.  THE EIGENVALUE COUNT
C  OBTAINED BY EACH FACTORIZATION IS USED TO OCCASIONALLY OVERRIDE
C  THE COMPUTED RAYLEIGH QUOTIENT WITH A DIFFERENT SHIFT TO 
C  INSURE CONVERGENCE TO THE DESIRED EIGENVALUES.
C
C  FORMAL PARAMETERS.
C
      INTEGER N, NBAND, NL, NR, LDE
      DOUBLE PRECISION A(NBAND,1), EIGVAL(1),
     1  EIGVEC(LDE,1), ATOL, ARTOL, BOUND(2,1), ATEMP(1),
     2  D(1), VTEMP(1)
C
C
C  LOCAL VARIABLES
C
      LOGICAL FLAG
      INTEGER I, J, L, M, NUML, NUMVEC, NVAL
      DOUBLE PRECISION ERRB, GAP, RESID, RQ, SIGMA, VNORM
C
C
C  FUNCTIONS CALLED
C
      INTEGER MIN0
      DOUBLE PRECISION DMAX1, DMIN1, DDOT, DNRM2
C
C  SUBROUTINES CALLED
C
C     DLABAX, DLABFC, DLARAN, DAXPY, DCOPY, DSCAL
C
C  REPLACE ZERO VECTORS BY RANDOM
C
      NVAL = NR - NL + 1
      FLAG = .FALSE.
      DO 5 I = 1, NVAL
         IF(DDOT(N, EIGVEC(1,I), 1, EIGVEC(1,I), 1) .EQ. 0.0D0)
     1      CALL DLARAN(N,EIGVEC(1,I))
    5 CONTINUE
C
C  LOOP OVER EIGENVALUES
C
      SIGMA = BOUND(2,NVAL+1)
      DO 400 J = 1, NVAL
         NUML = J
C
C  PREPARE TO COMPUTE FIRST RAYLEIGH QUOTIENT
C
   10    CALL DLABAX(N, NBAND, A, EIGVEC(1,J), VTEMP)
         VNORM = DNRM2(N, VTEMP, 1)
         IF(VNORM .EQ. 0.0D0) GO TO 20
         CALL DSCAL(N, 1.0D0/VNORM, VTEMP, 1)
         CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
         CALL DAXPY(N, -SIGMA, EIGVEC(1,J), 1, VTEMP, 1)
C
C  LOOP OVER SHIFTS
C
C  COMPUTE RAYLEIGH QUOTIENT, RESIDUAL NORM, AND CURRENT TOLERANCE
C
   20       VNORM = DNRM2(N, EIGVEC(1,J), 1)
            IF(VNORM .NE. 0.0D0) GO TO 30
            CALL DLARAN(N, EIGVEC(1,J))
            GO TO 10
C
   30       RQ = SIGMA + DDOT(N, EIGVEC(1,J), 1, VTEMP, 1)
     1        /VNORM/VNORM
            CALL DAXPY(N, SIGMA-RQ, EIGVEC(1,J), 1, VTEMP, 1)
            RESID = DMAX1(ATOL, DNRM2(N, VTEMP, 1)/VNORM)
            CALL DSCAL(N, 1.0/VNORM, EIGVEC(1,J), 1)
C
C  ACCEPT EIGENVALUE IF THE INTERVAL IS SMALL ENOUGH
C
            IF(BOUND(2,J+1) - BOUND(1,J+1) .LT. 3.0D0*ATOL) GO TO 300
C
C  COMPUTE MINIMAL ERROR BOUND
C
            ERRB = RESID
            GAP = DMIN1(BOUND(1,J+2) - RQ, RQ - BOUND(2,J))
            IF(GAP .GT. RESID) ERRB = DMAX1(ATOL, RESID*RESID/GAP)
C
C  TENTATIVE NEW SHIFT
C
            SIGMA = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1))
C
C  CHECK FOR TERMINALTION
C
            IF(RESID .GT. 2.0D0*ATOL) GO TO 40
            IF(RQ - ERRB .GT. BOUND(2,J) .AND. 
     1        RQ + ERRB .LT. BOUND(1,J+2)) GO TO 310
C
C  RQ IS TO THE LEFT OF THE INTERVAL
C
   40       IF(RQ .GE. BOUND(1,J+1)) GO TO 50
            IF(RQ - ERRB .GT. BOUND(2,J)) GO TO 100
            IF(RQ + ERRB .LT. BOUND(1,J+1)) CALL DLARAN(N,EIGVEC(1,J))
            GO TO 200
C
C  RQ IS TO THE RIGHT OF THE INTERVAL
C
   50       IF(RQ .LE. BOUND(2,J+1)) GO TO 100
            IF(RQ + ERRB .LT. BOUND(1,J+2)) GO TO 100
C
C  SAVE THE REJECTED VECTOR IF INDICATED
C
            IF(RQ - ERRB .LE. BOUND(2,J+1)) GO TO 200
            DO 60 I = J, NVAL
               IF(BOUND(2,I+1) .GT. RQ) GO TO 70
   60       CONTINUE
            GO TO 80
C
   70       CALL DCOPY(N, EIGVEC(1,J), 1, EIGVEC(1,I), 1)
C
   80       CALL DLARAN(N, EIGVEC(1,J))
            GO TO 200
C
C  PERTURB RQ TOWARD THE MIDDLE
C
  100       IF(SIGMA .LT. RQ) SIGMA = DMAX1(SIGMA, RQ-ERRB)
            IF(SIGMA .GT. RQ) SIGMA = DMIN1(SIGMA, RQ+ERRB)
C
C  FACTOR AND SOLVE
C
  200       DO 210 I = J, NVAL
               IF(SIGMA .LT. BOUND(1,I+1)) GO TO 220
  210       CONTINUE
            I = NVAL + 1
  220       NUMVEC = I - J
            NUMVEC = MIN0(NUMVEC, NBAND + 2)
            IF(RESID .LT. ARTOL) NUMVEC = MIN0(1,NUMVEC)
            CALL DCOPY(N, EIGVEC(1,J), 1, VTEMP, 1)
            CALL DLABFC(N, NBAND, A, SIGMA, NUMVEC, LDE, 
     1        EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL)
C
C  PARTIALLY SCALE EXTRA VECTORS TO PREVENT UNDERFLOW OR OVERFLOW
C
            IF(NUMVEC .EQ. 1) GO TO 227
            L = NUMVEC - 1 
            DO 225 I = 1,L
               M = J + I
               CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,M), 1)
  225       CONTINUE
C
C  UPDATE INTERVALS
C
  227       NUML = NUML - NL + 1
            IF(NUML .GE. 0) BOUND(2,1) = DMIN1(BOUND(2,1), SIGMA)
            DO 230 I = J, NVAL
               IF(SIGMA .LT. BOUND(1,I+1)) GO TO 20
               IF(NUML .LT. I) BOUND(1,I+1) = SIGMA
               IF(NUML .GE. I) BOUND(2,I+1) = SIGMA
  230       CONTINUE
            IF(NUML .LT. NVAL + 1) BOUND(1,NVAL+2) = DMAX1(SIGMA,
     1        BOUND(1,NVAL+2))
            GO TO 20
C
C  ACCEPT AN EIGENPAIR
C
  300    CALL DLARAN(N, EIGVEC(1,J))
         FLAG = .TRUE.
         GO TO 310
C
  305    FLAG = .FALSE.
         RQ = 0.5D0*(BOUND(1,J+1) + BOUND(2,J+1))
         CALL DLABFC(N, NBAND, A, RQ, NUMVEC, LDE, 
     1    EIGVEC(1,J), NUML, 2*NBAND-1, ATEMP, D, ATOL)
         VNORM = DNRM2(N, EIGVEC(1,J), 1)
         IF(VNORM .NE. 0.0) CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
C
C  ORTHOGONALIZE THE NEW EIGENVECTOR AGAINST THE OLD ONES
C
  310    EIGVAL(J) = RQ
         IF(J .EQ. 1) GO TO 330
         M = J - 1
         DO 320 I = 1, M
            CALL DAXPY(N, -DDOT(N,EIGVEC(1,I),1,EIGVEC(1,J),1),
     1        EIGVEC(1,I), 1, EIGVEC(1,J), 1)
  320    CONTINUE
  330    VNORM = DNRM2(N, EIGVEC(1,J), 1)
         IF(VNORM .EQ. 0.0D0) GO TO 305
         CALL DSCAL(N, 1.0D0/VNORM, EIGVEC(1,J), 1)
C
C   ORTHOGONALIZE LATER VECTORS AGAINST THE CONVERGED ONE
C
         IF(FLAG) GO TO 305
         IF(J .EQ. NVAL) RETURN
         M = J + 1
         DO 340 I = M, NVAL
            CALL DAXPY(N, -DDOT(N,EIGVEC(1,J),1,EIGVEC(1,I),1),
     1        EIGVEC(1,J), 1, EIGVEC(1,I), 1)
  340    CONTINUE
  400 CONTINUE
      RETURN
C
  500 CONTINUE
      END

.