zstemr.f

Go to the documentation of this file.

*> \brief \b ZSTEMR
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> Download ZSTEMR + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f">
*> [TXT]</a>
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
*                          M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
*                          IWORK, LIWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          JOBZ, RANGE
*       LOGICAL            TRYRAC
*       INTEGER            IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
*       DOUBLE PRECISION VL, VU
*       ..
*       .. Array Arguments ..
*       INTEGER            ISUPPZ( * ), IWORK( * )
*       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
*       COMPLEX*16         Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZSTEMR computes selected eigenvalues and, optionally, eigenvectors
*> of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
*> a well defined set of pairwise different real eigenvalues, the corresponding
*> real eigenvectors are pairwise orthogonal.
*>
*> The spectrum may be computed either completely or partially by specifying
*> either an interval (VL,VU] or a range of indices IL:IU for the desired
*> eigenvalues.
*>
*> Depending on the number of desired eigenvalues, these are computed either
*> by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are
*> computed by the use of various suitable L D L^T factorizations near clusters
*> of close eigenvalues (referred to as RRRs, Relatively Robust
*> Representations). An informal sketch of the algorithm follows.
*>
*> For each unreduced block (submatrix) of T,
*>    (a) Compute T - sigma I  = L D L^T, so that L and D
*>        define all the wanted eigenvalues to high relative accuracy.
*>        This means that small relative changes in the entries of D and L
*>        cause only small relative changes in the eigenvalues and
*>        eigenvectors. The standard (unfactored) representation of the
*>        tridiagonal matrix T does not have this property in general.
*>    (b) Compute the eigenvalues to suitable accuracy.
*>        If the eigenvectors are desired, the algorithm attains full
*>        accuracy of the computed eigenvalues only right before
*>        the corresponding vectors have to be computed, see steps c) and d).
*>    (c) For each cluster of close eigenvalues, select a new
*>        shift close to the cluster, find a new factorization, and refine
*>        the shifted eigenvalues to suitable accuracy.
*>    (d) For each eigenvalue with a large enough relative separation compute
*>        the corresponding eigenvector by forming a rank revealing twisted
*>        factorization. Go back to (c) for any clusters that remain.
*>
*> For more details, see:
*> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
*>   to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
*>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
*> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
*>   Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
*>   2004.  Also LAPACK Working Note 154.
*> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
*>   tridiagonal eigenvalue/eigenvector problem",
*>   Computer Science Division Technical Report No. UCB/CSD-97-971,
*>   UC Berkeley, May 1997.
*>
*> Further Details
*> 1.ZSTEMR works only on machines which follow IEEE-754
*> floating-point standard in their handling of infinities and NaNs.
*> This permits the use of efficient inner loops avoiding a check for
*> zero divisors.
*>
*> 2. LAPACK routines can be used to reduce a complex Hermitean matrix to
*> real symmetric tridiagonal form.
*>
*> (Any complex Hermitean tridiagonal matrix has real values on its diagonal
*> and potentially complex numbers on its off-diagonals. By applying a
*> similarity transform with an appropriate diagonal matrix
*> diag(1,e^{i \phy_1}, ... , e^{i \phy_{n-1}}), the complex Hermitean
*> matrix can be transformed into a real symmetric matrix and complex
*> arithmetic can be entirely avoided.)
*>
*> While the eigenvectors of the real symmetric tridiagonal matrix are real,
*> the eigenvectors of original complex Hermitean matrix have complex entries
*> in general.
*> Since LAPACK drivers overwrite the matrix data with the eigenvectors,
*> ZSTEMR accepts complex workspace to facilitate interoperability
*> with ZUNMTR or ZUPMTR.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOBZ
*> \verbatim
*>          JOBZ is CHARACTER*1
*>          = 'N':  Compute eigenvalues only;
*>          = 'V':  Compute eigenvalues and eigenvectors.
*> \endverbatim
*>
*> \param[in] RANGE
*> \verbatim
*>          RANGE is CHARACTER*1
*>          = 'A': all eigenvalues will be found.
*>          = 'V': all eigenvalues in the half-open interval (VL,VU]
*>                 will be found.
*>          = 'I': the IL-th through IU-th eigenvalues will be found.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] D
*> \verbatim
*>          D is DOUBLE PRECISION array, dimension (N)
*>          On entry, the N diagonal elements of the tridiagonal matrix
*>          T. On exit, D is overwritten.
*> \endverbatim
*>
*> \param[in,out] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (N)
*>          On entry, the (N-1) subdiagonal elements of the tridiagonal
*>          matrix T in elements 1 to N-1 of E. E(N) need not be set on
*>          input, but is used internally as workspace.
*>          On exit, E is overwritten.
*> \endverbatim
*>
*> \param[in] VL
*> \verbatim
*>          VL is DOUBLE PRECISION
*>
*>          If RANGE='V', the lower bound of the interval to
*>          be searched for eigenvalues. VL < VU.
*>          Not referenced if RANGE = 'A' or 'I'.
*> \endverbatim
*>
*> \param[in] VU
*> \verbatim
*>          VU is DOUBLE PRECISION
*>
*>          If RANGE='V', the upper bound of the interval to
*>          be searched for eigenvalues. VL < VU.
*>          Not referenced if RANGE = 'A' or 'I'.
*> \endverbatim
*>
*> \param[in] IL
*> \verbatim
*>          IL is INTEGER
*>
*>          If RANGE='I', the index of the
*>          smallest eigenvalue to be returned.
*>          1 <= IL <= IU <= N, if N > 0.
*>          Not referenced if RANGE = 'A' or 'V'.
*> \endverbatim
*>
*> \param[in] IU
*> \verbatim
*>          IU is INTEGER
*>
*>          If RANGE='I', the index of the
*>          largest eigenvalue to be returned.
*>          1 <= IL <= IU <= N, if N > 0.
*>          Not referenced if RANGE = 'A' or 'V'.
*> \endverbatim
*>
*> \param[out] M
*> \verbatim
*>          M is INTEGER
*>          The total number of eigenvalues found.  0 <= M <= N.
*>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
*> \endverbatim
*>
*> \param[out] W
*> \verbatim
*>          W is DOUBLE PRECISION array, dimension (N)
*>          The first M elements contain the selected eigenvalues in
*>          ascending order.
*> \endverbatim
*>
*> \param[out] Z
*> \verbatim
*>          Z is COMPLEX*16 array, dimension (LDZ, max(1,M) )
*>          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
*>          contain the orthonormal eigenvectors of the matrix T
*>          corresponding to the selected eigenvalues, with the i-th
*>          column of Z holding the eigenvector associated with W(i).
*>          If JOBZ = 'N', then Z is not referenced.
*>          Note: the user must ensure that at least max(1,M) columns are
*>          supplied in the array Z; if RANGE = 'V', the exact value of M
*>          is not known in advance and can be computed with a workspace
*>          query by setting NZC = -1, see below.
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>          The leading dimension of the array Z.  LDZ >= 1, and if
*>          JOBZ = 'V', then LDZ >= max(1,N).
*> \endverbatim
*>
*> \param[in] NZC
*> \verbatim
*>          NZC is INTEGER
*>          The number of eigenvectors to be held in the array Z.
*>          If RANGE = 'A', then NZC >= max(1,N).
*>          If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
*>          If RANGE = 'I', then NZC >= IU-IL+1.
*>          If NZC = -1, then a workspace query is assumed; the
*>          routine calculates the number of columns of the array Z that
*>          are needed to hold the eigenvectors.
*>          This value is returned as the first entry of the Z array, and
*>          no error message related to NZC is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] ISUPPZ
*> \verbatim
*>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
*>          The support of the eigenvectors in Z, i.e., the indices
*>          indicating the nonzero elements in Z. The i-th computed eigenvector
*>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
*>          ISUPPZ( 2*i ). This is relevant in the case when the matrix
*>          is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
*> \endverbatim
*>
*> \param[in,out] TRYRAC
*> \verbatim
*>          TRYRAC is LOGICAL
*>          If TRYRAC = .TRUE., indicates that the code should check whether
*>          the tridiagonal matrix defines its eigenvalues to high relative
*>          accuracy.  If so, the code uses relative-accuracy preserving
*>          algorithms that might be (a bit) slower depending on the matrix.
*>          If the matrix does not define its eigenvalues to high relative
*>          accuracy, the code can uses possibly faster algorithms.
*>          If TRYRAC = .FALSE., the code is not required to guarantee
*>          relatively accurate eigenvalues and can use the fastest possible
*>          techniques.
*>          On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
*>          does not define its eigenvalues to high relative accuracy.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (LWORK)
*>          On exit, if INFO = 0, WORK(1) returns the optimal
*>          (and minimal) LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK. LWORK >= max(1,18*N)
*>          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates the optimal size of the WORK array, returns
*>          this value as the first entry of the WORK array, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (LIWORK)
*>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
*> \endverbatim
*>
*> \param[in] LIWORK
*> \verbatim
*>          LIWORK is INTEGER
*>          The dimension of the array IWORK.  LIWORK >= max(1,10*N)
*>          if the eigenvectors are desired, and LIWORK >= max(1,8*N)
*>          if only the eigenvalues are to be computed.
*>          If LIWORK = -1, then a workspace query is assumed; the
*>          routine only calculates the optimal size of the IWORK array,
*>          returns this value as the first entry of the IWORK array, and
*>          no error message related to LIWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          On exit, INFO
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*>          > 0:  if INFO = 1X, internal error in DLARRE,
*>                if INFO = 2X, internal error in ZLARRV.
*>                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
*>                the nonzero error code returned by DLARRE or
*>                ZLARRV, respectively.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup stemr
*
*> \par Contributors:
*  ==================
*>
*> Beresford Parlett, University of California, Berkeley, USA \n
*> Jim Demmel, University of California, Berkeley, USA \n
*> Inderjit Dhillon, University of Texas, Austin, USA \n
*> Osni Marques, LBNL/NERSC, USA \n
*> Christof Voemel, University of California, Berkeley, USA \n
*> Aravindh Krishnamoorthy, FAU, Erlangen, Germany \n
*
*  =====================================================================
      SUBROUTINE zstemr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
     $                   M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
     $                   IWORK, LIWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ, RANGE
      LOGICAL            TRYRAC
      INTEGER            IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
      DOUBLE PRECISION VL, VU
*     ..
*     .. Array Arguments ..
      INTEGER            ISUPPZ( * ), IWORK( * )
      DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
      COMPLEX*16         Z( LDZ, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, FOUR, MINRGP
      PARAMETER          ( ZERO = 0.0d0, one = 1.0d0,
     $                     four = 4.0d0,
     $                     minrgp = 1.0d-3 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
     $                   LAESWAP
      INTEGER            I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
     $                   IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
     $                   inde2, inderr, indgp, indgrs, indwrk, itmp,
     $                   itmp2, j, jblk, jj, liwmin, lwmin, nsplit,
     $                   nzcmin, offset, wbegin, wend
      DOUBLE PRECISION   BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN,
     $                   RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN,
     $                   thresh, tmp, tnrm, wl, wu
*     ..
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLANST
      EXTERNAL           lsame, dlamch, dlanst
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlae2, dlaev2, dlarrc, dlarre,
     $                   dlarrj,
     $                   dlarrr, dlasrt, dscal, xerbla, zlarrv, zswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min, sqrt
 
 
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      wantz = lsame( jobz, 'V' )
      alleig = lsame( range, 'A' )
      valeig = lsame( range, 'V' )
      indeig = lsame( range, 'I' )
*
      lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
      zquery = ( nzc.EQ.-1 )
      laeswap = .false.
 
*     DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
*     In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
*     Furthermore, ZLARRV needs WORK of size 12*N, IWORK of size 7*N.
      IF( wantz ) THEN
         lwmin = 18*n
         liwmin = 10*n
      ELSE
*        need less workspace if only the eigenvalues are wanted
         lwmin = 12*n
         liwmin = 8*n
      ENDIF
 
      wl = zero
      wu = zero
      iil = 0
      iiu = 0
      nsplit = 0
 
      IF( valeig ) THEN
*        We do not reference VL, VU in the cases RANGE = 'I','A'
*        The interval (WL, WU] contains all the wanted eigenvalues.
*        It is either given by the user or computed in DLARRE.
         wl = vl
         wu = vu
      ELSEIF( indeig ) THEN
*        We do not reference IL, IU in the cases RANGE = 'V','A'
         iil = il
         iiu = iu
      ENDIF
*
      info = 0
      IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
         info = -1
      ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
         info = -2
      ELSE IF( n.LT.0 ) THEN
         info = -3
      ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
         info = -7
      ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
         info = -8
      ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
         info = -9
      ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
         info = -13
      ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
         info = -17
      ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
         info = -19
      END IF
*
*     Get machine constants.
*
      safmin = dlamch( 'Safe minimum' )
      eps = dlamch( 'Precision' )
      smlnum = safmin / eps
      bignum = one / smlnum
      rmin = sqrt( smlnum )
      rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
*
      IF( info.EQ.0 ) THEN
         work( 1 ) = lwmin
         iwork( 1 ) = liwmin
*
         IF( wantz .AND. alleig ) THEN
            nzcmin = n
         ELSE IF( wantz .AND. valeig ) THEN
            CALL dlarrc( 'T', n, vl, vu, d, e, safmin,
     $                            nzcmin, itmp, itmp2, info )
         ELSE IF( wantz .AND. indeig ) THEN
            nzcmin = iiu-iil+1
         ELSE
*           WANTZ .EQ. FALSE.
            nzcmin = 0
         ENDIF
         IF( zquery .AND. info.EQ.0 ) THEN
            z( 1,1 ) = nzcmin
         ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
            info = -14
         END IF
      END IF
 
      IF( info.NE.0 ) THEN
*
         CALL xerbla( 'ZSTEMR', -info )
*
         RETURN
      ELSE IF( lquery .OR. zquery ) THEN
         RETURN
      END IF
*
*     Handle N = 0, 1, and 2 cases immediately
*
      m = 0
      IF( n.EQ.0 )
     $   RETURN
*
      IF( n.EQ.1 ) THEN
         IF( alleig .OR. indeig ) THEN
            m = 1
            w( 1 ) = d( 1 )
         ELSE
            IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
               m = 1
               w( 1 ) = d( 1 )
            END IF
         END IF
         IF( wantz.AND.(.NOT.zquery) ) THEN
            z( 1, 1 ) = one
            isuppz(1) = 1
            isuppz(2) = 1
         END IF
         RETURN
      END IF
*
      IF( n.EQ.2 ) THEN
         IF( .NOT.wantz ) THEN
            CALL dlae2( d(1), e(1), d(2), r1, r2 )
         ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
            CALL dlaev2( d(1), e(1), d(2), r1, r2, cs, sn )
         END IF
*        D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
*        the following code requires R1 >= R2. Hence, we correct
*        the order of R1, R2, CS, SN if R1 < R2 before further processing.
         IF( r1.LT.r2 ) THEN
            e(2) = r1
            r1 = r2
            r2 = e(2)
            laeswap = .true.
         ENDIF
         IF( alleig.OR.
     $      (valeig.AND.(r2.GT.wl).AND.
     $                  (r2.LE.wu)).OR.
     $      (indeig.AND.(iil.EQ.1)) ) THEN
            m = m+1
            w( m ) = r2
            IF( wantz.AND.(.NOT.zquery) ) THEN
               IF( laeswap ) THEN
                  z( 1, m ) = cs
                  z( 2, m ) = sn
               ELSE
                  z( 1, m ) = -sn
                  z( 2, m ) = cs
               ENDIF
*              Note: At most one of SN and CS can be zero.
               IF (sn.NE.zero) THEN
                  IF (cs.NE.zero) THEN
                     isuppz(2*m-1) = 1
                     isuppz(2*m) = 2
                  ELSE
                     isuppz(2*m-1) = 1
                     isuppz(2*m) = 1
                  END IF
               ELSE
                  isuppz(2*m-1) = 2
                  isuppz(2*m) = 2
               END IF
            ENDIF
         ENDIF
         IF( alleig.OR.
     $      (valeig.AND.(r1.GT.wl).AND.
     $                  (r1.LE.wu)).OR.
     $      (indeig.AND.(iiu.EQ.2)) ) THEN
            m = m+1
            w( m ) = r1
            IF( wantz.AND.(.NOT.zquery) ) THEN
               IF( laeswap ) THEN
                  z( 1, m ) = -sn
                  z( 2, m ) = cs
               ELSE
                  z( 1, m ) = cs
                  z( 2, m ) = sn
               ENDIF
*              Note: At most one of SN and CS can be zero.
               IF (sn.NE.zero) THEN
                  IF (cs.NE.zero) THEN
                     isuppz(2*m-1) = 1
                     isuppz(2*m) = 2
                  ELSE
                     isuppz(2*m-1) = 1
                     isuppz(2*m) = 1
                  END IF
               ELSE
                  isuppz(2*m-1) = 2
                  isuppz(2*m) = 2
               END IF
            ENDIF
         ENDIF
      ELSE
 
*        Continue with general N
 
         indgrs = 1
         inderr = 2*n + 1
         indgp = 3*n + 1
         indd = 4*n + 1
         inde2 = 5*n + 1
         indwrk = 6*n + 1
*
         iinspl = 1
         iindbl = n + 1
         iindw = 2*n + 1
         iindwk = 3*n + 1
*
*        Scale matrix to allowable range, if necessary.
*        The allowable range is related to the PIVMIN parameter; see the
*        comments in DLARRD.  The preference for scaling small values
*        up is heuristic; we expect users' matrices not to be close to the
*        RMAX threshold.
*
         scale = one
         tnrm = dlanst( 'M', n, d, e )
         IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
            scale = rmin / tnrm
         ELSE IF( tnrm.GT.rmax ) THEN
            scale = rmax / tnrm
         END IF
         IF( scale.NE.one ) THEN
            CALL dscal( n, scale, d, 1 )
            CALL dscal( n-1, scale, e, 1 )
            tnrm = tnrm*scale
            IF( valeig ) THEN
*              If eigenvalues in interval have to be found,
*              scale (WL, WU] accordingly
               wl = wl*scale
               wu = wu*scale
            ENDIF
         END IF
*
*        Compute the desired eigenvalues of the tridiagonal after splitting
*        into smaller subblocks if the corresponding off-diagonal elements
*        are small
*        THRESH is the splitting parameter for DLARRE
*        A negative THRESH forces the old splitting criterion based on the
*        size of the off-diagonal. A positive THRESH switches to splitting
*        which preserves relative accuracy.
*
         IF( tryrac ) THEN
*           Test whether the matrix warrants the more expensive relative approach.
            CALL dlarrr( n, d, e, iinfo )
         ELSE
*           The user does not care about relative accurately eigenvalues
            iinfo = -1
         ENDIF
*        Set the splitting criterion
         IF (iinfo.EQ.0) THEN
            thresh = eps
         ELSE
            thresh = -eps
*           relative accuracy is desired but T does not guarantee it
            tryrac = .false.
         ENDIF
*
         IF( tryrac ) THEN
*           Copy original diagonal, needed to guarantee relative accuracy
            CALL dcopy(n,d,1,work(indd),1)
         ENDIF
*        Store the squares of the offdiagonal values of T
         DO 5 j = 1, n-1
            work( inde2+j-1 ) = e(j)**2
 5    CONTINUE
 
*        Set the tolerance parameters for bisection
         IF( .NOT.wantz ) THEN
*           DLARRE computes the eigenvalues to full precision.
            rtol1 = four * eps
            rtol2 = four * eps
         ELSE
*           DLARRE computes the eigenvalues to less than full precision.
*           ZLARRV will refine the eigenvalue approximations, and we only
*           need less accurate initial bisection in DLARRE.
*           Note: these settings do only affect the subset case and DLARRE
            rtol1 = sqrt(eps)
            rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
         ENDIF
         CALL dlarre( range, n, wl, wu, iil, iiu, d, e,
     $             work(inde2), rtol1, rtol2, thresh, nsplit,
     $             iwork( iinspl ), m, w, work( inderr ),
     $             work( indgp ), iwork( iindbl ),
     $             iwork( iindw ), work( indgrs ), pivmin,
     $             work( indwrk ), iwork( iindwk ), iinfo )
         IF( iinfo.NE.0 ) THEN
            info = 10 + abs( iinfo )
            RETURN
         END IF
*        Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
*        part of the spectrum. All desired eigenvalues are contained in
*        (WL,WU]
 
 
         IF( wantz ) THEN
*
*           Compute the desired eigenvectors corresponding to the computed
*           eigenvalues
*
            CALL zlarrv( n, wl, wu, d, e,
     $                pivmin, iwork( iinspl ), m,
     $                1, m, minrgp, rtol1, rtol2,
     $                w, work( inderr ), work( indgp ), iwork( iindbl ),
     $                iwork( iindw ), work( indgrs ), z, ldz,
     $                isuppz, work( indwrk ), iwork( iindwk ), iinfo )
            IF( iinfo.NE.0 ) THEN
               info = 20 + abs( iinfo )
               RETURN
            END IF
         ELSE
*           DLARRE computes eigenvalues of the (shifted) root representation
*           ZLARRV returns the eigenvalues of the unshifted matrix.
*           However, if the eigenvectors are not desired by the user, we need
*           to apply the corresponding shifts from DLARRE to obtain the
*           eigenvalues of the original matrix.
            DO 20 j = 1, m
               itmp = iwork( iindbl+j-1 )
               w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
 20      CONTINUE
         END IF
*
 
         IF ( tryrac ) THEN
*           Refine computed eigenvalues so that they are relatively accurate
*           with respect to the original matrix T.
            ibegin = 1
            wbegin = 1
            DO 39  jblk = 1, iwork( iindbl+m-1 )
               iend = iwork( iinspl+jblk-1 )
               in = iend - ibegin + 1
               wend = wbegin - 1
*              check if any eigenvalues have to be refined in this block
 36         CONTINUE
               IF( wend.LT.m ) THEN
                  IF( iwork( iindbl+wend ).EQ.jblk ) THEN
                     wend = wend + 1
                     GO TO 36
                  END IF
               END IF
               IF( wend.LT.wbegin ) THEN
                  ibegin = iend + 1
                  GO TO 39
               END IF
 
               offset = iwork(iindw+wbegin-1)-1
               ifirst = iwork(iindw+wbegin-1)
               ilast = iwork(iindw+wend-1)
               rtol2 = four * eps
               CALL dlarrj( in,
     $                   work(indd+ibegin-1), work(inde2+ibegin-1),
     $                   ifirst, ilast, rtol2, offset, w(wbegin),
     $                   work( inderr+wbegin-1 ),
     $                   work( indwrk ), iwork( iindwk ), pivmin,
     $                   tnrm, iinfo )
               ibegin = iend + 1
               wbegin = wend + 1
 39      CONTINUE
         ENDIF
*
*        If matrix was scaled, then rescale eigenvalues appropriately.
*
         IF( scale.NE.one ) THEN
            CALL dscal( m, one / scale, w, 1 )
         END IF
      END IF
*
*     If eigenvalues are not in increasing order, then sort them,
*     possibly along with eigenvectors.
*
      IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
         IF( .NOT. wantz ) THEN
            CALL dlasrt( 'I', m, w, iinfo )
            IF( iinfo.NE.0 ) THEN
               info = 3
               RETURN
            END IF
         ELSE
            DO 60 j = 1, m - 1
               i = 0
               tmp = w( j )
               DO 50 jj = j + 1, m
                  IF( w( jj ).LT.tmp ) THEN
                     i = jj
                     tmp = w( jj )
                  END IF
 50            CONTINUE
               IF( i.NE.0 ) THEN
                  w( i ) = w( j )
                  w( j ) = tmp
                  IF( wantz ) THEN
                     CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
                     itmp = isuppz( 2*i-1 )
                     isuppz( 2*i-1 ) = isuppz( 2*j-1 )
                     isuppz( 2*j-1 ) = itmp
                     itmp = isuppz( 2*i )
                     isuppz( 2*i ) = isuppz( 2*j )
                     isuppz( 2*j ) = itmp
                  END IF
               END IF
 60         CONTINUE
         END IF
      ENDIF
*
*
      work( 1 ) = lwmin
      iwork( 1 ) = liwmin
      RETURN
*
*     End of ZSTEMR
*
      SUBROUTINE zstemr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, …
      END