SRC/zstemr.f

   1 *> \brief \b ZSTEMR
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZSTEMR + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
  22 *                          M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
  23 *                          IWORK, LIWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBZ, RANGE
  27 *       LOGICAL            TRYRAC
  28 *       INTEGER            IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
  29 *       DOUBLE PRECISION VL, VU
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            ISUPPZ( * ), IWORK( * )
  33 *       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
  34 *       COMPLEX*16         Z( LDZ, * )
  35 *       ..
  36 *
  37 *
  38 *> \par Purpose:
  39 *  =============
  40 *>
  41 *> \verbatim
  42 *>
  43 *> ZSTEMR computes selected eigenvalues and, optionally, eigenvectors
  44 *> of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
  45 *> a well defined set of pairwise different real eigenvalues, the corresponding
  46 *> real eigenvectors are pairwise orthogonal.
  47 *>
  48 *> The spectrum may be computed either completely or partially by specifying
  49 *> either an interval (VL,VU] or a range of indices IL:IU for the desired
  50 *> eigenvalues.
  51 *>
  52 *> Depending on the number of desired eigenvalues, these are computed either
  53 *> by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are
  54 *> computed by the use of various suitable L D L^T factorizations near clusters
  55 *> of close eigenvalues (referred to as RRRs, Relatively Robust
  56 *> Representations). An informal sketch of the algorithm follows.
  57 *>
  58 *> For each unreduced block (submatrix) of T,
  59 *>    (a) Compute T - sigma I  = L D L^T, so that L and D
  60 *>        define all the wanted eigenvalues to high relative accuracy.
  61 *>        This means that small relative changes in the entries of D and L
  62 *>        cause only small relative changes in the eigenvalues and
  63 *>        eigenvectors. The standard (unfactored) representation of the
  64 *>        tridiagonal matrix T does not have this property in general.
  65 *>    (b) Compute the eigenvalues to suitable accuracy.
  66 *>        If the eigenvectors are desired, the algorithm attains full
  67 *>        accuracy of the computed eigenvalues only right before
  68 *>        the corresponding vectors have to be computed, see steps c) and d).
  69 *>    (c) For each cluster of close eigenvalues, select a new
  70 *>        shift close to the cluster, find a new factorization, and refine
  71 *>        the shifted eigenvalues to suitable accuracy.
  72 *>    (d) For each eigenvalue with a large enough relative separation compute
  73 *>        the corresponding eigenvector by forming a rank revealing twisted
  74 *>        factorization. Go back to (c) for any clusters that remain.
  75 *>
  76 *> For more details, see:
  77 *> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
  78 *>   to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
  79 *>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
  80 *> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
  81 *>   Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
  82 *>   2004.  Also LAPACK Working Note 154.
  83 *> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
  84 *>   tridiagonal eigenvalue/eigenvector problem",
  85 *>   Computer Science Division Technical Report No. UCB/CSD-97-971,
  86 *>   UC Berkeley, May 1997.
  87 *>
  88 *> Further Details
  89 *> 1.ZSTEMR works only on machines which follow IEEE-754
  90 *> floating-point standard in their handling of infinities and NaNs.
  91 *> This permits the use of efficient inner loops avoiding a check for
  92 *> zero divisors.
  93 *>
  94 *> 2. LAPACK routines can be used to reduce a complex Hermitean matrix to
  95 *> real symmetric tridiagonal form.
  96 *>
  97 *> (Any complex Hermitean tridiagonal matrix has real values on its diagonal
  98 *> and potentially complex numbers on its off-diagonals. By applying a
  99 *> similarity transform with an appropriate diagonal matrix
 100 *> diag(1,e^{i \phy_1}, ... , e^{i \phy_{n-1}}), the complex Hermitean
 101 *> matrix can be transformed into a real symmetric matrix and complex
 102 *> arithmetic can be entirely avoided.)
 103 *>
 104 *> While the eigenvectors of the real symmetric tridiagonal matrix are real,
 105 *> the eigenvectors of original complex Hermitean matrix have complex entries
 106 *> in general.
 107 *> Since LAPACK drivers overwrite the matrix data with the eigenvectors,
 108 *> ZSTEMR accepts complex workspace to facilitate interoperability
 109 *> with ZUNMTR or ZUPMTR.
 110 *> \endverbatim
 111 *
 112 *  Arguments:
 113 *  ==========
 114 *
 115 *> \param[in] JOBZ
 116 *> \verbatim
 117 *>          JOBZ is CHARACTER*1
 118 *>          = 'N':  Compute eigenvalues only;
 119 *>          = 'V':  Compute eigenvalues and eigenvectors.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] RANGE
 123 *> \verbatim
 124 *>          RANGE is CHARACTER*1
 125 *>          = 'A': all eigenvalues will be found.
 126 *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
 127 *>                 will be found.
 128 *>          = 'I': the IL-th through IU-th eigenvalues will be found.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] N
 132 *> \verbatim
 133 *>          N is INTEGER
 134 *>          The order of the matrix.  N >= 0.
 135 *> \endverbatim
 136 *>
 137 *> \param[in,out] D
 138 *> \verbatim
 139 *>          D is DOUBLE PRECISION array, dimension (N)
 140 *>          On entry, the N diagonal elements of the tridiagonal matrix
 141 *>          T. On exit, D is overwritten.
 142 *> \endverbatim
 143 *>
 144 *> \param[in,out] E
 145 *> \verbatim
 146 *>          E is DOUBLE PRECISION array, dimension (N)
 147 *>          On entry, the (N-1) subdiagonal elements of the tridiagonal
 148 *>          matrix T in elements 1 to N-1 of E. E(N) need not be set on
 149 *>          input, but is used internally as workspace.
 150 *>          On exit, E is overwritten.
 151 *> \endverbatim
 152 *>
 153 *> \param[in] VL
 154 *> \verbatim
 155 *>          VL is DOUBLE PRECISION
 156 *>
 157 *>          If RANGE='V', the lower bound of the interval to
 158 *>          be searched for eigenvalues. VL < VU.
 159 *>          Not referenced if RANGE = 'A' or 'I'.
 160 *> \endverbatim
 161 *>
 162 *> \param[in] VU
 163 *> \verbatim
 164 *>          VU is DOUBLE PRECISION
 165 *>
 166 *>          If RANGE='V', the upper bound of the interval to
 167 *>          be searched for eigenvalues. VL < VU.
 168 *>          Not referenced if RANGE = 'A' or 'I'.
 169 *> \endverbatim
 170 *>
 171 *> \param[in] IL
 172 *> \verbatim
 173 *>          IL is INTEGER
 174 *>
 175 *>          If RANGE='I', the index of the
 176 *>          smallest eigenvalue to be returned.
 177 *>          1 <= IL <= IU <= N, if N > 0.
 178 *>          Not referenced if RANGE = 'A' or 'V'.
 179 *> \endverbatim
 180 *>
 181 *> \param[in] IU
 182 *> \verbatim
 183 *>          IU is INTEGER
 184 *>
 185 *>          If RANGE='I', the index of the
 186 *>          largest eigenvalue to be returned.
 187 *>          1 <= IL <= IU <= N, if N > 0.
 188 *>          Not referenced if RANGE = 'A' or 'V'.
 189 *> \endverbatim
 190 *>
 191 *> \param[out] M
 192 *> \verbatim
 193 *>          M is INTEGER
 194 *>          The total number of eigenvalues found.  0 <= M <= N.
 195 *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
 196 *> \endverbatim
 197 *>
 198 *> \param[out] W
 199 *> \verbatim
 200 *>          W is DOUBLE PRECISION array, dimension (N)
 201 *>          The first M elements contain the selected eigenvalues in
 202 *>          ascending order.
 203 *> \endverbatim
 204 *>
 205 *> \param[out] Z
 206 *> \verbatim
 207 *>          Z is COMPLEX*16 array, dimension (LDZ, max(1,M) )
 208 *>          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
 209 *>          contain the orthonormal eigenvectors of the matrix T
 210 *>          corresponding to the selected eigenvalues, with the i-th
 211 *>          column of Z holding the eigenvector associated with W(i).
 212 *>          If JOBZ = 'N', then Z is not referenced.
 213 *>          Note: the user must ensure that at least max(1,M) columns are
 214 *>          supplied in the array Z; if RANGE = 'V', the exact value of M
 215 *>          is not known in advance and can be computed with a workspace
 216 *>          query by setting NZC = -1, see below.
 217 *> \endverbatim
 218 *>
 219 *> \param[in] LDZ
 220 *> \verbatim
 221 *>          LDZ is INTEGER
 222 *>          The leading dimension of the array Z.  LDZ >= 1, and if
 223 *>          JOBZ = 'V', then LDZ >= max(1,N).
 224 *> \endverbatim
 225 *>
 226 *> \param[in] NZC
 227 *> \verbatim
 228 *>          NZC is INTEGER
 229 *>          The number of eigenvectors to be held in the array Z.
 230 *>          If RANGE = 'A', then NZC >= max(1,N).
 231 *>          If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
 232 *>          If RANGE = 'I', then NZC >= IU-IL+1.
 233 *>          If NZC = -1, then a workspace query is assumed; the
 234 *>          routine calculates the number of columns of the array Z that
 235 *>          are needed to hold the eigenvectors.
 236 *>          This value is returned as the first entry of the Z array, and
 237 *>          no error message related to NZC is issued by XERBLA.
 238 *> \endverbatim
 239 *>
 240 *> \param[out] ISUPPZ
 241 *> \verbatim
 242 *>          ISUPPZ is INTEGER ARRAY, dimension ( 2*max(1,M) )
 243 *>          The support of the eigenvectors in Z, i.e., the indices
 244 *>          indicating the nonzero elements in Z. The i-th computed eigenvector
 245 *>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
 246 *>          ISUPPZ( 2*i ). This is relevant in the case when the matrix
 247 *>          is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
 248 *> \endverbatim
 249 *>
 250 *> \param[in,out] TRYRAC
 251 *> \verbatim
 252 *>          TRYRAC is LOGICAL
 253 *>          If TRYRAC.EQ..TRUE., indicates that the code should check whether
 254 *>          the tridiagonal matrix defines its eigenvalues to high relative
 255 *>          accuracy.  If so, the code uses relative-accuracy preserving
 256 *>          algorithms that might be (a bit) slower depending on the matrix.
 257 *>          If the matrix does not define its eigenvalues to high relative
 258 *>          accuracy, the code can uses possibly faster algorithms.
 259 *>          If TRYRAC.EQ..FALSE., the code is not required to guarantee
 260 *>          relatively accurate eigenvalues and can use the fastest possible
 261 *>          techniques.
 262 *>          On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
 263 *>          does not define its eigenvalues to high relative accuracy.
 264 *> \endverbatim
 265 *>
 266 *> \param[out] WORK
 267 *> \verbatim
 268 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 269 *>          On exit, if INFO = 0, WORK(1) returns the optimal
 270 *>          (and minimal) LWORK.
 271 *> \endverbatim
 272 *>
 273 *> \param[in] LWORK
 274 *> \verbatim
 275 *>          LWORK is INTEGER
 276 *>          The dimension of the array WORK. LWORK >= max(1,18*N)
 277 *>          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
 278 *>          If LWORK = -1, then a workspace query is assumed; the routine
 279 *>          only calculates the optimal size of the WORK array, returns
 280 *>          this value as the first entry of the WORK array, and no error
 281 *>          message related to LWORK is issued by XERBLA.
 282 *> \endverbatim
 283 *>
 284 *> \param[out] IWORK
 285 *> \verbatim
 286 *>          IWORK is INTEGER array, dimension (LIWORK)
 287 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 288 *> \endverbatim
 289 *>
 290 *> \param[in] LIWORK
 291 *> \verbatim
 292 *>          LIWORK is INTEGER
 293 *>          The dimension of the array IWORK.  LIWORK >= max(1,10*N)
 294 *>          if the eigenvectors are desired, and LIWORK >= max(1,8*N)
 295 *>          if only the eigenvalues are to be computed.
 296 *>          If LIWORK = -1, then a workspace query is assumed; the
 297 *>          routine only calculates the optimal size of the IWORK array,
 298 *>          returns this value as the first entry of the IWORK array, and
 299 *>          no error message related to LIWORK is issued by XERBLA.
 300 *> \endverbatim
 301 *>
 302 *> \param[out] INFO
 303 *> \verbatim
 304 *>          INFO is INTEGER
 305 *>          On exit, INFO
 306 *>          = 0:  successful exit
 307 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 308 *>          > 0:  if INFO = 1X, internal error in DLARRE,
 309 *>                if INFO = 2X, internal error in ZLARRV.
 310 *>                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
 311 *>                the nonzero error code returned by DLARRE or
 312 *>                ZLARRV, respectively.
 313 *> \endverbatim
 314 *
 315 *  Authors:
 316 *  ========
 317 *
 318 *> \author Univ. of Tennessee
 319 *> \author Univ. of California Berkeley
 320 *> \author Univ. of Colorado Denver
 321 *> \author NAG Ltd.
 322 *
 323 *> \date June 2016
 324 *
 325 *> \ingroup complex16OTHERcomputational
 326 *
 327 *> \par Contributors:
 328 *  ==================
 329 *>
 330 *> Beresford Parlett, University of California, Berkeley, USA \n
 331 *> Jim Demmel, University of California, Berkeley, USA \n
 332 *> Inderjit Dhillon, University of Texas, Austin, USA \n
 333 *> Osni Marques, LBNL/NERSC, USA \n
 334 *> Christof Voemel, University of California, Berkeley, USA \n
 335 *
 336 *  =====================================================================
 337       SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
 338      $                   M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
 339      $                   IWORK, LIWORK, INFO )
 340 *
 341 *  -- LAPACK computational routine (version 3.6.1) --
 342 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 343 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 344 *     June 2016
 345 *
 346 *     .. Scalar Arguments ..
 347       CHARACTER          JOBZ, RANGE
 348       LOGICAL            TRYRAC
 349       INTEGER            IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
 350       DOUBLE PRECISION VL, VU
 351 *     ..
 352 *     .. Array Arguments ..
 353       INTEGER            ISUPPZ( * ), IWORK( * )
 354       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
 355       COMPLEX*16         Z( LDZ, * )
 356 *     ..
 357 *
 358 *  =====================================================================
 359 *
 360 *     .. Parameters ..
 361       DOUBLE PRECISION   ZERO, ONE, FOUR, MINRGP
 362       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0,
 363      $                     FOUR = 4.0D0,
 364      $                     MINRGP = 1.0D-3 )
 365 *     ..
 366 *     .. Local Scalars ..
 367       LOGICAL            ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
 368       INTEGER            I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
 369      $                   IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
 370      $                   INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
 371      $                   ITMP2, J, JBLK, JJ, LIWMIN, LWMIN, NSPLIT,
 372      $                   NZCMIN, OFFSET, WBEGIN, WEND
 373       DOUBLE PRECISION   BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN,
 374      $                   RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN,
 375      $                   THRESH, TMP, TNRM, WL, WU
 376 *     ..
 377 *     ..
 378 *     .. External Functions ..
 379       LOGICAL            LSAME
 380       DOUBLE PRECISION   DLAMCH, DLANST
 381       EXTERNAL           LSAME, DLAMCH, DLANST
 382 *     ..
 383 *     .. External Subroutines ..
 384       EXTERNAL           DCOPY, DLAE2, DLAEV2, DLARRC, DLARRE, DLARRJ,
 385      $                   DLARRR, DLASRT, DSCAL, XERBLA, ZLARRV, ZSWAP
 386 *     ..
 387 *     .. Intrinsic Functions ..
 388       INTRINSIC          MAX, MIN, SQRT
 389
 390
 391 *     ..
 392 *     .. Executable Statements ..
 393 *
 394 *     Test the input parameters.
 395 *
 396       WANTZ = LSAME( JOBZ, 'V' )
 397       ALLEIG = LSAME( RANGE, 'A' )
 398       VALEIG = LSAME( RANGE, 'V' )
 399       INDEIG = LSAME( RANGE, 'I' )
 400 *
 401       LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) )
 402       ZQUERY = ( NZC.EQ.-1 )
 403
 404 *     DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
 405 *     In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
 406 *     Furthermore, ZLARRV needs WORK of size 12*N, IWORK of size 7*N.
 407       IF( WANTZ ) THEN
 408          LWMIN = 18*N
 409          LIWMIN = 10*N
 410       ELSE
 411 *        need less workspace if only the eigenvalues are wanted
 412          LWMIN = 12*N
 413          LIWMIN = 8*N
 414       ENDIF
 415
 416       WL = ZERO
 417       WU = ZERO
 418       IIL = 0
 419       IIU = 0
 420       NSPLIT = 0
 421
 422       IF( VALEIG ) THEN
 423 *        We do not reference VL, VU in the cases RANGE = 'I','A'
 424 *        The interval (WL, WU] contains all the wanted eigenvalues.
 425 *        It is either given by the user or computed in DLARRE.
 426          WL = VL
 427          WU = VU
 428       ELSEIF( INDEIG ) THEN
 429 *        We do not reference IL, IU in the cases RANGE = 'V','A'
 430          IIL = IL
 431          IIU = IU
 432       ENDIF
 433 *
 434       INFO = 0
 435       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 436          INFO = -1
 437       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
 438          INFO = -2
 439       ELSE IF( N.LT.0 ) THEN
 440          INFO = -3
 441       ELSE IF( VALEIG .AND. N.GT.0 .AND. WU.LE.WL ) THEN
 442          INFO = -7
 443       ELSE IF( INDEIG .AND. ( IIL.LT.1 .OR. IIL.GT.N ) ) THEN
 444          INFO = -8
 445       ELSE IF( INDEIG .AND. ( IIU.LT.IIL .OR. IIU.GT.N ) ) THEN
 446          INFO = -9
 447       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 448          INFO = -13
 449       ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 450          INFO = -17
 451       ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 452          INFO = -19
 453       END IF
 454 *
 455 *     Get machine constants.
 456 *
 457       SAFMIN = DLAMCH( 'Safe minimum' )
 458       EPS = DLAMCH( 'Precision' )
 459       SMLNUM = SAFMIN / EPS
 460       BIGNUM = ONE / SMLNUM
 461       RMIN = SQRT( SMLNUM )
 462       RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
 463 *
 464       IF( INFO.EQ.0 ) THEN
 465          WORK( 1 ) = LWMIN
 466          IWORK( 1 ) = LIWMIN
 467 *
 468          IF( WANTZ .AND. ALLEIG ) THEN
 469             NZCMIN = N
 470          ELSE IF( WANTZ .AND. VALEIG ) THEN
 471             CALL DLARRC( 'T', N, VL, VU, D, E, SAFMIN,
 472      $                            NZCMIN, ITMP, ITMP2, INFO )
 473          ELSE IF( WANTZ .AND. INDEIG ) THEN
 474             NZCMIN = IIU-IIL+1
 475          ELSE
 476 *           WANTZ .EQ. FALSE.
 477             NZCMIN = 0
 478          ENDIF
 479          IF( ZQUERY .AND. INFO.EQ.0 ) THEN
 480             Z( 1,1 ) = NZCMIN
 481          ELSE IF( NZC.LT.NZCMIN .AND. .NOT.ZQUERY ) THEN
 482             INFO = -14
 483          END IF
 484       END IF
 485
 486       IF( INFO.NE.0 ) THEN
 487 *
 488          CALL XERBLA( 'ZSTEMR', -INFO )
 489 *
 490          RETURN
 491       ELSE IF( LQUERY .OR. ZQUERY ) THEN
 492          RETURN
 493       END IF
 494 *
 495 *     Handle N = 0, 1, and 2 cases immediately
 496 *
 497       M = 0
 498       IF( N.EQ.0 )
 499      $   RETURN
 500 *
 501       IF( N.EQ.1 ) THEN
 502          IF( ALLEIG .OR. INDEIG ) THEN
 503             M = 1
 504             W( 1 ) = D( 1 )
 505          ELSE
 506             IF( WL.LT.D( 1 ) .AND. WU.GE.D( 1 ) ) THEN
 507                M = 1
 508                W( 1 ) = D( 1 )
 509             END IF
 510          END IF
 511          IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
 512             Z( 1, 1 ) = ONE
 513             ISUPPZ(1) = 1
 514             ISUPPZ(2) = 1
 515          END IF
 516          RETURN
 517       END IF
 518 *
 519       IF( N.EQ.2 ) THEN
 520          IF( .NOT.WANTZ ) THEN
 521             CALL DLAE2( D(1), E(1), D(2), R1, R2 )
 522          ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
 523             CALL DLAEV2( D(1), E(1), D(2), R1, R2, CS, SN )
 524          END IF
 525          IF( ALLEIG.OR.
 526      $      (VALEIG.AND.(R2.GT.WL).AND.
 527      $                  (R2.LE.WU)).OR.
 528      $      (INDEIG.AND.(IIL.EQ.1)) ) THEN
 529             M = M+1
 530             W( M ) = R2
 531             IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
 532                Z( 1, M ) = -SN
 533                Z( 2, M ) = CS
 534 *              Note: At most one of SN and CS can be zero.
 535                IF (SN.NE.ZERO) THEN
 536                   IF (CS.NE.ZERO) THEN
 537                      ISUPPZ(2*M-1) = 1
 538                      ISUPPZ(2*M) = 2
 539                   ELSE
 540                      ISUPPZ(2*M-1) = 1
 541                      ISUPPZ(2*M) = 1
 542                   END IF
 543                ELSE
 544                   ISUPPZ(2*M-1) = 2
 545                   ISUPPZ(2*M) = 2
 546                END IF
 547             ENDIF
 548          ENDIF
 549          IF( ALLEIG.OR.
 550      $      (VALEIG.AND.(R1.GT.WL).AND.
 551      $                  (R1.LE.WU)).OR.
 552      $      (INDEIG.AND.(IIU.EQ.2)) ) THEN
 553             M = M+1
 554             W( M ) = R1
 555             IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN
 556                Z( 1, M ) = CS
 557                Z( 2, M ) = SN
 558 *              Note: At most one of SN and CS can be zero.
 559                IF (SN.NE.ZERO) THEN
 560                   IF (CS.NE.ZERO) THEN
 561                      ISUPPZ(2*M-1) = 1
 562                      ISUPPZ(2*M) = 2
 563                   ELSE
 564                      ISUPPZ(2*M-1) = 1
 565                      ISUPPZ(2*M) = 1
 566                   END IF
 567                ELSE
 568                   ISUPPZ(2*M-1) = 2
 569                   ISUPPZ(2*M) = 2
 570                END IF
 571             ENDIF
 572          ENDIF
 573       ELSE
 574
 575 *        Continue with general N
 576
 577          INDGRS = 1
 578          INDERR = 2*N + 1
 579          INDGP = 3*N + 1
 580          INDD = 4*N + 1
 581          INDE2 = 5*N + 1
 582          INDWRK = 6*N + 1
 583 *
 584          IINSPL = 1
 585          IINDBL = N + 1
 586          IINDW = 2*N + 1
 587          IINDWK = 3*N + 1
 588 *
 589 *        Scale matrix to allowable range, if necessary.
 590 *        The allowable range is related to the PIVMIN parameter; see the
 591 *        comments in DLARRD.  The preference for scaling small values
 592 *        up is heuristic; we expect users' matrices not to be close to the
 593 *        RMAX threshold.
 594 *
 595          SCALE = ONE
 596          TNRM = DLANST( 'M', N, D, E )
 597          IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
 598             SCALE = RMIN / TNRM
 599          ELSE IF( TNRM.GT.RMAX ) THEN
 600             SCALE = RMAX / TNRM
 601          END IF
 602          IF( SCALE.NE.ONE ) THEN
 603             CALL DSCAL( N, SCALE, D, 1 )
 604             CALL DSCAL( N-1, SCALE, E, 1 )
 605             TNRM = TNRM*SCALE
 606             IF( VALEIG ) THEN
 607 *              If eigenvalues in interval have to be found,
 608 *              scale (WL, WU] accordingly
 609                WL = WL*SCALE
 610                WU = WU*SCALE
 611             ENDIF
 612          END IF
 613 *
 614 *        Compute the desired eigenvalues of the tridiagonal after splitting
 615 *        into smaller subblocks if the corresponding off-diagonal elements
 616 *        are small
 617 *        THRESH is the splitting parameter for DLARRE
 618 *        A negative THRESH forces the old splitting criterion based on the
 619 *        size of the off-diagonal. A positive THRESH switches to splitting
 620 *        which preserves relative accuracy.
 621 *
 622          IF( TRYRAC ) THEN
 623 *           Test whether the matrix warrants the more expensive relative approach.
 624             CALL DLARRR( N, D, E, IINFO )
 625          ELSE
 626 *           The user does not care about relative accurately eigenvalues
 627             IINFO = -1
 628          ENDIF
 629 *        Set the splitting criterion
 630          IF (IINFO.EQ.0) THEN
 631             THRESH = EPS
 632          ELSE
 633             THRESH = -EPS
 634 *           relative accuracy is desired but T does not guarantee it
 635             TRYRAC = .FALSE.
 636          ENDIF
 637 *
 638          IF( TRYRAC ) THEN
 639 *           Copy original diagonal, needed to guarantee relative accuracy
 640             CALL DCOPY(N,D,1,WORK(INDD),1)
 641          ENDIF
 642 *        Store the squares of the offdiagonal values of T
 643          DO 5 J = 1, N-1
 644             WORK( INDE2+J-1 ) = E(J)**2
 645  5    CONTINUE
 646
 647 *        Set the tolerance parameters for bisection
 648          IF( .NOT.WANTZ ) THEN
 649 *           DLARRE computes the eigenvalues to full precision.
 650             RTOL1 = FOUR * EPS
 651             RTOL2 = FOUR * EPS
 652          ELSE
 653 *           DLARRE computes the eigenvalues to less than full precision.
 654 *           ZLARRV will refine the eigenvalue approximations, and we only
 655 *           need less accurate initial bisection in DLARRE.
 656 *           Note: these settings do only affect the subset case and DLARRE
 657             RTOL1 = SQRT(EPS)
 658             RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS )
 659          ENDIF
 660          CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E,
 661      $             WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT,
 662      $             IWORK( IINSPL ), M, W, WORK( INDERR ),
 663      $             WORK( INDGP ), IWORK( IINDBL ),
 664      $             IWORK( IINDW ), WORK( INDGRS ), PIVMIN,
 665      $             WORK( INDWRK ), IWORK( IINDWK ), IINFO )
 666          IF( IINFO.NE.0 ) THEN
 667             INFO = 10 + ABS( IINFO )
 668             RETURN
 669          END IF
 670 *        Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
 671 *        part of the spectrum. All desired eigenvalues are contained in
 672 *        (WL,WU]
 673
 674
 675          IF( WANTZ ) THEN
 676 *
 677 *           Compute the desired eigenvectors corresponding to the computed
 678 *           eigenvalues
 679 *
 680             CALL ZLARRV( N, WL, WU, D, E,
 681      $                PIVMIN, IWORK( IINSPL ), M,
 682      $                1, M, MINRGP, RTOL1, RTOL2,
 683      $                W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ),
 684      $                IWORK( IINDW ), WORK( INDGRS ), Z, LDZ,
 685      $                ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO )
 686             IF( IINFO.NE.0 ) THEN
 687                INFO = 20 + ABS( IINFO )
 688                RETURN
 689             END IF
 690          ELSE
 691 *           DLARRE computes eigenvalues of the (shifted) root representation
 692 *           ZLARRV returns the eigenvalues of the unshifted matrix.
 693 *           However, if the eigenvectors are not desired by the user, we need
 694 *           to apply the corresponding shifts from DLARRE to obtain the
 695 *           eigenvalues of the original matrix.
 696             DO 20 J = 1, M
 697                ITMP = IWORK( IINDBL+J-1 )
 698                W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) )
 699  20      CONTINUE
 700          END IF
 701 *
 702
 703          IF ( TRYRAC ) THEN
 704 *           Refine computed eigenvalues so that they are relatively accurate
 705 *           with respect to the original matrix T.
 706             IBEGIN = 1
 707             WBEGIN = 1
 708             DO 39  JBLK = 1, IWORK( IINDBL+M-1 )
 709                IEND = IWORK( IINSPL+JBLK-1 )
 710                IN = IEND - IBEGIN + 1
 711                WEND = WBEGIN - 1
 712 *              check if any eigenvalues have to be refined in this block
 713  36         CONTINUE
 714                IF( WEND.LT.M ) THEN
 715                   IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN
 716                      WEND = WEND + 1
 717                      GO TO 36
 718                   END IF
 719                END IF
 720                IF( WEND.LT.WBEGIN ) THEN
 721                   IBEGIN = IEND + 1
 722                   GO TO 39
 723                END IF
 724
 725                OFFSET = IWORK(IINDW+WBEGIN-1)-1
 726                IFIRST = IWORK(IINDW+WBEGIN-1)
 727                ILAST = IWORK(IINDW+WEND-1)
 728                RTOL2 = FOUR * EPS
 729                CALL DLARRJ( IN,
 730      $                   WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1),
 731      $                   IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN),
 732      $                   WORK( INDERR+WBEGIN-1 ),
 733      $                   WORK( INDWRK ), IWORK( IINDWK ), PIVMIN,
 734      $                   TNRM, IINFO )
 735                IBEGIN = IEND + 1
 736                WBEGIN = WEND + 1
 737  39      CONTINUE
 738          ENDIF
 739 *
 740 *        If matrix was scaled, then rescale eigenvalues appropriately.
 741 *
 742          IF( SCALE.NE.ONE ) THEN
 743             CALL DSCAL( M, ONE / SCALE, W, 1 )
 744          END IF
 745       END IF
 746 *
 747 *     If eigenvalues are not in increasing order, then sort them,
 748 *     possibly along with eigenvectors.
 749 *
 750       IF( NSPLIT.GT.1 .OR. N.EQ.2 ) THEN
 751          IF( .NOT. WANTZ ) THEN
 752             CALL DLASRT( 'I', M, W, IINFO )
 753             IF( IINFO.NE.0 ) THEN
 754                INFO = 3
 755                RETURN
 756             END IF
 757          ELSE
 758             DO 60 J = 1, M - 1
 759                I = 0
 760                TMP = W( J )
 761                DO 50 JJ = J + 1, M
 762                   IF( W( JJ ).LT.TMP ) THEN
 763                      I = JJ
 764                      TMP = W( JJ )
 765                   END IF
 766  50            CONTINUE
 767                IF( I.NE.0 ) THEN
 768                   W( I ) = W( J )
 769                   W( J ) = TMP
 770                   IF( WANTZ ) THEN
 771                      CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
 772                      ITMP = ISUPPZ( 2*I-1 )
 773                      ISUPPZ( 2*I-1 ) = ISUPPZ( 2*J-1 )
 774                      ISUPPZ( 2*J-1 ) = ITMP
 775                      ITMP = ISUPPZ( 2*I )
 776                      ISUPPZ( 2*I ) = ISUPPZ( 2*J )
 777                      ISUPPZ( 2*J ) = ITMP
 778                   END IF
 779                END IF
 780  60         CONTINUE
 781          END IF
 782       ENDIF
 783 *
 784 *
 785       WORK( 1 ) = LWMIN
 786       IWORK( 1 ) = LIWMIN
 787       RETURN
 788 *
 789 *     End of ZSTEMR
 790 *
 791       END