SRC/dlaebz.f

   1 *> \brief \b DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than or equal to a given value, and performs other tasks required by the routine sstebz.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLAEBZ + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaebz.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaebz.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaebz.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
  22 *                          RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
  23 *                          NAB, WORK, IWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
  27 *       DOUBLE PRECISION   ABSTOL, PIVMIN, RELTOL
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
  31 *       DOUBLE PRECISION   AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
  32 *      $                   WORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> DLAEBZ contains the iteration loops which compute and use the
  42 *> function N(w), which is the count of eigenvalues of a symmetric
  43 *> tridiagonal matrix T less than or equal to its argument  w.  It
  44 *> performs a choice of two types of loops:
  45 *>
  46 *> IJOB=1, followed by
  47 *> IJOB=2: It takes as input a list of intervals and returns a list of
  48 *>         sufficiently small intervals whose union contains the same
  49 *>         eigenvalues as the union of the original intervals.
  50 *>         The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
  51 *>         The output interval (AB(j,1),AB(j,2)] will contain
  52 *>         eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
  53 *>
  54 *> IJOB=3: It performs a binary search in each input interval
  55 *>         (AB(j,1),AB(j,2)] for a point  w(j)  such that
  56 *>         N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
  57 *>         the search.  If such a w(j) is found, then on output
  58 *>         AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
  59 *>         (AB(j,1),AB(j,2)] will be a small interval containing the
  60 *>         point where N(w) jumps through NVAL(j), unless that point
  61 *>         lies outside the initial interval.
  62 *>
  63 *> Note that the intervals are in all cases half-open intervals,
  64 *> i.e., of the form  (a,b] , which includes  b  but not  a .
  65 *>
  66 *> To avoid underflow, the matrix should be scaled so that its largest
  67 *> element is no greater than  overflow**(1/2) * underflow**(1/4)
  68 *> in absolute value.  To assure the most accurate computation
  69 *> of small eigenvalues, the matrix should be scaled to be
  70 *> not much smaller than that, either.
  71 *>
  72 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
  73 *> Matrix", Report CS41, Computer Science Dept., Stanford
  74 *> University, July 21, 1966
  75 *>
  76 *> Note: the arguments are, in general, *not* checked for unreasonable
  77 *> values.
  78 *> \endverbatim
  79 *
  80 *  Arguments:
  81 *  ==========
  82 *
  83 *> \param[in] IJOB
  84 *> \verbatim
  85 *>          IJOB is INTEGER
  86 *>          Specifies what is to be done:
  87 *>          = 1:  Compute NAB for the initial intervals.
  88 *>          = 2:  Perform bisection iteration to find eigenvalues of T.
  89 *>          = 3:  Perform bisection iteration to invert N(w), i.e.,
  90 *>                to find a point which has a specified number of
  91 *>                eigenvalues of T to its left.
  92 *>          Other values will cause DLAEBZ to return with INFO=-1.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] NITMAX
  96 *> \verbatim
  97 *>          NITMAX is INTEGER
  98 *>          The maximum number of "levels" of bisection to be
  99 *>          performed, i.e., an interval of width W will not be made
 100 *>          smaller than 2^(-NITMAX) * W.  If not all intervals
 101 *>          have converged after NITMAX iterations, then INFO is set
 102 *>          to the number of non-converged intervals.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] N
 106 *> \verbatim
 107 *>          N is INTEGER
 108 *>          The dimension n of the tridiagonal matrix T.  It must be at
 109 *>          least 1.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] MMAX
 113 *> \verbatim
 114 *>          MMAX is INTEGER
 115 *>          The maximum number of intervals.  If more than MMAX intervals
 116 *>          are generated, then DLAEBZ will quit with INFO=MMAX+1.
 117 *> \endverbatim
 118 *>
 119 *> \param[in] MINP
 120 *> \verbatim
 121 *>          MINP is INTEGER
 122 *>          The initial number of intervals.  It may not be greater than
 123 *>          MMAX.
 124 *> \endverbatim
 125 *>
 126 *> \param[in] NBMIN
 127 *> \verbatim
 128 *>          NBMIN is INTEGER
 129 *>          The smallest number of intervals that should be processed
 130 *>          using a vector loop.  If zero, then only the scalar loop
 131 *>          will be used.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] ABSTOL
 135 *> \verbatim
 136 *>          ABSTOL is DOUBLE PRECISION
 137 *>          The minimum (absolute) width of an interval.  When an
 138 *>          interval is narrower than ABSTOL, or than RELTOL times the
 139 *>          larger (in magnitude) endpoint, then it is considered to be
 140 *>          sufficiently small, i.e., converged.  This must be at least
 141 *>          zero.
 142 *> \endverbatim
 143 *>
 144 *> \param[in] RELTOL
 145 *> \verbatim
 146 *>          RELTOL is DOUBLE PRECISION
 147 *>          The minimum relative width of an interval.  When an interval
 148 *>          is narrower than ABSTOL, or than RELTOL times the larger (in
 149 *>          magnitude) endpoint, then it is considered to be
 150 *>          sufficiently small, i.e., converged.  Note: this should
 151 *>          always be at least radix*machine epsilon.
 152 *> \endverbatim
 153 *>
 154 *> \param[in] PIVMIN
 155 *> \verbatim
 156 *>          PIVMIN is DOUBLE PRECISION
 157 *>          The minimum absolute value of a "pivot" in the Sturm
 158 *>          sequence loop.
 159 *>          This must be at least  max |e(j)**2|*safe_min  and at
 160 *>          least safe_min, where safe_min is at least
 161 *>          the smallest number that can divide one without overflow.
 162 *> \endverbatim
 163 *>
 164 *> \param[in] D
 165 *> \verbatim
 166 *>          D is DOUBLE PRECISION array, dimension (N)
 167 *>          The diagonal elements of the tridiagonal matrix T.
 168 *> \endverbatim
 169 *>
 170 *> \param[in] E
 171 *> \verbatim
 172 *>          E is DOUBLE PRECISION array, dimension (N)
 173 *>          The offdiagonal elements of the tridiagonal matrix T in
 174 *>          positions 1 through N-1.  E(N) is arbitrary.
 175 *> \endverbatim
 176 *>
 177 *> \param[in] E2
 178 *> \verbatim
 179 *>          E2 is DOUBLE PRECISION array, dimension (N)
 180 *>          The squares of the offdiagonal elements of the tridiagonal
 181 *>          matrix T.  E2(N) is ignored.
 182 *> \endverbatim
 183 *>
 184 *> \param[in,out] NVAL
 185 *> \verbatim
 186 *>          NVAL is INTEGER array, dimension (MINP)
 187 *>          If IJOB=1 or 2, not referenced.
 188 *>          If IJOB=3, the desired values of N(w).  The elements of NVAL
 189 *>          will be reordered to correspond with the intervals in AB.
 190 *>          Thus, NVAL(j) on output will not, in general be the same as
 191 *>          NVAL(j) on input, but it will correspond with the interval
 192 *>          (AB(j,1),AB(j,2)] on output.
 193 *> \endverbatim
 194 *>
 195 *> \param[in,out] AB
 196 *> \verbatim
 197 *>          AB is DOUBLE PRECISION array, dimension (MMAX,2)
 198 *>          The endpoints of the intervals.  AB(j,1) is  a(j), the left
 199 *>          endpoint of the j-th interval, and AB(j,2) is b(j), the
 200 *>          right endpoint of the j-th interval.  The input intervals
 201 *>          will, in general, be modified, split, and reordered by the
 202 *>          calculation.
 203 *> \endverbatim
 204 *>
 205 *> \param[in,out] C
 206 *> \verbatim
 207 *>          C is DOUBLE PRECISION array, dimension (MMAX)
 208 *>          If IJOB=1, ignored.
 209 *>          If IJOB=2, workspace.
 210 *>          If IJOB=3, then on input C(j) should be initialized to the
 211 *>          first search point in the binary search.
 212 *> \endverbatim
 213 *>
 214 *> \param[out] MOUT
 215 *> \verbatim
 216 *>          MOUT is INTEGER
 217 *>          If IJOB=1, the number of eigenvalues in the intervals.
 218 *>          If IJOB=2 or 3, the number of intervals output.
 219 *>          If IJOB=3, MOUT will equal MINP.
 220 *> \endverbatim
 221 *>
 222 *> \param[in,out] NAB
 223 *> \verbatim
 224 *>          NAB is INTEGER array, dimension (MMAX,2)
 225 *>          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
 226 *>          If IJOB=2, then on input, NAB(i,j) should be set.  It must
 227 *>             satisfy the condition:
 228 *>             N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),
 229 *>             which means that in interval i only eigenvalues
 230 *>             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
 231 *>             NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with
 232 *>             IJOB=1.
 233 *>             On output, NAB(i,j) will contain
 234 *>             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
 235 *>             the input interval that the output interval
 236 *>             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
 237 *>             the input values of NAB(k,1) and NAB(k,2).
 238 *>          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
 239 *>             unless N(w) > NVAL(i) for all search points  w , in which
 240 *>             case NAB(i,1) will not be modified, i.e., the output
 241 *>             value will be the same as the input value (modulo
 242 *>             reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)
 243 *>             for all search points  w , in which case NAB(i,2) will
 244 *>             not be modified.  Normally, NAB should be set to some
 245 *>             distinctive value(s) before DLAEBZ is called.
 246 *> \endverbatim
 247 *>
 248 *> \param[out] WORK
 249 *> \verbatim
 250 *>          WORK is DOUBLE PRECISION array, dimension (MMAX)
 251 *>          Workspace.
 252 *> \endverbatim
 253 *>
 254 *> \param[out] IWORK
 255 *> \verbatim
 256 *>          IWORK is INTEGER array, dimension (MMAX)
 257 *>          Workspace.
 258 *> \endverbatim
 259 *>
 260 *> \param[out] INFO
 261 *> \verbatim
 262 *>          INFO is INTEGER
 263 *>          = 0:       All intervals converged.
 264 *>          = 1--MMAX: The last INFO intervals did not converge.
 265 *>          = MMAX+1:  More than MMAX intervals were generated.
 266 *> \endverbatim
 267 *
 268 *  Authors:
 269 *  ========
 270 *
 271 *> \author Univ. of Tennessee
 272 *> \author Univ. of California Berkeley
 273 *> \author Univ. of Colorado Denver
 274 *> \author NAG Ltd.
 275 *
 276 *> \date September 2012
 277 *
 278 *> \ingroup auxOTHERauxiliary
 279 *
 280 *> \par Further Details:
 281 *  =====================
 282 *>
 283 *> \verbatim
 284 *>
 285 *>      This routine is intended to be called only by other LAPACK
 286 *>  routines, thus the interface is less user-friendly.  It is intended
 287 *>  for two purposes:
 288 *>
 289 *>  (a) finding eigenvalues.  In this case, DLAEBZ should have one or
 290 *>      more initial intervals set up in AB, and DLAEBZ should be called
 291 *>      with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
 292 *>      Intervals with no eigenvalues would usually be thrown out at
 293 *>      this point.  Also, if not all the eigenvalues in an interval i
 294 *>      are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
 295 *>      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
 296 *>      eigenvalue.  DLAEBZ is then called with IJOB=2 and MMAX
 297 *>      no smaller than the value of MOUT returned by the call with
 298 *>      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
 299 *>      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
 300 *>      tolerance specified by ABSTOL and RELTOL.
 301 *>
 302 *>  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
 303 *>      In this case, start with a Gershgorin interval  (a,b).  Set up
 304 *>      AB to contain 2 search intervals, both initially (a,b).  One
 305 *>      NVAL element should contain  f-1  and the other should contain  l
 306 *>      , while C should contain a and b, resp.  NAB(i,1) should be -1
 307 *>      and NAB(i,2) should be N+1, to flag an error if the desired
 308 *>      interval does not lie in (a,b).  DLAEBZ is then called with
 309 *>      IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --
 310 *>      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
 311 *>      if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r
 312 *>      >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
 313 *>      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and
 314 *>      w(l-r)=...=w(l+k) are handled similarly.
 315 *> \endverbatim
 316 *>
 317 *  =====================================================================
 318       SUBROUTINE DLAEBZ( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
 319      $                   RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
 320      $                   NAB, WORK, IWORK, INFO )
 321 *
 322 *  -- LAPACK auxiliary routine (version 3.4.2) --
 323 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 324 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 325 *     September 2012
 326 *
 327 *     .. Scalar Arguments ..
 328       INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
 329       DOUBLE PRECISION   ABSTOL, PIVMIN, RELTOL
 330 *     ..
 331 *     .. Array Arguments ..
 332       INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
 333       DOUBLE PRECISION   AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
 334      $                   WORK( * )
 335 *     ..
 336 *
 337 *  =====================================================================
 338 *
 339 *     .. Parameters ..
 340       DOUBLE PRECISION   ZERO, TWO, HALF
 341       PARAMETER          ( ZERO = 0.0D0, TWO = 2.0D0,
 342      $                   HALF = 1.0D0 / TWO )
 343 *     ..
 344 *     .. Local Scalars ..
 345       INTEGER            ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
 346      $                   KLNEW
 347       DOUBLE PRECISION   TMP1, TMP2
 348 *     ..
 349 *     .. Intrinsic Functions ..
 350       INTRINSIC          ABS, MAX, MIN
 351 *     ..
 352 *     .. Executable Statements ..
 353 *
 354 *     Check for Errors
 355 *
 356       INFO = 0
 357       IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
 358          INFO = -1
 359          RETURN
 360       END IF
 361 *
 362 *     Initialize NAB
 363 *
 364       IF( IJOB.EQ.1 ) THEN
 365 *
 366 *        Compute the number of eigenvalues in the initial intervals.
 367 *
 368          MOUT = 0
 369          DO 30 JI = 1, MINP
 370             DO 20 JP = 1, 2
 371                TMP1 = D( 1 ) - AB( JI, JP )
 372                IF( ABS( TMP1 ).LT.PIVMIN )
 373      $            TMP1 = -PIVMIN
 374                NAB( JI, JP ) = 0
 375                IF( TMP1.LE.ZERO )
 376      $            NAB( JI, JP ) = 1
 377 *
 378                DO 10 J = 2, N
 379                   TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
 380                   IF( ABS( TMP1 ).LT.PIVMIN )
 381      $               TMP1 = -PIVMIN
 382                   IF( TMP1.LE.ZERO )
 383      $               NAB( JI, JP ) = NAB( JI, JP ) + 1
 384    10          CONTINUE
 385    20       CONTINUE
 386             MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
 387    30    CONTINUE
 388          RETURN
 389       END IF
 390 *
 391 *     Initialize for loop
 392 *
 393 *     KF and KL have the following meaning:
 394 *        Intervals 1,...,KF-1 have converged.
 395 *        Intervals KF,...,KL  still need to be refined.
 396 *
 397       KF = 1
 398       KL = MINP
 399 *
 400 *     If IJOB=2, initialize C.
 401 *     If IJOB=3, use the user-supplied starting point.
 402 *
 403       IF( IJOB.EQ.2 ) THEN
 404          DO 40 JI = 1, MINP
 405             C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
 406    40    CONTINUE
 407       END IF
 408 *
 409 *     Iteration loop
 410 *
 411       DO 130 JIT = 1, NITMAX
 412 *
 413 *        Loop over intervals
 414 *
 415          IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
 416 *
 417 *           Begin of Parallel Version of the loop
 418 *
 419             DO 60 JI = KF, KL
 420 *
 421 *              Compute N(c), the number of eigenvalues less than c
 422 *
 423                WORK( JI ) = D( 1 ) - C( JI )
 424                IWORK( JI ) = 0
 425                IF( WORK( JI ).LE.PIVMIN ) THEN
 426                   IWORK( JI ) = 1
 427                   WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
 428                END IF
 429 *
 430                DO 50 J = 2, N
 431                   WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
 432                   IF( WORK( JI ).LE.PIVMIN ) THEN
 433                      IWORK( JI ) = IWORK( JI ) + 1
 434                      WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
 435                   END IF
 436    50          CONTINUE
 437    60       CONTINUE
 438 *
 439             IF( IJOB.LE.2 ) THEN
 440 *
 441 *              IJOB=2: Choose all intervals containing eigenvalues.
 442 *
 443                KLNEW = KL
 444                DO 70 JI = KF, KL
 445 *
 446 *                 Insure that N(w) is monotone
 447 *
 448                   IWORK( JI ) = MIN( NAB( JI, 2 ),
 449      $                          MAX( NAB( JI, 1 ), IWORK( JI ) ) )
 450 *
 451 *                 Update the Queue -- add intervals if both halves
 452 *                 contain eigenvalues.
 453 *
 454                   IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
 455 *
 456 *                    No eigenvalue in the upper interval:
 457 *                    just use the lower interval.
 458 *
 459                      AB( JI, 2 ) = C( JI )
 460 *
 461                   ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
 462 *
 463 *                    No eigenvalue in the lower interval:
 464 *                    just use the upper interval.
 465 *
 466                      AB( JI, 1 ) = C( JI )
 467                   ELSE
 468                      KLNEW = KLNEW + 1
 469                      IF( KLNEW.LE.MMAX ) THEN
 470 *
 471 *                       Eigenvalue in both intervals -- add upper to
 472 *                       queue.
 473 *
 474                         AB( KLNEW, 2 ) = AB( JI, 2 )
 475                         NAB( KLNEW, 2 ) = NAB( JI, 2 )
 476                         AB( KLNEW, 1 ) = C( JI )
 477                         NAB( KLNEW, 1 ) = IWORK( JI )
 478                         AB( JI, 2 ) = C( JI )
 479                         NAB( JI, 2 ) = IWORK( JI )
 480                      ELSE
 481                         INFO = MMAX + 1
 482                      END IF
 483                   END IF
 484    70          CONTINUE
 485                IF( INFO.NE.0 )
 486      $            RETURN
 487                KL = KLNEW
 488             ELSE
 489 *
 490 *              IJOB=3: Binary search.  Keep only the interval containing
 491 *                      w   s.t. N(w) = NVAL
 492 *
 493                DO 80 JI = KF, KL
 494                   IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
 495                      AB( JI, 1 ) = C( JI )
 496                      NAB( JI, 1 ) = IWORK( JI )
 497                   END IF
 498                   IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
 499                      AB( JI, 2 ) = C( JI )
 500                      NAB( JI, 2 ) = IWORK( JI )
 501                   END IF
 502    80          CONTINUE
 503             END IF
 504 *
 505          ELSE
 506 *
 507 *           End of Parallel Version of the loop
 508 *
 509 *           Begin of Serial Version of the loop
 510 *
 511             KLNEW = KL
 512             DO 100 JI = KF, KL
 513 *
 514 *              Compute N(w), the number of eigenvalues less than w
 515 *
 516                TMP1 = C( JI )
 517                TMP2 = D( 1 ) - TMP1
 518                ITMP1 = 0
 519                IF( TMP2.LE.PIVMIN ) THEN
 520                   ITMP1 = 1
 521                   TMP2 = MIN( TMP2, -PIVMIN )
 522                END IF
 523 *
 524                DO 90 J = 2, N
 525                   TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
 526                   IF( TMP2.LE.PIVMIN ) THEN
 527                      ITMP1 = ITMP1 + 1
 528                      TMP2 = MIN( TMP2, -PIVMIN )
 529                   END IF
 530    90          CONTINUE
 531 *
 532                IF( IJOB.LE.2 ) THEN
 533 *
 534 *                 IJOB=2: Choose all intervals containing eigenvalues.
 535 *
 536 *                 Insure that N(w) is monotone
 537 *
 538                   ITMP1 = MIN( NAB( JI, 2 ),
 539      $                    MAX( NAB( JI, 1 ), ITMP1 ) )
 540 *
 541 *                 Update the Queue -- add intervals if both halves
 542 *                 contain eigenvalues.
 543 *
 544                   IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
 545 *
 546 *                    No eigenvalue in the upper interval:
 547 *                    just use the lower interval.
 548 *
 549                      AB( JI, 2 ) = TMP1
 550 *
 551                   ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
 552 *
 553 *                    No eigenvalue in the lower interval:
 554 *                    just use the upper interval.
 555 *
 556                      AB( JI, 1 ) = TMP1
 557                   ELSE IF( KLNEW.LT.MMAX ) THEN
 558 *
 559 *                    Eigenvalue in both intervals -- add upper to queue.
 560 *
 561                      KLNEW = KLNEW + 1
 562                      AB( KLNEW, 2 ) = AB( JI, 2 )
 563                      NAB( KLNEW, 2 ) = NAB( JI, 2 )
 564                      AB( KLNEW, 1 ) = TMP1
 565                      NAB( KLNEW, 1 ) = ITMP1
 566                      AB( JI, 2 ) = TMP1
 567                      NAB( JI, 2 ) = ITMP1
 568                   ELSE
 569                      INFO = MMAX + 1
 570                      RETURN
 571                   END IF
 572                ELSE
 573 *
 574 *                 IJOB=3: Binary search.  Keep only the interval
 575 *                         containing  w  s.t. N(w) = NVAL
 576 *
 577                   IF( ITMP1.LE.NVAL( JI ) ) THEN
 578                      AB( JI, 1 ) = TMP1
 579                      NAB( JI, 1 ) = ITMP1
 580                   END IF
 581                   IF( ITMP1.GE.NVAL( JI ) ) THEN
 582                      AB( JI, 2 ) = TMP1
 583                      NAB( JI, 2 ) = ITMP1
 584                   END IF
 585                END IF
 586   100       CONTINUE
 587             KL = KLNEW
 588 *
 589          END IF
 590 *
 591 *        Check for convergence
 592 *
 593          KFNEW = KF
 594          DO 110 JI = KF, KL
 595             TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
 596             TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
 597             IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
 598      $          NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
 599 *
 600 *              Converged -- Swap with position KFNEW,
 601 *                           then increment KFNEW
 602 *
 603                IF( JI.GT.KFNEW ) THEN
 604                   TMP1 = AB( JI, 1 )
 605                   TMP2 = AB( JI, 2 )
 606                   ITMP1 = NAB( JI, 1 )
 607                   ITMP2 = NAB( JI, 2 )
 608                   AB( JI, 1 ) = AB( KFNEW, 1 )
 609                   AB( JI, 2 ) = AB( KFNEW, 2 )
 610                   NAB( JI, 1 ) = NAB( KFNEW, 1 )
 611                   NAB( JI, 2 ) = NAB( KFNEW, 2 )
 612                   AB( KFNEW, 1 ) = TMP1
 613                   AB( KFNEW, 2 ) = TMP2
 614                   NAB( KFNEW, 1 ) = ITMP1
 615                   NAB( KFNEW, 2 ) = ITMP2
 616                   IF( IJOB.EQ.3 ) THEN
 617                      ITMP1 = NVAL( JI )
 618                      NVAL( JI ) = NVAL( KFNEW )
 619                      NVAL( KFNEW ) = ITMP1
 620                   END IF
 621                END IF
 622                KFNEW = KFNEW + 1
 623             END IF
 624   110    CONTINUE
 625          KF = KFNEW
 626 *
 627 *        Choose Midpoints
 628 *
 629          DO 120 JI = KF, KL
 630             C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
 631   120    CONTINUE
 632 *
 633 *        If no more intervals to refine, quit.
 634 *
 635          IF( KF.GT.KL )
 636      $      GO TO 140
 637   130 CONTINUE
 638 *
 639 *     Converged
 640 *
 641   140 CONTINUE
 642       INFO = MAX( KL+1-KF, 0 )
 643       MOUT = KL
 644 *
 645       RETURN
 646 *
 647 *     End of DLAEBZ
 648 *
 649       END