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