Annotation of rpl/lapack/lapack/dstedc.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
! 2: $ LIWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER COMPZ
! 11: INTEGER INFO, LDZ, LIWORK, LWORK, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: INTEGER IWORK( * )
! 15: DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
! 22: * symmetric tridiagonal matrix using the divide and conquer method.
! 23: * The eigenvectors of a full or band real symmetric matrix can also be
! 24: * found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
! 25: * matrix to tridiagonal form.
! 26: *
! 27: * This code makes very mild assumptions about floating point
! 28: * arithmetic. It will work on machines with a guard digit in
! 29: * add/subtract, or on those binary machines without guard digits
! 30: * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
! 31: * It could conceivably fail on hexadecimal or decimal machines
! 32: * without guard digits, but we know of none. See DLAED3 for details.
! 33: *
! 34: * Arguments
! 35: * =========
! 36: *
! 37: * COMPZ (input) CHARACTER*1
! 38: * = 'N': Compute eigenvalues only.
! 39: * = 'I': Compute eigenvectors of tridiagonal matrix also.
! 40: * = 'V': Compute eigenvectors of original dense symmetric
! 41: * matrix also. On entry, Z contains the orthogonal
! 42: * matrix used to reduce the original matrix to
! 43: * tridiagonal form.
! 44: *
! 45: * N (input) INTEGER
! 46: * The dimension of the symmetric tridiagonal matrix. N >= 0.
! 47: *
! 48: * D (input/output) DOUBLE PRECISION array, dimension (N)
! 49: * On entry, the diagonal elements of the tridiagonal matrix.
! 50: * On exit, if INFO = 0, the eigenvalues in ascending order.
! 51: *
! 52: * E (input/output) DOUBLE PRECISION array, dimension (N-1)
! 53: * On entry, the subdiagonal elements of the tridiagonal matrix.
! 54: * On exit, E has been destroyed.
! 55: *
! 56: * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
! 57: * On entry, if COMPZ = 'V', then Z contains the orthogonal
! 58: * matrix used in the reduction to tridiagonal form.
! 59: * On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
! 60: * orthonormal eigenvectors of the original symmetric matrix,
! 61: * and if COMPZ = 'I', Z contains the orthonormal eigenvectors
! 62: * of the symmetric tridiagonal matrix.
! 63: * If COMPZ = 'N', then Z is not referenced.
! 64: *
! 65: * LDZ (input) INTEGER
! 66: * The leading dimension of the array Z. LDZ >= 1.
! 67: * If eigenvectors are desired, then LDZ >= max(1,N).
! 68: *
! 69: * WORK (workspace/output) DOUBLE PRECISION array,
! 70: * dimension (LWORK)
! 71: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 72: *
! 73: * LWORK (input) INTEGER
! 74: * The dimension of the array WORK.
! 75: * If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
! 76: * If COMPZ = 'V' and N > 1 then LWORK must be at least
! 77: * ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
! 78: * where lg( N ) = smallest integer k such
! 79: * that 2**k >= N.
! 80: * If COMPZ = 'I' and N > 1 then LWORK must be at least
! 81: * ( 1 + 4*N + N**2 ).
! 82: * Note that for COMPZ = 'I' or 'V', then if N is less than or
! 83: * equal to the minimum divide size, usually 25, then LWORK need
! 84: * only be max(1,2*(N-1)).
! 85: *
! 86: * If LWORK = -1, then a workspace query is assumed; the routine
! 87: * only calculates the optimal size of the WORK array, returns
! 88: * this value as the first entry of the WORK array, and no error
! 89: * message related to LWORK is issued by XERBLA.
! 90: *
! 91: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
! 92: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
! 93: *
! 94: * LIWORK (input) INTEGER
! 95: * The dimension of the array IWORK.
! 96: * If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
! 97: * If COMPZ = 'V' and N > 1 then LIWORK must be at least
! 98: * ( 6 + 6*N + 5*N*lg N ).
! 99: * If COMPZ = 'I' and N > 1 then LIWORK must be at least
! 100: * ( 3 + 5*N ).
! 101: * Note that for COMPZ = 'I' or 'V', then if N is less than or
! 102: * equal to the minimum divide size, usually 25, then LIWORK
! 103: * need only be 1.
! 104: *
! 105: * If LIWORK = -1, then a workspace query is assumed; the
! 106: * routine only calculates the optimal size of the IWORK array,
! 107: * returns this value as the first entry of the IWORK array, and
! 108: * no error message related to LIWORK is issued by XERBLA.
! 109: *
! 110: * INFO (output) INTEGER
! 111: * = 0: successful exit.
! 112: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 113: * > 0: The algorithm failed to compute an eigenvalue while
! 114: * working on the submatrix lying in rows and columns
! 115: * INFO/(N+1) through mod(INFO,N+1).
! 116: *
! 117: * Further Details
! 118: * ===============
! 119: *
! 120: * Based on contributions by
! 121: * Jeff Rutter, Computer Science Division, University of California
! 122: * at Berkeley, USA
! 123: * Modified by Francoise Tisseur, University of Tennessee.
! 124: *
! 125: * =====================================================================
! 126: *
! 127: * .. Parameters ..
! 128: DOUBLE PRECISION ZERO, ONE, TWO
! 129: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
! 130: * ..
! 131: * .. Local Scalars ..
! 132: LOGICAL LQUERY
! 133: INTEGER FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
! 134: $ LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
! 135: DOUBLE PRECISION EPS, ORGNRM, P, TINY
! 136: * ..
! 137: * .. External Functions ..
! 138: LOGICAL LSAME
! 139: INTEGER ILAENV
! 140: DOUBLE PRECISION DLAMCH, DLANST
! 141: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
! 142: * ..
! 143: * .. External Subroutines ..
! 144: EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
! 145: $ DSTEQR, DSTERF, DSWAP, XERBLA
! 146: * ..
! 147: * .. Intrinsic Functions ..
! 148: INTRINSIC ABS, DBLE, INT, LOG, MAX, MOD, SQRT
! 149: * ..
! 150: * .. Executable Statements ..
! 151: *
! 152: * Test the input parameters.
! 153: *
! 154: INFO = 0
! 155: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
! 156: *
! 157: IF( LSAME( COMPZ, 'N' ) ) THEN
! 158: ICOMPZ = 0
! 159: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
! 160: ICOMPZ = 1
! 161: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
! 162: ICOMPZ = 2
! 163: ELSE
! 164: ICOMPZ = -1
! 165: END IF
! 166: IF( ICOMPZ.LT.0 ) THEN
! 167: INFO = -1
! 168: ELSE IF( N.LT.0 ) THEN
! 169: INFO = -2
! 170: ELSE IF( ( LDZ.LT.1 ) .OR.
! 171: $ ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
! 172: INFO = -6
! 173: END IF
! 174: *
! 175: IF( INFO.EQ.0 ) THEN
! 176: *
! 177: * Compute the workspace requirements
! 178: *
! 179: SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
! 180: IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
! 181: LIWMIN = 1
! 182: LWMIN = 1
! 183: ELSE IF( N.LE.SMLSIZ ) THEN
! 184: LIWMIN = 1
! 185: LWMIN = 2*( N - 1 )
! 186: ELSE
! 187: LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
! 188: IF( 2**LGN.LT.N )
! 189: $ LGN = LGN + 1
! 190: IF( 2**LGN.LT.N )
! 191: $ LGN = LGN + 1
! 192: IF( ICOMPZ.EQ.1 ) THEN
! 193: LWMIN = 1 + 3*N + 2*N*LGN + 3*N**2
! 194: LIWMIN = 6 + 6*N + 5*N*LGN
! 195: ELSE IF( ICOMPZ.EQ.2 ) THEN
! 196: LWMIN = 1 + 4*N + N**2
! 197: LIWMIN = 3 + 5*N
! 198: END IF
! 199: END IF
! 200: WORK( 1 ) = LWMIN
! 201: IWORK( 1 ) = LIWMIN
! 202: *
! 203: IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
! 204: INFO = -8
! 205: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
! 206: INFO = -10
! 207: END IF
! 208: END IF
! 209: *
! 210: IF( INFO.NE.0 ) THEN
! 211: CALL XERBLA( 'DSTEDC', -INFO )
! 212: RETURN
! 213: ELSE IF (LQUERY) THEN
! 214: RETURN
! 215: END IF
! 216: *
! 217: * Quick return if possible
! 218: *
! 219: IF( N.EQ.0 )
! 220: $ RETURN
! 221: IF( N.EQ.1 ) THEN
! 222: IF( ICOMPZ.NE.0 )
! 223: $ Z( 1, 1 ) = ONE
! 224: RETURN
! 225: END IF
! 226: *
! 227: * If the following conditional clause is removed, then the routine
! 228: * will use the Divide and Conquer routine to compute only the
! 229: * eigenvalues, which requires (3N + 3N**2) real workspace and
! 230: * (2 + 5N + 2N lg(N)) integer workspace.
! 231: * Since on many architectures DSTERF is much faster than any other
! 232: * algorithm for finding eigenvalues only, it is used here
! 233: * as the default. If the conditional clause is removed, then
! 234: * information on the size of workspace needs to be changed.
! 235: *
! 236: * If COMPZ = 'N', use DSTERF to compute the eigenvalues.
! 237: *
! 238: IF( ICOMPZ.EQ.0 ) THEN
! 239: CALL DSTERF( N, D, E, INFO )
! 240: GO TO 50
! 241: END IF
! 242: *
! 243: * If N is smaller than the minimum divide size (SMLSIZ+1), then
! 244: * solve the problem with another solver.
! 245: *
! 246: IF( N.LE.SMLSIZ ) THEN
! 247: *
! 248: CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
! 249: *
! 250: ELSE
! 251: *
! 252: * If COMPZ = 'V', the Z matrix must be stored elsewhere for later
! 253: * use.
! 254: *
! 255: IF( ICOMPZ.EQ.1 ) THEN
! 256: STOREZ = 1 + N*N
! 257: ELSE
! 258: STOREZ = 1
! 259: END IF
! 260: *
! 261: IF( ICOMPZ.EQ.2 ) THEN
! 262: CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
! 263: END IF
! 264: *
! 265: * Scale.
! 266: *
! 267: ORGNRM = DLANST( 'M', N, D, E )
! 268: IF( ORGNRM.EQ.ZERO )
! 269: $ GO TO 50
! 270: *
! 271: EPS = DLAMCH( 'Epsilon' )
! 272: *
! 273: START = 1
! 274: *
! 275: * while ( START <= N )
! 276: *
! 277: 10 CONTINUE
! 278: IF( START.LE.N ) THEN
! 279: *
! 280: * Let FINISH be the position of the next subdiagonal entry
! 281: * such that E( FINISH ) <= TINY or FINISH = N if no such
! 282: * subdiagonal exists. The matrix identified by the elements
! 283: * between START and FINISH constitutes an independent
! 284: * sub-problem.
! 285: *
! 286: FINISH = START
! 287: 20 CONTINUE
! 288: IF( FINISH.LT.N ) THEN
! 289: TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
! 290: $ SQRT( ABS( D( FINISH+1 ) ) )
! 291: IF( ABS( E( FINISH ) ).GT.TINY ) THEN
! 292: FINISH = FINISH + 1
! 293: GO TO 20
! 294: END IF
! 295: END IF
! 296: *
! 297: * (Sub) Problem determined. Compute its size and solve it.
! 298: *
! 299: M = FINISH - START + 1
! 300: IF( M.EQ.1 ) THEN
! 301: START = FINISH + 1
! 302: GO TO 10
! 303: END IF
! 304: IF( M.GT.SMLSIZ ) THEN
! 305: *
! 306: * Scale.
! 307: *
! 308: ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
! 309: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
! 310: $ INFO )
! 311: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
! 312: $ M-1, INFO )
! 313: *
! 314: IF( ICOMPZ.EQ.1 ) THEN
! 315: STRTRW = 1
! 316: ELSE
! 317: STRTRW = START
! 318: END IF
! 319: CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
! 320: $ Z( STRTRW, START ), LDZ, WORK( 1 ), N,
! 321: $ WORK( STOREZ ), IWORK, INFO )
! 322: IF( INFO.NE.0 ) THEN
! 323: INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
! 324: $ MOD( INFO, ( M+1 ) ) + START - 1
! 325: GO TO 50
! 326: END IF
! 327: *
! 328: * Scale back.
! 329: *
! 330: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
! 331: $ INFO )
! 332: *
! 333: ELSE
! 334: IF( ICOMPZ.EQ.1 ) THEN
! 335: *
! 336: * Since QR won't update a Z matrix which is larger than
! 337: * the length of D, we must solve the sub-problem in a
! 338: * workspace and then multiply back into Z.
! 339: *
! 340: CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
! 341: $ WORK( M*M+1 ), INFO )
! 342: CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
! 343: $ WORK( STOREZ ), N )
! 344: CALL DGEMM( 'N', 'N', N, M, M, ONE,
! 345: $ WORK( STOREZ ), N, WORK, M, ZERO,
! 346: $ Z( 1, START ), LDZ )
! 347: ELSE IF( ICOMPZ.EQ.2 ) THEN
! 348: CALL DSTEQR( 'I', M, D( START ), E( START ),
! 349: $ Z( START, START ), LDZ, WORK, INFO )
! 350: ELSE
! 351: CALL DSTERF( M, D( START ), E( START ), INFO )
! 352: END IF
! 353: IF( INFO.NE.0 ) THEN
! 354: INFO = START*( N+1 ) + FINISH
! 355: GO TO 50
! 356: END IF
! 357: END IF
! 358: *
! 359: START = FINISH + 1
! 360: GO TO 10
! 361: END IF
! 362: *
! 363: * endwhile
! 364: *
! 365: * If the problem split any number of times, then the eigenvalues
! 366: * will not be properly ordered. Here we permute the eigenvalues
! 367: * (and the associated eigenvectors) into ascending order.
! 368: *
! 369: IF( M.NE.N ) THEN
! 370: IF( ICOMPZ.EQ.0 ) THEN
! 371: *
! 372: * Use Quick Sort
! 373: *
! 374: CALL DLASRT( 'I', N, D, INFO )
! 375: *
! 376: ELSE
! 377: *
! 378: * Use Selection Sort to minimize swaps of eigenvectors
! 379: *
! 380: DO 40 II = 2, N
! 381: I = II - 1
! 382: K = I
! 383: P = D( I )
! 384: DO 30 J = II, N
! 385: IF( D( J ).LT.P ) THEN
! 386: K = J
! 387: P = D( J )
! 388: END IF
! 389: 30 CONTINUE
! 390: IF( K.NE.I ) THEN
! 391: D( K ) = D( I )
! 392: D( I ) = P
! 393: CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
! 394: END IF
! 395: 40 CONTINUE
! 396: END IF
! 397: END IF
! 398: END IF
! 399: *
! 400: 50 CONTINUE
! 401: WORK( 1 ) = LWMIN
! 402: IWORK( 1 ) = LIWMIN
! 403: *
! 404: RETURN
! 405: *
! 406: * End of DSTEDC
! 407: *
! 408: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dstedc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack