SRC/sstebz.f

   1 *> \brief \b SSTEBZ
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSTEBZ + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstebz.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstebz.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstebz.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
  22 *                          M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
  23 *                          INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          ORDER, RANGE
  27 *       INTEGER            IL, INFO, IU, M, N, NSPLIT
  28 *       REAL               ABSTOL, VL, VU
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
  32 *       REAL               D( * ), E( * ), W( * ), WORK( * )
  33 *       ..
  34 *
  35 *
  36 *> \par Purpose:
  37 *  =============
  38 *>
  39 *> \verbatim
  40 *>
  41 *> SSTEBZ computes the eigenvalues of a symmetric tridiagonal
  42 *> matrix T.  The user may ask for all eigenvalues, all eigenvalues
  43 *> in the half-open interval (VL, VU], or the IL-th through IU-th
  44 *> eigenvalues.
  45 *>
  46 *> To avoid overflow, the matrix must be scaled so that its
  47 *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
  48 *> accuracy, it should not be much smaller than that.
  49 *>
  50 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
  51 *> Matrix", Report CS41, Computer Science Dept., Stanford
  52 *> University, July 21, 1966.
  53 *> \endverbatim
  54 *
  55 *  Arguments:
  56 *  ==========
  57 *
  58 *> \param[in] RANGE
  59 *> \verbatim
  60 *>          RANGE is CHARACTER*1
  61 *>          = 'A': ("All")   all eigenvalues will be found.
  62 *>          = 'V': ("Value") all eigenvalues in the half-open interval
  63 *>                           (VL, VU] will be found.
  64 *>          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
  65 *>                           entire matrix) will be found.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] ORDER
  69 *> \verbatim
  70 *>          ORDER is CHARACTER*1
  71 *>          = 'B': ("By Block") the eigenvalues will be grouped by
  72 *>                              split-off block (see IBLOCK, ISPLIT) and
  73 *>                              ordered from smallest to largest within
  74 *>                              the block.
  75 *>          = 'E': ("Entire matrix")
  76 *>                              the eigenvalues for the entire matrix
  77 *>                              will be ordered from smallest to
  78 *>                              largest.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] N
  82 *> \verbatim
  83 *>          N is INTEGER
  84 *>          The order of the tridiagonal matrix T.  N >= 0.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] VL
  88 *> \verbatim
  89 *>          VL is REAL
  90 *>
  91 *>          If RANGE='V', the lower bound of the interval to
  92 *>          be searched for eigenvalues.  Eigenvalues less than or equal
  93 *>          to VL, or greater than VU, will not be returned.  VL < VU.
  94 *>          Not referenced if RANGE = 'A' or 'I'.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] VU
  98 *> \verbatim
  99 *>          VU is REAL
 100 *>
 101 *>          If RANGE='V', the upper bound of the interval to
 102 *>          be searched for eigenvalues.  Eigenvalues less than or equal
 103 *>          to VL, or greater than VU, will not be returned.  VL < VU.
 104 *>          Not referenced if RANGE = 'A' or 'I'.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] IL
 108 *> \verbatim
 109 *>          IL is INTEGER
 110 *>
 111 *>          If RANGE='I', the index of the
 112 *>          smallest eigenvalue to be returned.
 113 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
 114 *>          Not referenced if RANGE = 'A' or 'V'.
 115 *> \endverbatim
 116 *>
 117 *> \param[in] IU
 118 *> \verbatim
 119 *>          IU is INTEGER
 120 *>
 121 *>          If RANGE='I', the index of the
 122 *>          largest eigenvalue to be returned.
 123 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
 124 *>          Not referenced if RANGE = 'A' or 'V'.
 125 *> \endverbatim
 126 *>
 127 *> \param[in] ABSTOL
 128 *> \verbatim
 129 *>          ABSTOL is REAL
 130 *>          The absolute tolerance for the eigenvalues.  An eigenvalue
 131 *>          (or cluster) is considered to be located if it has been
 132 *>          determined to lie in an interval whose width is ABSTOL or
 133 *>          less.  If ABSTOL is less than or equal to zero, then ULP*|T|
 134 *>          will be used, where |T| means the 1-norm of T.
 135 *>
 136 *>          Eigenvalues will be computed most accurately when ABSTOL is
 137 *>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
 138 *> \endverbatim
 139 *>
 140 *> \param[in] D
 141 *> \verbatim
 142 *>          D is REAL array, dimension (N)
 143 *>          The n diagonal elements of the tridiagonal matrix T.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] E
 147 *> \verbatim
 148 *>          E is REAL array, dimension (N-1)
 149 *>          The (n-1) off-diagonal elements of the tridiagonal matrix T.
 150 *> \endverbatim
 151 *>
 152 *> \param[out] M
 153 *> \verbatim
 154 *>          M is INTEGER
 155 *>          The actual number of eigenvalues found. 0 <= M <= N.
 156 *>          (See also the description of INFO=2,3.)
 157 *> \endverbatim
 158 *>
 159 *> \param[out] NSPLIT
 160 *> \verbatim
 161 *>          NSPLIT is INTEGER
 162 *>          The number of diagonal blocks in the matrix T.
 163 *>          1 <= NSPLIT <= N.
 164 *> \endverbatim
 165 *>
 166 *> \param[out] W
 167 *> \verbatim
 168 *>          W is REAL array, dimension (N)
 169 *>          On exit, the first M elements of W will contain the
 170 *>          eigenvalues.  (SSTEBZ may use the remaining N-M elements as
 171 *>          workspace.)
 172 *> \endverbatim
 173 *>
 174 *> \param[out] IBLOCK
 175 *> \verbatim
 176 *>          IBLOCK is INTEGER array, dimension (N)
 177 *>          At each row/column j where E(j) is zero or small, the
 178 *>          matrix T is considered to split into a block diagonal
 179 *>          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which
 180 *>          block (from 1 to the number of blocks) the eigenvalue W(i)
 181 *>          belongs.  (SSTEBZ may use the remaining N-M elements as
 182 *>          workspace.)
 183 *> \endverbatim
 184 *>
 185 *> \param[out] ISPLIT
 186 *> \verbatim
 187 *>          ISPLIT is INTEGER array, dimension (N)
 188 *>          The splitting points, at which T breaks up into submatrices.
 189 *>          The first submatrix consists of rows/columns 1 to ISPLIT(1),
 190 *>          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
 191 *>          etc., and the NSPLIT-th consists of rows/columns
 192 *>          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
 193 *>          (Only the first NSPLIT elements will actually be used, but
 194 *>          since the user cannot know a priori what value NSPLIT will
 195 *>          have, N words must be reserved for ISPLIT.)
 196 *> \endverbatim
 197 *>
 198 *> \param[out] WORK
 199 *> \verbatim
 200 *>          WORK is REAL array, dimension (4*N)
 201 *> \endverbatim
 202 *>
 203 *> \param[out] IWORK
 204 *> \verbatim
 205 *>          IWORK is INTEGER array, dimension (3*N)
 206 *> \endverbatim
 207 *>
 208 *> \param[out] INFO
 209 *> \verbatim
 210 *>          INFO is INTEGER
 211 *>          = 0:  successful exit
 212 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 213 *>          > 0:  some or all of the eigenvalues failed to converge or
 214 *>                were not computed:
 215 *>                =1 or 3: Bisection failed to converge for some
 216 *>                        eigenvalues; these eigenvalues are flagged by a
 217 *>                        negative block number.  The effect is that the
 218 *>                        eigenvalues may not be as accurate as the
 219 *>                        absolute and relative tolerances.  This is
 220 *>                        generally caused by unexpectedly inaccurate
 221 *>                        arithmetic.
 222 *>                =2 or 3: RANGE='I' only: Not all of the eigenvalues
 223 *>                        IL:IU were found.
 224 *>                        Effect: M < IU+1-IL
 225 *>                        Cause:  non-monotonic arithmetic, causing the
 226 *>                                Sturm sequence to be non-monotonic.
 227 *>                        Cure:   recalculate, using RANGE='A', and pick
 228 *>                                out eigenvalues IL:IU.  In some cases,
 229 *>                                increasing the PARAMETER "FUDGE" may
 230 *>                                make things work.
 231 *>                = 4:    RANGE='I', and the Gershgorin interval
 232 *>                        initially used was too small.  No eigenvalues
 233 *>                        were computed.
 234 *>                        Probable cause: your machine has sloppy
 235 *>                                        floating-point arithmetic.
 236 *>                        Cure: Increase the PARAMETER "FUDGE",
 237 *>                              recompile, and try again.
 238 *> \endverbatim
 239 *
 240 *> \par Internal Parameters:
 241 *  =========================
 242 *>
 243 *> \verbatim
 244 *>  RELFAC  REAL, default = 2.0e0
 245 *>          The relative tolerance.  An interval (a,b] lies within
 246 *>          "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|),
 247 *>          where "ulp" is the machine precision (distance from 1 to
 248 *>          the next larger floating point number.)
 249 *>
 250 *>  FUDGE   REAL, default = 2
 251 *>          A "fudge factor" to widen the Gershgorin intervals.  Ideally,
 252 *>          a value of 1 should work, but on machines with sloppy
 253 *>          arithmetic, this needs to be larger.  The default for
 254 *>          publicly released versions should be large enough to handle
 255 *>          the worst machine around.  Note that this has no effect
 256 *>          on accuracy of the solution.
 257 *> \endverbatim
 258 *
 259 *  Authors:
 260 *  ========
 261 *
 262 *> \author Univ. of Tennessee
 263 *> \author Univ. of California Berkeley
 264 *> \author Univ. of Colorado Denver
 265 *> \author NAG Ltd.
 266 *
 267 *> \date June 2016
 268 *
 269 *> \ingroup auxOTHERcomputational
 270 *
 271 *  =====================================================================
 272       SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E,
 273      $                   M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK,
 274      $                   INFO )
 275 *
 276 *  -- LAPACK computational routine (version 3.6.1) --
 277 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 278 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 279 *     June 2016
 280 *
 281 *     .. Scalar Arguments ..
 282       CHARACTER          ORDER, RANGE
 283       INTEGER            IL, INFO, IU, M, N, NSPLIT
 284       REAL               ABSTOL, VL, VU
 285 *     ..
 286 *     .. Array Arguments ..
 287       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * )
 288       REAL               D( * ), E( * ), W( * ), WORK( * )
 289 *     ..
 290 *
 291 *  =====================================================================
 292 *
 293 *     .. Parameters ..
 294       REAL               ZERO, ONE, TWO, HALF
 295       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
 296      $                   HALF = 1.0E0 / TWO )
 297       REAL               FUDGE, RELFAC
 298       PARAMETER          ( FUDGE = 2.1E0, RELFAC = 2.0E0 )
 299 *     ..
 300 *     .. Local Scalars ..
 301       LOGICAL            NCNVRG, TOOFEW
 302       INTEGER            IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
 303      $                   IM, IN, IOFF, IORDER, IOUT, IRANGE, ITMAX,
 304      $                   ITMP1, IW, IWOFF, J, JB, JDISC, JE, NB, NWL,
 305      $                   NWU
 306       REAL               ATOLI, BNORM, GL, GU, PIVMIN, RTOLI, SAFEMN,
 307      $                   TMP1, TMP2, TNORM, ULP, WKILL, WL, WLU, WU, WUL
 308 *     ..
 309 *     .. Local Arrays ..
 310       INTEGER            IDUMMA( 1 )
 311 *     ..
 312 *     .. External Functions ..
 313       LOGICAL            LSAME
 314       INTEGER            ILAENV
 315       REAL               SLAMCH
 316       EXTERNAL           LSAME, ILAENV, SLAMCH
 317 *     ..
 318 *     .. External Subroutines ..
 319       EXTERNAL           SLAEBZ, XERBLA
 320 *     ..
 321 *     .. Intrinsic Functions ..
 322       INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT
 323 *     ..
 324 *     .. Executable Statements ..
 325 *
 326       INFO = 0
 327 *
 328 *     Decode RANGE
 329 *
 330       IF( LSAME( RANGE, 'A' ) ) THEN
 331          IRANGE = 1
 332       ELSE IF( LSAME( RANGE, 'V' ) ) THEN
 333          IRANGE = 2
 334       ELSE IF( LSAME( RANGE, 'I' ) ) THEN
 335          IRANGE = 3
 336       ELSE
 337          IRANGE = 0
 338       END IF
 339 *
 340 *     Decode ORDER
 341 *
 342       IF( LSAME( ORDER, 'B' ) ) THEN
 343          IORDER = 2
 344       ELSE IF( LSAME( ORDER, 'E' ) ) THEN
 345          IORDER = 1
 346       ELSE
 347          IORDER = 0
 348       END IF
 349 *
 350 *     Check for Errors
 351 *
 352       IF( IRANGE.LE.0 ) THEN
 353          INFO = -1
 354       ELSE IF( IORDER.LE.0 ) THEN
 355          INFO = -2
 356       ELSE IF( N.LT.0 ) THEN
 357          INFO = -3
 358       ELSE IF( IRANGE.EQ.2 ) THEN
 359          IF( VL.GE.VU ) INFO = -5
 360       ELSE IF( IRANGE.EQ.3 .AND. ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) )
 361      $          THEN
 362          INFO = -6
 363       ELSE IF( IRANGE.EQ.3 .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) )
 364      $          THEN
 365          INFO = -7
 366       END IF
 367 *
 368       IF( INFO.NE.0 ) THEN
 369          CALL XERBLA( 'SSTEBZ', -INFO )
 370          RETURN
 371       END IF
 372 *
 373 *     Initialize error flags
 374 *
 375       INFO = 0
 376       NCNVRG = .FALSE.
 377       TOOFEW = .FALSE.
 378 *
 379 *     Quick return if possible
 380 *
 381       M = 0
 382       IF( N.EQ.0 )
 383      $   RETURN
 384 *
 385 *     Simplifications:
 386 *
 387       IF( IRANGE.EQ.3 .AND. IL.EQ.1 .AND. IU.EQ.N )
 388      $   IRANGE = 1
 389 *
 390 *     Get machine constants
 391 *     NB is the minimum vector length for vector bisection, or 0
 392 *     if only scalar is to be done.
 393 *
 394       SAFEMN = SLAMCH( 'S' )
 395       ULP = SLAMCH( 'P' )
 396       RTOLI = ULP*RELFAC
 397       NB = ILAENV( 1, 'SSTEBZ', ' ', N, -1, -1, -1 )
 398       IF( NB.LE.1 )
 399      $   NB = 0
 400 *
 401 *     Special Case when N=1
 402 *
 403       IF( N.EQ.1 ) THEN
 404          NSPLIT = 1
 405          ISPLIT( 1 ) = 1
 406          IF( IRANGE.EQ.2 .AND. ( VL.GE.D( 1 ) .OR. VU.LT.D( 1 ) ) ) THEN
 407             M = 0
 408          ELSE
 409             W( 1 ) = D( 1 )
 410             IBLOCK( 1 ) = 1
 411             M = 1
 412          END IF
 413          RETURN
 414       END IF
 415 *
 416 *     Compute Splitting Points
 417 *
 418       NSPLIT = 1
 419       WORK( N ) = ZERO
 420       PIVMIN = ONE
 421 *
 422       DO 10 J = 2, N
 423          TMP1 = E( J-1 )**2
 424          IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
 425             ISPLIT( NSPLIT ) = J - 1
 426             NSPLIT = NSPLIT + 1
 427             WORK( J-1 ) = ZERO
 428          ELSE
 429             WORK( J-1 ) = TMP1
 430             PIVMIN = MAX( PIVMIN, TMP1 )
 431          END IF
 432    10 CONTINUE
 433       ISPLIT( NSPLIT ) = N
 434       PIVMIN = PIVMIN*SAFEMN
 435 *
 436 *     Compute Interval and ATOLI
 437 *
 438       IF( IRANGE.EQ.3 ) THEN
 439 *
 440 *        RANGE='I': Compute the interval containing eigenvalues
 441 *                   IL through IU.
 442 *
 443 *        Compute Gershgorin interval for entire (split) matrix
 444 *        and use it as the initial interval
 445 *
 446          GU = D( 1 )
 447          GL = D( 1 )
 448          TMP1 = ZERO
 449 *
 450          DO 20 J = 1, N - 1
 451             TMP2 = SQRT( WORK( J ) )
 452             GU = MAX( GU, D( J )+TMP1+TMP2 )
 453             GL = MIN( GL, D( J )-TMP1-TMP2 )
 454             TMP1 = TMP2
 455    20    CONTINUE
 456 *
 457          GU = MAX( GU, D( N )+TMP1 )
 458          GL = MIN( GL, D( N )-TMP1 )
 459          TNORM = MAX( ABS( GL ), ABS( GU ) )
 460          GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
 461          GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
 462 *
 463 *        Compute Iteration parameters
 464 *
 465          ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
 466      $           LOG( TWO ) ) + 2
 467          IF( ABSTOL.LE.ZERO ) THEN
 468             ATOLI = ULP*TNORM
 469          ELSE
 470             ATOLI = ABSTOL
 471          END IF
 472 *
 473          WORK( N+1 ) = GL
 474          WORK( N+2 ) = GL
 475          WORK( N+3 ) = GU
 476          WORK( N+4 ) = GU
 477          WORK( N+5 ) = GL
 478          WORK( N+6 ) = GU
 479          IWORK( 1 ) = -1
 480          IWORK( 2 ) = -1
 481          IWORK( 3 ) = N + 1
 482          IWORK( 4 ) = N + 1
 483          IWORK( 5 ) = IL - 1
 484          IWORK( 6 ) = IU
 485 *
 486          CALL SLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
 487      $                WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
 488      $                IWORK, W, IBLOCK, IINFO )
 489 *
 490          IF( IWORK( 6 ).EQ.IU ) THEN
 491             WL = WORK( N+1 )
 492             WLU = WORK( N+3 )
 493             NWL = IWORK( 1 )
 494             WU = WORK( N+4 )
 495             WUL = WORK( N+2 )
 496             NWU = IWORK( 4 )
 497          ELSE
 498             WL = WORK( N+2 )
 499             WLU = WORK( N+4 )
 500             NWL = IWORK( 2 )
 501             WU = WORK( N+3 )
 502             WUL = WORK( N+1 )
 503             NWU = IWORK( 3 )
 504          END IF
 505 *
 506          IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN
 507             INFO = 4
 508             RETURN
 509          END IF
 510       ELSE
 511 *
 512 *        RANGE='A' or 'V' -- Set ATOLI
 513 *
 514          TNORM = MAX( ABS( D( 1 ) )+ABS( E( 1 ) ),
 515      $           ABS( D( N ) )+ABS( E( N-1 ) ) )
 516 *
 517          DO 30 J = 2, N - 1
 518             TNORM = MAX( TNORM, ABS( D( J ) )+ABS( E( J-1 ) )+
 519      $              ABS( E( J ) ) )
 520    30    CONTINUE
 521 *
 522          IF( ABSTOL.LE.ZERO ) THEN
 523             ATOLI = ULP*TNORM
 524          ELSE
 525             ATOLI = ABSTOL
 526          END IF
 527 *
 528          IF( IRANGE.EQ.2 ) THEN
 529             WL = VL
 530             WU = VU
 531          ELSE
 532             WL = ZERO
 533             WU = ZERO
 534          END IF
 535       END IF
 536 *
 537 *     Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
 538 *     NWL accumulates the number of eigenvalues .le. WL,
 539 *     NWU accumulates the number of eigenvalues .le. WU
 540 *
 541       M = 0
 542       IEND = 0
 543       INFO = 0
 544       NWL = 0
 545       NWU = 0
 546 *
 547       DO 70 JB = 1, NSPLIT
 548          IOFF = IEND
 549          IBEGIN = IOFF + 1
 550          IEND = ISPLIT( JB )
 551          IN = IEND - IOFF
 552 *
 553          IF( IN.EQ.1 ) THEN
 554 *
 555 *           Special Case -- IN=1
 556 *
 557             IF( IRANGE.EQ.1 .OR. WL.GE.D( IBEGIN )-PIVMIN )
 558      $         NWL = NWL + 1
 559             IF( IRANGE.EQ.1 .OR. WU.GE.D( IBEGIN )-PIVMIN )
 560      $         NWU = NWU + 1
 561             IF( IRANGE.EQ.1 .OR. ( WL.LT.D( IBEGIN )-PIVMIN .AND. WU.GE.
 562      $          D( IBEGIN )-PIVMIN ) ) THEN
 563                M = M + 1
 564                W( M ) = D( IBEGIN )
 565                IBLOCK( M ) = JB
 566             END IF
 567          ELSE
 568 *
 569 *           General Case -- IN > 1
 570 *
 571 *           Compute Gershgorin Interval
 572 *           and use it as the initial interval
 573 *
 574             GU = D( IBEGIN )
 575             GL = D( IBEGIN )
 576             TMP1 = ZERO
 577 *
 578             DO 40 J = IBEGIN, IEND - 1
 579                TMP2 = ABS( E( J ) )
 580                GU = MAX( GU, D( J )+TMP1+TMP2 )
 581                GL = MIN( GL, D( J )-TMP1-TMP2 )
 582                TMP1 = TMP2
 583    40       CONTINUE
 584 *
 585             GU = MAX( GU, D( IEND )+TMP1 )
 586             GL = MIN( GL, D( IEND )-TMP1 )
 587             BNORM = MAX( ABS( GL ), ABS( GU ) )
 588             GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN
 589             GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN
 590 *
 591 *           Compute ATOLI for the current submatrix
 592 *
 593             IF( ABSTOL.LE.ZERO ) THEN
 594                ATOLI = ULP*MAX( ABS( GL ), ABS( GU ) )
 595             ELSE
 596                ATOLI = ABSTOL
 597             END IF
 598 *
 599             IF( IRANGE.GT.1 ) THEN
 600                IF( GU.LT.WL ) THEN
 601                   NWL = NWL + IN
 602                   NWU = NWU + IN
 603                   GO TO 70
 604                END IF
 605                GL = MAX( GL, WL )
 606                GU = MIN( GU, WU )
 607                IF( GL.GE.GU )
 608      $            GO TO 70
 609             END IF
 610 *
 611 *           Set Up Initial Interval
 612 *
 613             WORK( N+1 ) = GL
 614             WORK( N+IN+1 ) = GU
 615             CALL SLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
 616      $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
 617      $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
 618      $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
 619 *
 620             NWL = NWL + IWORK( 1 )
 621             NWU = NWU + IWORK( IN+1 )
 622             IWOFF = M - IWORK( 1 )
 623 *
 624 *           Compute Eigenvalues
 625 *
 626             ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) /
 627      $              LOG( TWO ) ) + 2
 628             CALL SLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
 629      $                   D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
 630      $                   IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
 631      $                   IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
 632 *
 633 *           Copy Eigenvalues Into W and IBLOCK
 634 *           Use -JB for block number for unconverged eigenvalues.
 635 *
 636             DO 60 J = 1, IOUT
 637                TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) )
 638 *
 639 *              Flag non-convergence.
 640 *
 641                IF( J.GT.IOUT-IINFO ) THEN
 642                   NCNVRG = .TRUE.
 643                   IB = -JB
 644                ELSE
 645                   IB = JB
 646                END IF
 647                DO 50 JE = IWORK( J ) + 1 + IWOFF,
 648      $                 IWORK( J+IN ) + IWOFF
 649                   W( JE ) = TMP1
 650                   IBLOCK( JE ) = IB
 651    50          CONTINUE
 652    60       CONTINUE
 653 *
 654             M = M + IM
 655          END IF
 656    70 CONTINUE
 657 *
 658 *     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
 659 *     If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
 660 *
 661       IF( IRANGE.EQ.3 ) THEN
 662          IM = 0
 663          IDISCL = IL - 1 - NWL
 664          IDISCU = NWU - IU
 665 *
 666          IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
 667             DO 80 JE = 1, M
 668                IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN
 669                   IDISCL = IDISCL - 1
 670                ELSE IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN
 671                   IDISCU = IDISCU - 1
 672                ELSE
 673                   IM = IM + 1
 674                   W( IM ) = W( JE )
 675                   IBLOCK( IM ) = IBLOCK( JE )
 676                END IF
 677    80       CONTINUE
 678             M = IM
 679          END IF
 680          IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN
 681 *
 682 *           Code to deal with effects of bad arithmetic:
 683 *           Some low eigenvalues to be discarded are not in (WL,WLU],
 684 *           or high eigenvalues to be discarded are not in (WUL,WU]
 685 *           so just kill off the smallest IDISCL/largest IDISCU
 686 *           eigenvalues, by simply finding the smallest/largest
 687 *           eigenvalue(s).
 688 *
 689 *           (If N(w) is monotone non-decreasing, this should never
 690 *               happen.)
 691 *
 692             IF( IDISCL.GT.0 ) THEN
 693                WKILL = WU
 694                DO 100 JDISC = 1, IDISCL
 695                   IW = 0
 696                   DO 90 JE = 1, M
 697                      IF( IBLOCK( JE ).NE.0 .AND.
 698      $                   ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN
 699                         IW = JE
 700                         WKILL = W( JE )
 701                      END IF
 702    90             CONTINUE
 703                   IBLOCK( IW ) = 0
 704   100          CONTINUE
 705             END IF
 706             IF( IDISCU.GT.0 ) THEN
 707 *
 708                WKILL = WL
 709                DO 120 JDISC = 1, IDISCU
 710                   IW = 0
 711                   DO 110 JE = 1, M
 712                      IF( IBLOCK( JE ).NE.0 .AND.
 713      $                   ( W( JE ).GT.WKILL .OR. IW.EQ.0 ) ) THEN
 714                         IW = JE
 715                         WKILL = W( JE )
 716                      END IF
 717   110             CONTINUE
 718                   IBLOCK( IW ) = 0
 719   120          CONTINUE
 720             END IF
 721             IM = 0
 722             DO 130 JE = 1, M
 723                IF( IBLOCK( JE ).NE.0 ) THEN
 724                   IM = IM + 1
 725                   W( IM ) = W( JE )
 726                   IBLOCK( IM ) = IBLOCK( JE )
 727                END IF
 728   130       CONTINUE
 729             M = IM
 730          END IF
 731          IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN
 732             TOOFEW = .TRUE.
 733          END IF
 734       END IF
 735 *
 736 *     If ORDER='B', do nothing -- the eigenvalues are already sorted
 737 *        by block.
 738 *     If ORDER='E', sort the eigenvalues from smallest to largest
 739 *
 740       IF( IORDER.EQ.1 .AND. NSPLIT.GT.1 ) THEN
 741          DO 150 JE = 1, M - 1
 742             IE = 0
 743             TMP1 = W( JE )
 744             DO 140 J = JE + 1, M
 745                IF( W( J ).LT.TMP1 ) THEN
 746                   IE = J
 747                   TMP1 = W( J )
 748                END IF
 749   140       CONTINUE
 750 *
 751             IF( IE.NE.0 ) THEN
 752                ITMP1 = IBLOCK( IE )
 753                W( IE ) = W( JE )
 754                IBLOCK( IE ) = IBLOCK( JE )
 755                W( JE ) = TMP1
 756                IBLOCK( JE ) = ITMP1
 757             END IF
 758   150    CONTINUE
 759       END IF
 760 *
 761       INFO = 0
 762       IF( NCNVRG )
 763      $   INFO = INFO + 1
 764       IF( TOOFEW )
 765      $   INFO = INFO + 2
 766       RETURN
 767 *
 768 *     End of SSTEBZ
 769 *
 770       END