![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZHETRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO ) 2: * 3: * -- LAPACK routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * .. Scalar Arguments .. 9: CHARACTER UPLO 10: INTEGER INFO, LDA, LWORK, N 11: * .. 12: * .. Array Arguments .. 13: DOUBLE PRECISION D( * ), E( * ) 14: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * ZHETRD reduces a complex Hermitian matrix A to real symmetric 21: * tridiagonal form T by a unitary similarity transformation: 22: * Q**H * A * Q = T. 23: * 24: * Arguments 25: * ========= 26: * 27: * UPLO (input) CHARACTER*1 28: * = 'U': Upper triangle of A is stored; 29: * = 'L': Lower triangle of A is stored. 30: * 31: * N (input) INTEGER 32: * The order of the matrix A. N >= 0. 33: * 34: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 35: * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 36: * N-by-N upper triangular part of A contains the upper 37: * triangular part of the matrix A, and the strictly lower 38: * triangular part of A is not referenced. If UPLO = 'L', the 39: * leading N-by-N lower triangular part of A contains the lower 40: * triangular part of the matrix A, and the strictly upper 41: * triangular part of A is not referenced. 42: * On exit, if UPLO = 'U', the diagonal and first superdiagonal 43: * of A are overwritten by the corresponding elements of the 44: * tridiagonal matrix T, and the elements above the first 45: * superdiagonal, with the array TAU, represent the unitary 46: * matrix Q as a product of elementary reflectors; if UPLO 47: * = 'L', the diagonal and first subdiagonal of A are over- 48: * written by the corresponding elements of the tridiagonal 49: * matrix T, and the elements below the first subdiagonal, with 50: * the array TAU, represent the unitary matrix Q as a product 51: * of elementary reflectors. See Further Details. 52: * 53: * LDA (input) INTEGER 54: * The leading dimension of the array A. LDA >= max(1,N). 55: * 56: * D (output) DOUBLE PRECISION array, dimension (N) 57: * The diagonal elements of the tridiagonal matrix T: 58: * D(i) = A(i,i). 59: * 60: * E (output) DOUBLE PRECISION array, dimension (N-1) 61: * The off-diagonal elements of the tridiagonal matrix T: 62: * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 63: * 64: * TAU (output) COMPLEX*16 array, dimension (N-1) 65: * The scalar factors of the elementary reflectors (see Further 66: * Details). 67: * 68: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 69: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 70: * 71: * LWORK (input) INTEGER 72: * The dimension of the array WORK. LWORK >= 1. 73: * For optimum performance LWORK >= N*NB, where NB is the 74: * optimal blocksize. 75: * 76: * If LWORK = -1, then a workspace query is assumed; the routine 77: * only calculates the optimal size of the WORK array, returns 78: * this value as the first entry of the WORK array, and no error 79: * message related to LWORK is issued by XERBLA. 80: * 81: * INFO (output) INTEGER 82: * = 0: successful exit 83: * < 0: if INFO = -i, the i-th argument had an illegal value 84: * 85: * Further Details 86: * =============== 87: * 88: * If UPLO = 'U', the matrix Q is represented as a product of elementary 89: * reflectors 90: * 91: * Q = H(n-1) . . . H(2) H(1). 92: * 93: * Each H(i) has the form 94: * 95: * H(i) = I - tau * v * v' 96: * 97: * where tau is a complex scalar, and v is a complex vector with 98: * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 99: * A(1:i-1,i+1), and tau in TAU(i). 100: * 101: * If UPLO = 'L', the matrix Q is represented as a product of elementary 102: * reflectors 103: * 104: * Q = H(1) H(2) . . . H(n-1). 105: * 106: * Each H(i) has the form 107: * 108: * H(i) = I - tau * v * v' 109: * 110: * where tau is a complex scalar, and v is a complex vector with 111: * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 112: * and tau in TAU(i). 113: * 114: * The contents of A on exit are illustrated by the following examples 115: * with n = 5: 116: * 117: * if UPLO = 'U': if UPLO = 'L': 118: * 119: * ( d e v2 v3 v4 ) ( d ) 120: * ( d e v3 v4 ) ( e d ) 121: * ( d e v4 ) ( v1 e d ) 122: * ( d e ) ( v1 v2 e d ) 123: * ( d ) ( v1 v2 v3 e d ) 124: * 125: * where d and e denote diagonal and off-diagonal elements of T, and vi 126: * denotes an element of the vector defining H(i). 127: * 128: * ===================================================================== 129: * 130: * .. Parameters .. 131: DOUBLE PRECISION ONE 132: PARAMETER ( ONE = 1.0D+0 ) 133: COMPLEX*16 CONE 134: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 135: * .. 136: * .. Local Scalars .. 137: LOGICAL LQUERY, UPPER 138: INTEGER I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB, 139: $ NBMIN, NX 140: * .. 141: * .. External Subroutines .. 142: EXTERNAL XERBLA, ZHER2K, ZHETD2, ZLATRD 143: * .. 144: * .. Intrinsic Functions .. 145: INTRINSIC MAX 146: * .. 147: * .. External Functions .. 148: LOGICAL LSAME 149: INTEGER ILAENV 150: EXTERNAL LSAME, ILAENV 151: * .. 152: * .. Executable Statements .. 153: * 154: * Test the input parameters 155: * 156: INFO = 0 157: UPPER = LSAME( UPLO, 'U' ) 158: LQUERY = ( LWORK.EQ.-1 ) 159: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 160: INFO = -1 161: ELSE IF( N.LT.0 ) THEN 162: INFO = -2 163: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 164: INFO = -4 165: ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN 166: INFO = -9 167: END IF 168: * 169: IF( INFO.EQ.0 ) THEN 170: * 171: * Determine the block size. 172: * 173: NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) 174: LWKOPT = N*NB 175: WORK( 1 ) = LWKOPT 176: END IF 177: * 178: IF( INFO.NE.0 ) THEN 179: CALL XERBLA( 'ZHETRD', -INFO ) 180: RETURN 181: ELSE IF( LQUERY ) THEN 182: RETURN 183: END IF 184: * 185: * Quick return if possible 186: * 187: IF( N.EQ.0 ) THEN 188: WORK( 1 ) = 1 189: RETURN 190: END IF 191: * 192: NX = N 193: IWS = 1 194: IF( NB.GT.1 .AND. NB.LT.N ) THEN 195: * 196: * Determine when to cross over from blocked to unblocked code 197: * (last block is always handled by unblocked code). 198: * 199: NX = MAX( NB, ILAENV( 3, 'ZHETRD', UPLO, N, -1, -1, -1 ) ) 200: IF( NX.LT.N ) THEN 201: * 202: * Determine if workspace is large enough for blocked code. 203: * 204: LDWORK = N 205: IWS = LDWORK*NB 206: IF( LWORK.LT.IWS ) THEN 207: * 208: * Not enough workspace to use optimal NB: determine the 209: * minimum value of NB, and reduce NB or force use of 210: * unblocked code by setting NX = N. 211: * 212: NB = MAX( LWORK / LDWORK, 1 ) 213: NBMIN = ILAENV( 2, 'ZHETRD', UPLO, N, -1, -1, -1 ) 214: IF( NB.LT.NBMIN ) 215: $ NX = N 216: END IF 217: ELSE 218: NX = N 219: END IF 220: ELSE 221: NB = 1 222: END IF 223: * 224: IF( UPPER ) THEN 225: * 226: * Reduce the upper triangle of A. 227: * Columns 1:kk are handled by the unblocked method. 228: * 229: KK = N - ( ( N-NX+NB-1 ) / NB )*NB 230: DO 20 I = N - NB + 1, KK + 1, -NB 231: * 232: * Reduce columns i:i+nb-1 to tridiagonal form and form the 233: * matrix W which is needed to update the unreduced part of 234: * the matrix 235: * 236: CALL ZLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK, 237: $ LDWORK ) 238: * 239: * Update the unreduced submatrix A(1:i-1,1:i-1), using an 240: * update of the form: A := A - V*W' - W*V' 241: * 242: CALL ZHER2K( UPLO, 'No transpose', I-1, NB, -CONE, 243: $ A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA ) 244: * 245: * Copy superdiagonal elements back into A, and diagonal 246: * elements into D 247: * 248: DO 10 J = I, I + NB - 1 249: A( J-1, J ) = E( J-1 ) 250: D( J ) = A( J, J ) 251: 10 CONTINUE 252: 20 CONTINUE 253: * 254: * Use unblocked code to reduce the last or only block 255: * 256: CALL ZHETD2( UPLO, KK, A, LDA, D, E, TAU, IINFO ) 257: ELSE 258: * 259: * Reduce the lower triangle of A 260: * 261: DO 40 I = 1, N - NX, NB 262: * 263: * Reduce columns i:i+nb-1 to tridiagonal form and form the 264: * matrix W which is needed to update the unreduced part of 265: * the matrix 266: * 267: CALL ZLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ), 268: $ TAU( I ), WORK, LDWORK ) 269: * 270: * Update the unreduced submatrix A(i+nb:n,i+nb:n), using 271: * an update of the form: A := A - V*W' - W*V' 272: * 273: CALL ZHER2K( UPLO, 'No transpose', N-I-NB+1, NB, -CONE, 274: $ A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE, 275: $ A( I+NB, I+NB ), LDA ) 276: * 277: * Copy subdiagonal elements back into A, and diagonal 278: * elements into D 279: * 280: DO 30 J = I, I + NB - 1 281: A( J+1, J ) = E( J ) 282: D( J ) = A( J, J ) 283: 30 CONTINUE 284: 40 CONTINUE 285: * 286: * Use unblocked code to reduce the last or only block 287: * 288: CALL ZHETD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ), 289: $ TAU( I ), IINFO ) 290: END IF 291: * 292: WORK( 1 ) = LWKOPT 293: RETURN 294: * 295: * End of ZHETRD 296: * 297: END