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