SRC/dlarrd.f

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