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