SRC/slarre.f

   1 *> \brief \b SLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLARRE + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarre.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarre.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarre.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
  22 *                           RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
  23 *                           W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
  24 *                           WORK, IWORK, INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       CHARACTER          RANGE
  28 *       INTEGER            IL, INFO, IU, M, N, NSPLIT
  29 *       REAL               PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
  33 *      $                   INDEXW( * )
  34 *       REAL               D( * ), E( * ), E2( * ), GERS( * ),
  35 *      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> To find the desired eigenvalues of a given real symmetric
  45 *> tridiagonal matrix T, SLARRE sets any "small" off-diagonal
  46 *> elements to zero, and for each unreduced block T_i, it finds
  47 *> (a) a suitable shift at one end of the block's spectrum,
  48 *> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
  49 *> (c) eigenvalues of each L_i D_i L_i^T.
  50 *> The representations and eigenvalues found are then used by
  51 *> SSTEMR to compute the eigenvectors of T.
  52 *> The accuracy varies depending on whether bisection is used to
  53 *> find a few eigenvalues or the dqds algorithm (subroutine SLASQ2) to
  54 *> conpute all and then discard any unwanted one.
  55 *> As an added benefit, SLARRE also outputs the n
  56 *> Gerschgorin intervals for the matrices L_i D_i L_i^T.
  57 *> \endverbatim
  58 *
  59 *  Arguments:
  60 *  ==========
  61 *
  62 *> \param[in] RANGE
  63 *> \verbatim
  64 *>          RANGE is CHARACTER*1
  65 *>          = 'A': ("All")   all eigenvalues will be found.
  66 *>          = 'V': ("Value") all eigenvalues in the half-open interval
  67 *>                           (VL, VU] will be found.
  68 *>          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
  69 *>                           entire matrix) will be found.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] N
  73 *> \verbatim
  74 *>          N is INTEGER
  75 *>          The order of the matrix. N > 0.
  76 *> \endverbatim
  77 *>
  78 *> \param[in,out] VL
  79 *> \verbatim
  80 *>          VL is REAL
  81 *>          If RANGE='V', the lower bound for the eigenvalues.
  82 *>          Eigenvalues less than or equal to VL, or greater than VU,
  83 *>          will not be returned.  VL < VU.
  84 *>          If RANGE='I' or ='A', SLARRE computes bounds on the desired
  85 *>          part of the spectrum.
  86 *> \endverbatim
  87 *>
  88 *> \param[in,out] VU
  89 *> \verbatim
  90 *>          VU is REAL
  91 *>          If RANGE='V', the upper bound for the eigenvalues.
  92 *>          Eigenvalues less than or equal to VL, or greater than VU,
  93 *>          will not be returned.  VL < VU.
  94 *>          If RANGE='I' or ='A', SLARRE computes bounds on the desired
  95 *>          part of the spectrum.
  96 *> \endverbatim
  97 *>
  98 *> \param[in] IL
  99 *> \verbatim
 100 *>          IL is INTEGER
 101 *>          If RANGE='I', the index of the
 102 *>          smallest eigenvalue to be returned.
 103 *>          1 <= IL <= IU <= N.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] IU
 107 *> \verbatim
 108 *>          IU is INTEGER
 109 *>          If RANGE='I', the index of the
 110 *>          largest eigenvalue to be returned.
 111 *>          1 <= IL <= IU <= N.
 112 *> \endverbatim
 113 *>
 114 *> \param[in,out] D
 115 *> \verbatim
 116 *>          D is REAL array, dimension (N)
 117 *>          On entry, the N diagonal elements of the tridiagonal
 118 *>          matrix T.
 119 *>          On exit, the N diagonal elements of the diagonal
 120 *>          matrices D_i.
 121 *> \endverbatim
 122 *>
 123 *> \param[in,out] E
 124 *> \verbatim
 125 *>          E is REAL array, dimension (N)
 126 *>          On entry, the first (N-1) entries contain the subdiagonal
 127 *>          elements of the tridiagonal matrix T; E(N) need not be set.
 128 *>          On exit, E contains the subdiagonal elements of the unit
 129 *>          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
 130 *>          1 <= I <= NSPLIT, contain the base points sigma_i on output.
 131 *> \endverbatim
 132 *>
 133 *> \param[in,out] E2
 134 *> \verbatim
 135 *>          E2 is REAL array, dimension (N)
 136 *>          On entry, the first (N-1) entries contain the SQUARES of the
 137 *>          subdiagonal elements of the tridiagonal matrix T;
 138 *>          E2(N) need not be set.
 139 *>          On exit, the entries E2( ISPLIT( I ) ),
 140 *>          1 <= I <= NSPLIT, have been set to zero
 141 *> \endverbatim
 142 *>
 143 *> \param[in] RTOL1
 144 *> \verbatim
 145 *>          RTOL1 is REAL
 146 *> \endverbatim
 147 *>
 148 *> \param[in] RTOL2
 149 *> \verbatim
 150 *>          RTOL2 is REAL
 151 *>           Parameters for bisection.
 152 *>           An interval [LEFT,RIGHT] has converged if
 153 *>           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
 154 *> \endverbatim
 155 *>
 156 *> \param[in] SPLTOL
 157 *> \verbatim
 158 *>          SPLTOL is REAL
 159 *>          The threshold for splitting.
 160 *> \endverbatim
 161 *>
 162 *> \param[out] NSPLIT
 163 *> \verbatim
 164 *>          NSPLIT is INTEGER
 165 *>          The number of blocks T splits into. 1 <= NSPLIT <= N.
 166 *> \endverbatim
 167 *>
 168 *> \param[out] ISPLIT
 169 *> \verbatim
 170 *>          ISPLIT is INTEGER array, dimension (N)
 171 *>          The splitting points, at which T breaks up into blocks.
 172 *>          The first block consists of rows/columns 1 to ISPLIT(1),
 173 *>          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
 174 *>          etc., and the NSPLIT-th consists of rows/columns
 175 *>          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
 176 *> \endverbatim
 177 *>
 178 *> \param[out] M
 179 *> \verbatim
 180 *>          M is INTEGER
 181 *>          The total number of eigenvalues (of all L_i D_i L_i^T)
 182 *>          found.
 183 *> \endverbatim
 184 *>
 185 *> \param[out] W
 186 *> \verbatim
 187 *>          W is REAL array, dimension (N)
 188 *>          The first M elements contain the eigenvalues. The
 189 *>          eigenvalues of each of the blocks, L_i D_i L_i^T, are
 190 *>          sorted in ascending order ( SLARRE may use the
 191 *>          remaining N-M elements as workspace).
 192 *> \endverbatim
 193 *>
 194 *> \param[out] WERR
 195 *> \verbatim
 196 *>          WERR is REAL array, dimension (N)
 197 *>          The error bound on the corresponding eigenvalue in W.
 198 *> \endverbatim
 199 *>
 200 *> \param[out] WGAP
 201 *> \verbatim
 202 *>          WGAP is REAL array, dimension (N)
 203 *>          The separation from the right neighbor eigenvalue in W.
 204 *>          The gap is only with respect to the eigenvalues of the same block
 205 *>          as each block has its own representation tree.
 206 *>          Exception: at the right end of a block we store the left gap
 207 *> \endverbatim
 208 *>
 209 *> \param[out] IBLOCK
 210 *> \verbatim
 211 *>          IBLOCK is INTEGER array, dimension (N)
 212 *>          The indices of the blocks (submatrices) associated with the
 213 *>          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
 214 *>          W(i) belongs to the first block from the top, =2 if W(i)
 215 *>          belongs to the second block, etc.
 216 *> \endverbatim
 217 *>
 218 *> \param[out] INDEXW
 219 *> \verbatim
 220 *>          INDEXW is INTEGER array, dimension (N)
 221 *>          The indices of the eigenvalues within each block (submatrix);
 222 *>          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
 223 *>          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
 224 *> \endverbatim
 225 *>
 226 *> \param[out] GERS
 227 *> \verbatim
 228 *>          GERS is REAL array, dimension (2*N)
 229 *>          The N Gerschgorin intervals (the i-th Gerschgorin interval
 230 *>          is (GERS(2*i-1), GERS(2*i)).
 231 *> \endverbatim
 232 *>
 233 *> \param[out] PIVMIN
 234 *> \verbatim
 235 *>          PIVMIN is REAL
 236 *>          The minimum pivot in the Sturm sequence for T.
 237 *> \endverbatim
 238 *>
 239 *> \param[out] WORK
 240 *> \verbatim
 241 *>          WORK is REAL array, dimension (6*N)
 242 *>          Workspace.
 243 *> \endverbatim
 244 *>
 245 *> \param[out] IWORK
 246 *> \verbatim
 247 *>          IWORK is INTEGER array, dimension (5*N)
 248 *>          Workspace.
 249 *> \endverbatim
 250 *>
 251 *> \param[out] INFO
 252 *> \verbatim
 253 *>          INFO is INTEGER
 254 *>          = 0:  successful exit
 255 *>          > 0:  A problem occurred in SLARRE.
 256 *>          < 0:  One of the called subroutines signaled an internal problem.
 257 *>                Needs inspection of the corresponding parameter IINFO
 258 *>                for further information.
 259 *>
 260 *>          =-1:  Problem in SLARRD.
 261 *>          = 2:  No base representation could be found in MAXTRY iterations.
 262 *>                Increasing MAXTRY and recompilation might be a remedy.
 263 *>          =-3:  Problem in SLARRB when computing the refined root
 264 *>                representation for SLASQ2.
 265 *>          =-4:  Problem in SLARRB when preforming bisection on the
 266 *>                desired part of the spectrum.
 267 *>          =-5:  Problem in SLASQ2.
 268 *>          =-6:  Problem in SLASQ2.
 269 *> \endverbatim
 270 *
 271 *  Authors:
 272 *  ========
 273 *
 274 *> \author Univ. of Tennessee
 275 *> \author Univ. of California Berkeley
 276 *> \author Univ. of Colorado Denver
 277 *> \author NAG Ltd.
 278 *
 279 *> \date June 2016
 280 *
 281 *> \ingroup OTHERauxiliary
 282 *
 283 *> \par Further Details:
 284 *  =====================
 285 *>
 286 *> \verbatim
 287 *>
 288 *>  The base representations are required to suffer very little
 289 *>  element growth and consequently define all their eigenvalues to
 290 *>  high relative accuracy.
 291 *> \endverbatim
 292 *
 293 *> \par Contributors:
 294 *  ==================
 295 *>
 296 *>     Beresford Parlett, University of California, Berkeley, USA \n
 297 *>     Jim Demmel, University of California, Berkeley, USA \n
 298 *>     Inderjit Dhillon, University of Texas, Austin, USA \n
 299 *>     Osni Marques, LBNL/NERSC, USA \n
 300 *>     Christof Voemel, University of California, Berkeley, USA \n
 301 *>
 302 *  =====================================================================
 303       SUBROUTINE SLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
 304      $                    RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
 305      $                    W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
 306      $                    WORK, IWORK, INFO )
 307 *
 308 *  -- LAPACK auxiliary routine (version 3.6.1) --
 309 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 310 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 311 *     June 2016
 312 *
 313 *     .. Scalar Arguments ..
 314       CHARACTER          RANGE
 315       INTEGER            IL, INFO, IU, M, N, NSPLIT
 316       REAL               PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
 317 *     ..
 318 *     .. Array Arguments ..
 319       INTEGER            IBLOCK( * ), ISPLIT( * ), IWORK( * ),
 320      $                   INDEXW( * )
 321       REAL               D( * ), E( * ), E2( * ), GERS( * ),
 322      $                   W( * ),WERR( * ), WGAP( * ), WORK( * )
 323 *     ..
 324 *
 325 *  =====================================================================
 326 *
 327 *     .. Parameters ..
 328       REAL               FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
 329      $                   MAXGROWTH, ONE, PERT, TWO, ZERO
 330       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0,
 331      $                     TWO = 2.0E0, FOUR=4.0E0,
 332      $                     HNDRD = 100.0E0,
 333      $                     PERT = 4.0E0,
 334      $                     HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF,
 335      $                     MAXGROWTH = 64.0E0, FUDGE = 2.0E0 )
 336       INTEGER            MAXTRY, ALLRNG, INDRNG, VALRNG
 337       PARAMETER          ( MAXTRY = 6, ALLRNG = 1, INDRNG = 2,
 338      $                     VALRNG = 3 )
 339 *     ..
 340 *     .. Local Scalars ..
 341       LOGICAL            FORCEB, NOREP, USEDQD
 342       INTEGER            CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
 343      $                   IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
 344      $                   WBEGIN, WEND
 345       REAL               AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
 346      $                   EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
 347      $                   RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
 348      $                   TAU, TMP, TMP1
 349
 350
 351 *     ..
 352 *     .. Local Arrays ..
 353       INTEGER            ISEED( 4 )
 354 *     ..
 355 *     .. External Functions ..
 356       LOGICAL            LSAME
 357       REAL                        SLAMCH
 358       EXTERNAL           SLAMCH, LSAME
 359
 360 *     ..
 361 *     .. External Subroutines ..
 362       EXTERNAL           SCOPY, SLARNV, SLARRA, SLARRB, SLARRC, SLARRD,
 363      $                   SLASQ2
 364 *     ..
 365 *     .. Intrinsic Functions ..
 366       INTRINSIC          ABS, MAX, MIN
 367
 368 *     ..
 369 *     .. Executable Statements ..
 370 *
 371
 372       INFO = 0
 373
 374 *
 375 *     Decode RANGE
 376 *
 377       IF( LSAME( RANGE, 'A' ) ) THEN
 378          IRANGE = ALLRNG
 379       ELSE IF( LSAME( RANGE, 'V' ) ) THEN
 380          IRANGE = VALRNG
 381       ELSE IF( LSAME( RANGE, 'I' ) ) THEN
 382          IRANGE = INDRNG
 383       END IF
 384
 385       M = 0
 386
 387 *     Get machine constants
 388       SAFMIN = SLAMCH( 'S' )
 389       EPS = SLAMCH( 'P' )
 390
 391 *     Set parameters
 392       RTL = HNDRD*EPS
 393 *     If one were ever to ask for less initial precision in BSRTOL,
 394 *     one should keep in mind that for the subset case, the extremal
 395 *     eigenvalues must be at least as accurate as the current setting
 396 *     (eigenvalues in the middle need not as much accuracy)
 397       BSRTOL = SQRT(EPS)*(0.5E-3)
 398
 399 *     Treat case of 1x1 matrix for quick return
 400       IF( N.EQ.1 ) THEN
 401          IF( (IRANGE.EQ.ALLRNG).OR.
 402      $       ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR.
 403      $       ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN
 404             M = 1
 405             W(1) = D(1)
 406 *           The computation error of the eigenvalue is zero
 407             WERR(1) = ZERO
 408             WGAP(1) = ZERO
 409             IBLOCK( 1 ) = 1
 410             INDEXW( 1 ) = 1
 411             GERS(1) = D( 1 )
 412             GERS(2) = D( 1 )
 413          ENDIF
 414 *        store the shift for the initial RRR, which is zero in this case
 415          E(1) = ZERO
 416          RETURN
 417       END IF
 418
 419 *     General case: tridiagonal matrix of order > 1
 420 *
 421 *     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
 422 *     Compute maximum off-diagonal entry and pivmin.
 423       GL = D(1)
 424       GU = D(1)
 425       EOLD = ZERO
 426       EMAX = ZERO
 427       E(N) = ZERO
 428       DO 5 I = 1,N
 429          WERR(I) = ZERO
 430          WGAP(I) = ZERO
 431          EABS = ABS( E(I) )
 432          IF( EABS .GE. EMAX ) THEN
 433             EMAX = EABS
 434          END IF
 435          TMP1 = EABS + EOLD
 436          GERS( 2*I-1) = D(I) - TMP1
 437          GL =  MIN( GL, GERS( 2*I - 1))
 438          GERS( 2*I ) = D(I) + TMP1
 439          GU = MAX( GU, GERS(2*I) )
 440          EOLD  = EABS
 441  5    CONTINUE
 442 *     The minimum pivot allowed in the Sturm sequence for T
 443       PIVMIN = SAFMIN * MAX( ONE, EMAX**2 )
 444 *     Compute spectral diameter. The Gerschgorin bounds give an
 445 *     estimate that is wrong by at most a factor of SQRT(2)
 446       SPDIAM = GU - GL
 447
 448 *     Compute splitting points
 449       CALL SLARRA( N, D, E, E2, SPLTOL, SPDIAM,
 450      $                    NSPLIT, ISPLIT, IINFO )
 451
 452 *     Can force use of bisection instead of faster DQDS.
 453 *     Option left in the code for future multisection work.
 454       FORCEB = .FALSE.
 455
 456 *     Initialize USEDQD, DQDS should be used for ALLRNG unless someone
 457 *     explicitly wants bisection.
 458       USEDQD = (( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB))
 459
 460       IF( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ) THEN
 461 *        Set interval [VL,VU] that contains all eigenvalues
 462          VL = GL
 463          VU = GU
 464       ELSE
 465 *        We call SLARRD to find crude approximations to the eigenvalues
 466 *        in the desired range. In case IRANGE = INDRNG, we also obtain the
 467 *        interval (VL,VU] that contains all the wanted eigenvalues.
 468 *        An interval [LEFT,RIGHT] has converged if
 469 *        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
 470 *        SLARRD needs a WORK of size 4*N, IWORK of size 3*N
 471          CALL SLARRD( RANGE, 'B', N, VL, VU, IL, IU, GERS,
 472      $                    BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
 473      $                    MM, W, WERR, VL, VU, IBLOCK, INDEXW,
 474      $                    WORK, IWORK, IINFO )
 475          IF( IINFO.NE.0 ) THEN
 476             INFO = -1
 477             RETURN
 478          ENDIF
 479 *        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
 480          DO 14 I = MM+1,N
 481             W( I ) = ZERO
 482             WERR( I ) = ZERO
 483             IBLOCK( I ) = 0
 484             INDEXW( I ) = 0
 485  14      CONTINUE
 486       END IF
 487
 488
 489 ***
 490 *     Loop over unreduced blocks
 491       IBEGIN = 1
 492       WBEGIN = 1
 493       DO 170 JBLK = 1, NSPLIT
 494          IEND = ISPLIT( JBLK )
 495          IN = IEND - IBEGIN + 1
 496
 497 *        1 X 1 block
 498          IF( IN.EQ.1 ) THEN
 499             IF( (IRANGE.EQ.ALLRNG).OR.( (IRANGE.EQ.VALRNG).AND.
 500      $         ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) )
 501      $        .OR. ( (IRANGE.EQ.INDRNG).AND.(IBLOCK(WBEGIN).EQ.JBLK))
 502      $        ) THEN
 503                M = M + 1
 504                W( M ) = D( IBEGIN )
 505                WERR(M) = ZERO
 506 *              The gap for a single block doesn't matter for the later
 507 *              algorithm and is assigned an arbitrary large value
 508                WGAP(M) = ZERO
 509                IBLOCK( M ) = JBLK
 510                INDEXW( M ) = 1
 511                WBEGIN = WBEGIN + 1
 512             ENDIF
 513 *           E( IEND ) holds the shift for the initial RRR
 514             E( IEND ) = ZERO
 515             IBEGIN = IEND + 1
 516             GO TO 170
 517          END IF
 518 *
 519 *        Blocks of size larger than 1x1
 520 *
 521 *        E( IEND ) will hold the shift for the initial RRR, for now set it =0
 522          E( IEND ) = ZERO
 523 *
 524 *        Find local outer bounds GL,GU for the block
 525          GL = D(IBEGIN)
 526          GU = D(IBEGIN)
 527          DO 15 I = IBEGIN , IEND
 528             GL = MIN( GERS( 2*I-1 ), GL )
 529             GU = MAX( GERS( 2*I ), GU )
 530  15      CONTINUE
 531          SPDIAM = GU - GL
 532
 533          IF(.NOT. ((IRANGE.EQ.ALLRNG).AND.(.NOT.FORCEB)) ) THEN
 534 *           Count the number of eigenvalues in the current block.
 535             MB = 0
 536             DO 20 I = WBEGIN,MM
 537                IF( IBLOCK(I).EQ.JBLK ) THEN
 538                   MB = MB+1
 539                ELSE
 540                   GOTO 21
 541                ENDIF
 542  20         CONTINUE
 543  21         CONTINUE
 544
 545             IF( MB.EQ.0) THEN
 546 *              No eigenvalue in the current block lies in the desired range
 547 *              E( IEND ) holds the shift for the initial RRR
 548                E( IEND ) = ZERO
 549                IBEGIN = IEND + 1
 550                GO TO 170
 551             ELSE
 552
 553 *              Decide whether dqds or bisection is more efficient
 554                USEDQD = ( (MB .GT. FAC*IN) .AND. (.NOT.FORCEB) )
 555                WEND = WBEGIN + MB - 1
 556 *              Calculate gaps for the current block
 557 *              In later stages, when representations for individual
 558 *              eigenvalues are different, we use SIGMA = E( IEND ).
 559                SIGMA = ZERO
 560                DO 30 I = WBEGIN, WEND - 1
 561                   WGAP( I ) = MAX( ZERO,
 562      $                        W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
 563  30            CONTINUE
 564                WGAP( WEND ) = MAX( ZERO,
 565      $                     VU - SIGMA - (W( WEND )+WERR( WEND )))
 566 *              Find local index of the first and last desired evalue.
 567                INDL = INDEXW(WBEGIN)
 568                INDU = INDEXW( WEND )
 569             ENDIF
 570          ENDIF
 571          IF(( (IRANGE.EQ.ALLRNG) .AND. (.NOT. FORCEB) ).OR.USEDQD) THEN
 572 *           Case of DQDS
 573 *           Find approximations to the extremal eigenvalues of the block
 574             CALL SLARRK( IN, 1, GL, GU, D(IBEGIN),
 575      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
 576             IF( IINFO.NE.0 ) THEN
 577                INFO = -1
 578                RETURN
 579             ENDIF
 580             ISLEFT = MAX(GL, TMP - TMP1
 581      $               - HNDRD * EPS* ABS(TMP - TMP1))
 582
 583             CALL SLARRK( IN, IN, GL, GU, D(IBEGIN),
 584      $               E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO )
 585             IF( IINFO.NE.0 ) THEN
 586                INFO = -1
 587                RETURN
 588             ENDIF
 589             ISRGHT = MIN(GU, TMP + TMP1
 590      $                 + HNDRD * EPS * ABS(TMP + TMP1))
 591 *           Improve the estimate of the spectral diameter
 592             SPDIAM = ISRGHT - ISLEFT
 593          ELSE
 594 *           Case of bisection
 595 *           Find approximations to the wanted extremal eigenvalues
 596             ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN)
 597      $                  - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) ))
 598             ISRGHT = MIN(GU,W(WEND) + WERR(WEND)
 599      $                  + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND)))
 600          ENDIF
 601
 602
 603 *        Decide whether the base representation for the current block
 604 *        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
 605 *        should be on the left or the right end of the current block.
 606 *        The strategy is to shift to the end which is "more populated"
 607 *        Furthermore, decide whether to use DQDS for the computation of
 608 *        the eigenvalue approximations at the end of SLARRE or bisection.
 609 *        dqds is chosen if all eigenvalues are desired or the number of
 610 *        eigenvalues to be computed is large compared to the blocksize.
 611          IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
 612 *           If all the eigenvalues have to be computed, we use dqd
 613             USEDQD = .TRUE.
 614 *           INDL is the local index of the first eigenvalue to compute
 615             INDL = 1
 616             INDU = IN
 617 *           MB =  number of eigenvalues to compute
 618             MB = IN
 619             WEND = WBEGIN + MB - 1
 620 *           Define 1/4 and 3/4 points of the spectrum
 621             S1 = ISLEFT + FOURTH * SPDIAM
 622             S2 = ISRGHT - FOURTH * SPDIAM
 623          ELSE
 624 *           SLARRD has computed IBLOCK and INDEXW for each eigenvalue
 625 *           approximation.
 626 *           choose sigma
 627             IF( USEDQD ) THEN
 628                S1 = ISLEFT + FOURTH * SPDIAM
 629                S2 = ISRGHT - FOURTH * SPDIAM
 630             ELSE
 631                TMP = MIN(ISRGHT,VU) -  MAX(ISLEFT,VL)
 632                S1 =  MAX(ISLEFT,VL) + FOURTH * TMP
 633                S2 =  MIN(ISRGHT,VU) - FOURTH * TMP
 634             ENDIF
 635          ENDIF
 636
 637 *        Compute the negcount at the 1/4 and 3/4 points
 638          IF(MB.GT.1) THEN
 639             CALL SLARRC( 'T', IN, S1, S2, D(IBEGIN),
 640      $                    E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO)
 641          ENDIF
 642
 643          IF(MB.EQ.1) THEN
 644             SIGMA = GL
 645             SGNDEF = ONE
 646          ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN
 647             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
 648                SIGMA = MAX(ISLEFT,GL)
 649             ELSEIF( USEDQD ) THEN
 650 *              use Gerschgorin bound as shift to get pos def matrix
 651 *              for dqds
 652                SIGMA = ISLEFT
 653             ELSE
 654 *              use approximation of the first desired eigenvalue of the
 655 *              block as shift
 656                SIGMA = MAX(ISLEFT,VL)
 657             ENDIF
 658             SGNDEF = ONE
 659          ELSE
 660             IF( ( IRANGE.EQ.ALLRNG ) .AND. (.NOT.FORCEB) ) THEN
 661                SIGMA = MIN(ISRGHT,GU)
 662             ELSEIF( USEDQD ) THEN
 663 *              use Gerschgorin bound as shift to get neg def matrix
 664 *              for dqds
 665                SIGMA = ISRGHT
 666             ELSE
 667 *              use approximation of the first desired eigenvalue of the
 668 *              block as shift
 669                SIGMA = MIN(ISRGHT,VU)
 670             ENDIF
 671             SGNDEF = -ONE
 672          ENDIF
 673
 674
 675 *        An initial SIGMA has been chosen that will be used for computing
 676 *        T - SIGMA I = L D L^T
 677 *        Define the increment TAU of the shift in case the initial shift
 678 *        needs to be refined to obtain a factorization with not too much
 679 *        element growth.
 680          IF( USEDQD ) THEN
 681 *           The initial SIGMA was to the outer end of the spectrum
 682 *           the matrix is definite and we need not retreat.
 683             TAU = SPDIAM*EPS*N + TWO*PIVMIN
 684             TAU = MAX( TAU,TWO*EPS*ABS(SIGMA) )
 685          ELSE
 686             IF(MB.GT.1) THEN
 687                CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN)
 688                AVGAP = ABS(CLWDTH / REAL(WEND-WBEGIN))
 689                IF( SGNDEF.EQ.ONE ) THEN
 690                   TAU = HALF*MAX(WGAP(WBEGIN),AVGAP)
 691                   TAU = MAX(TAU,WERR(WBEGIN))
 692                ELSE
 693                   TAU = HALF*MAX(WGAP(WEND-1),AVGAP)
 694                   TAU = MAX(TAU,WERR(WEND))
 695                ENDIF
 696             ELSE
 697                TAU = WERR(WBEGIN)
 698             ENDIF
 699          ENDIF
 700 *
 701          DO 80 IDUM = 1, MAXTRY
 702 *           Compute L D L^T factorization of tridiagonal matrix T - sigma I.
 703 *           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
 704 *           pivots in WORK(2*IN+1:3*IN)
 705             DPIVOT = D( IBEGIN ) - SIGMA
 706             WORK( 1 ) = DPIVOT
 707             DMAX = ABS( WORK(1) )
 708             J = IBEGIN
 709             DO 70 I = 1, IN - 1
 710                WORK( 2*IN+I ) = ONE / WORK( I )
 711                TMP = E( J )*WORK( 2*IN+I )
 712                WORK( IN+I ) = TMP
 713                DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J )
 714                WORK( I+1 ) = DPIVOT
 715                DMAX = MAX( DMAX, ABS(DPIVOT) )
 716                J = J + 1
 717  70         CONTINUE
 718 *           check for element growth
 719             IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN
 720                NOREP = .TRUE.
 721             ELSE
 722                NOREP = .FALSE.
 723             ENDIF
 724             IF( USEDQD .AND. .NOT.NOREP ) THEN
 725 *              Ensure the definiteness of the representation
 726 *              All entries of D (of L D L^T) must have the same sign
 727                DO 71 I = 1, IN
 728                   TMP = SGNDEF*WORK( I )
 729                   IF( TMP.LT.ZERO ) NOREP = .TRUE.
 730  71            CONTINUE
 731             ENDIF
 732             IF(NOREP) THEN
 733 *              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
 734 *              shift which makes the matrix definite. So we should end up
 735 *              here really only in the case of IRANGE = VALRNG or INDRNG.
 736                IF( IDUM.EQ.MAXTRY-1 ) THEN
 737                   IF( SGNDEF.EQ.ONE ) THEN
 738 *                    The fudged Gerschgorin shift should succeed
 739                      SIGMA =
 740      $                    GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN
 741                   ELSE
 742                      SIGMA =
 743      $                    GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN
 744                   END IF
 745                ELSE
 746                   SIGMA = SIGMA - SGNDEF * TAU
 747                   TAU = TWO * TAU
 748                END IF
 749             ELSE
 750 *              an initial RRR is found
 751                GO TO 83
 752             END IF
 753  80      CONTINUE
 754 *        if the program reaches this point, no base representation could be
 755 *        found in MAXTRY iterations.
 756          INFO = 2
 757          RETURN
 758
 759  83      CONTINUE
 760 *        At this point, we have found an initial base representation
 761 *        T - SIGMA I = L D L^T with not too much element growth.
 762 *        Store the shift.
 763          E( IEND ) = SIGMA
 764 *        Store D and L.
 765          CALL SCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
 766          CALL SCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
 767
 768
 769          IF(MB.GT.1 ) THEN
 770 *
 771 *           Perturb each entry of the base representation by a small
 772 *           (but random) relative amount to overcome difficulties with
 773 *           glued matrices.
 774 *
 775             DO 122 I = 1, 4
 776                ISEED( I ) = 1
 777  122        CONTINUE
 778
 779             CALL SLARNV(2, ISEED, 2*IN-1, WORK(1))
 780             DO 125 I = 1,IN-1
 781                D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(I))
 782                E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(IN+I))
 783  125        CONTINUE
 784             D(IEND) = D(IEND)*(ONE+EPS*FOUR*WORK(IN))
 785 *
 786          ENDIF
 787 *
 788 *        Don't update the Gerschgorin intervals because keeping track
 789 *        of the updates would be too much work in SLARRV.
 790 *        We update W instead and use it to locate the proper Gerschgorin
 791 *        intervals.
 792
 793 *        Compute the required eigenvalues of L D L' by bisection or dqds
 794          IF ( .NOT.USEDQD ) THEN
 795 *           If SLARRD has been used, shift the eigenvalue approximations
 796 *           according to their representation. This is necessary for
 797 *           a uniform SLARRV since dqds computes eigenvalues of the
 798 *           shifted representation. In SLARRV, W will always hold the
 799 *           UNshifted eigenvalue approximation.
 800             DO 134 J=WBEGIN,WEND
 801                W(J) = W(J) - SIGMA
 802                WERR(J) = WERR(J) + ABS(W(J)) * EPS
 803  134        CONTINUE
 804 *           call SLARRB to reduce eigenvalue error of the approximations
 805 *           from SLARRD
 806             DO 135 I = IBEGIN, IEND-1
 807                WORK( I ) = D( I ) * E( I )**2
 808  135        CONTINUE
 809 *           use bisection to find EV from INDL to INDU
 810             CALL SLARRB(IN, D(IBEGIN), WORK(IBEGIN),
 811      $                  INDL, INDU, RTOL1, RTOL2, INDL-1,
 812      $                  W(WBEGIN), WGAP(WBEGIN), WERR(WBEGIN),
 813      $                  WORK( 2*N+1 ), IWORK, PIVMIN, SPDIAM,
 814      $                  IN, IINFO )
 815             IF( IINFO .NE. 0 ) THEN
 816                INFO = -4
 817                RETURN
 818             END IF
 819 *           SLARRB computes all gaps correctly except for the last one
 820 *           Record distance to VU/GU
 821             WGAP( WEND ) = MAX( ZERO,
 822      $           ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) )
 823             DO 138 I = INDL, INDU
 824                M = M + 1
 825                IBLOCK(M) = JBLK
 826                INDEXW(M) = I
 827  138        CONTINUE
 828          ELSE
 829 *           Call dqds to get all eigs (and then possibly delete unwanted
 830 *           eigenvalues).
 831 *           Note that dqds finds the eigenvalues of the L D L^T representation
 832 *           of T to high relative accuracy. High relative accuracy
 833 *           might be lost when the shift of the RRR is subtracted to obtain
 834 *           the eigenvalues of T. However, T is not guaranteed to define its
 835 *           eigenvalues to high relative accuracy anyway.
 836 *           Set RTOL to the order of the tolerance used in SLASQ2
 837 *           This is an ESTIMATED error, the worst case bound is 4*N*EPS
 838 *           which is usually too large and requires unnecessary work to be
 839 *           done by bisection when computing the eigenvectors
 840             RTOL = LOG(REAL(IN)) * FOUR * EPS
 841             J = IBEGIN
 842             DO 140 I = 1, IN - 1
 843                WORK( 2*I-1 ) = ABS( D( J ) )
 844                WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
 845                J = J + 1
 846   140       CONTINUE
 847             WORK( 2*IN-1 ) = ABS( D( IEND ) )
 848             WORK( 2*IN ) = ZERO
 849             CALL SLASQ2( IN, WORK, IINFO )
 850             IF( IINFO .NE. 0 ) THEN
 851 *              If IINFO = -5 then an index is part of a tight cluster
 852 *              and should be changed. The index is in IWORK(1) and the
 853 *              gap is in WORK(N+1)
 854                INFO = -5
 855                RETURN
 856             ELSE
 857 *              Test that all eigenvalues are positive as expected
 858                DO 149 I = 1, IN
 859                   IF( WORK( I ).LT.ZERO ) THEN
 860                      INFO = -6
 861                      RETURN
 862                   ENDIF
 863  149           CONTINUE
 864             END IF
 865             IF( SGNDEF.GT.ZERO ) THEN
 866                DO 150 I = INDL, INDU
 867                   M = M + 1
 868                   W( M ) = WORK( IN-I+1 )
 869                   IBLOCK( M ) = JBLK
 870                   INDEXW( M ) = I
 871  150           CONTINUE
 872             ELSE
 873                DO 160 I = INDL, INDU
 874                   M = M + 1
 875                   W( M ) = -WORK( I )
 876                   IBLOCK( M ) = JBLK
 877                   INDEXW( M ) = I
 878  160           CONTINUE
 879             END IF
 880
 881             DO 165 I = M - MB + 1, M
 882 *              the value of RTOL below should be the tolerance in SLASQ2
 883                WERR( I ) = RTOL * ABS( W(I) )
 884  165        CONTINUE
 885             DO 166 I = M - MB + 1, M - 1
 886 *              compute the right gap between the intervals
 887                WGAP( I ) = MAX( ZERO,
 888      $                          W(I+1)-WERR(I+1) - (W(I)+WERR(I)) )
 889  166        CONTINUE
 890             WGAP( M ) = MAX( ZERO,
 891      $           ( VU-SIGMA ) - ( W( M ) + WERR( M ) ) )
 892          END IF
 893 *        proceed with next block
 894          IBEGIN = IEND + 1
 895          WBEGIN = WEND + 1
 896  170  CONTINUE
 897 *
 898
 899       RETURN
 900 *
 901 *     end of SLARRE
 902 *
 903       END