![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZHETD2( UPLO, N, A, LDA, D, E, TAU, 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, N 11: * .. 12: * .. Array Arguments .. 13: DOUBLE PRECISION D( * ), E( * ) 14: COMPLEX*16 A( LDA, * ), TAU( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * ZHETD2 reduces a complex Hermitian matrix A to real symmetric 21: * tridiagonal form T by a unitary similarity transformation: 22: * Q' * A * Q = T. 23: * 24: * Arguments 25: * ========= 26: * 27: * UPLO (input) CHARACTER*1 28: * Specifies whether the upper or lower triangular part of the 29: * Hermitian matrix A is stored: 30: * = 'U': Upper triangular 31: * = 'L': Lower triangular 32: * 33: * N (input) INTEGER 34: * The order of the matrix A. N >= 0. 35: * 36: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 37: * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 38: * n-by-n upper triangular part of A contains the upper 39: * triangular part of the matrix A, and the strictly lower 40: * triangular part of A is not referenced. If UPLO = 'L', the 41: * leading n-by-n lower triangular part of A contains the lower 42: * triangular part of the matrix A, and the strictly upper 43: * triangular part of A is not referenced. 44: * On exit, if UPLO = 'U', the diagonal and first superdiagonal 45: * of A are overwritten by the corresponding elements of the 46: * tridiagonal matrix T, and the elements above the first 47: * superdiagonal, with the array TAU, represent the unitary 48: * matrix Q as a product of elementary reflectors; if UPLO 49: * = 'L', the diagonal and first subdiagonal of A are over- 50: * written by the corresponding elements of the tridiagonal 51: * matrix T, and the elements below the first subdiagonal, with 52: * the array TAU, represent the unitary matrix Q as a product 53: * of elementary reflectors. See Further Details. 54: * 55: * LDA (input) INTEGER 56: * The leading dimension of the array A. LDA >= max(1,N). 57: * 58: * D (output) DOUBLE PRECISION array, dimension (N) 59: * The diagonal elements of the tridiagonal matrix T: 60: * D(i) = A(i,i). 61: * 62: * E (output) DOUBLE PRECISION array, dimension (N-1) 63: * The off-diagonal elements of the tridiagonal matrix T: 64: * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 65: * 66: * TAU (output) COMPLEX*16 array, dimension (N-1) 67: * The scalar factors of the elementary reflectors (see Further 68: * Details). 69: * 70: * INFO (output) INTEGER 71: * = 0: successful exit 72: * < 0: if INFO = -i, the i-th argument had an illegal value. 73: * 74: * Further Details 75: * =============== 76: * 77: * If UPLO = 'U', the matrix Q is represented as a product of elementary 78: * reflectors 79: * 80: * Q = H(n-1) . . . H(2) H(1). 81: * 82: * Each H(i) has the form 83: * 84: * H(i) = I - tau * v * v' 85: * 86: * where tau is a complex scalar, and v is a complex vector with 87: * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 88: * A(1:i-1,i+1), and tau in TAU(i). 89: * 90: * If UPLO = 'L', the matrix Q is represented as a product of elementary 91: * reflectors 92: * 93: * Q = H(1) H(2) . . . H(n-1). 94: * 95: * Each H(i) has the form 96: * 97: * H(i) = I - tau * v * v' 98: * 99: * where tau is a complex scalar, and v is a complex vector with 100: * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 101: * and tau in TAU(i). 102: * 103: * The contents of A on exit are illustrated by the following examples 104: * with n = 5: 105: * 106: * if UPLO = 'U': if UPLO = 'L': 107: * 108: * ( d e v2 v3 v4 ) ( d ) 109: * ( d e v3 v4 ) ( e d ) 110: * ( d e v4 ) ( v1 e d ) 111: * ( d e ) ( v1 v2 e d ) 112: * ( d ) ( v1 v2 v3 e d ) 113: * 114: * where d and e denote diagonal and off-diagonal elements of T, and vi 115: * denotes an element of the vector defining H(i). 116: * 117: * ===================================================================== 118: * 119: * .. Parameters .. 120: COMPLEX*16 ONE, ZERO, HALF 121: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 122: $ ZERO = ( 0.0D+0, 0.0D+0 ), 123: $ HALF = ( 0.5D+0, 0.0D+0 ) ) 124: * .. 125: * .. Local Scalars .. 126: LOGICAL UPPER 127: INTEGER I 128: COMPLEX*16 ALPHA, TAUI 129: * .. 130: * .. External Subroutines .. 131: EXTERNAL XERBLA, ZAXPY, ZHEMV, ZHER2, ZLARFG 132: * .. 133: * .. External Functions .. 134: LOGICAL LSAME 135: COMPLEX*16 ZDOTC 136: EXTERNAL LSAME, ZDOTC 137: * .. 138: * .. Intrinsic Functions .. 139: INTRINSIC DBLE, MAX, MIN 140: * .. 141: * .. Executable Statements .. 142: * 143: * Test the input parameters 144: * 145: INFO = 0 146: UPPER = LSAME( UPLO, 'U' ) 147: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 148: INFO = -1 149: ELSE IF( N.LT.0 ) THEN 150: INFO = -2 151: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 152: INFO = -4 153: END IF 154: IF( INFO.NE.0 ) THEN 155: CALL XERBLA( 'ZHETD2', -INFO ) 156: RETURN 157: END IF 158: * 159: * Quick return if possible 160: * 161: IF( N.LE.0 ) 162: $ RETURN 163: * 164: IF( UPPER ) THEN 165: * 166: * Reduce the upper triangle of A 167: * 168: A( N, N ) = DBLE( A( N, N ) ) 169: DO 10 I = N - 1, 1, -1 170: * 171: * Generate elementary reflector H(i) = I - tau * v * v' 172: * to annihilate A(1:i-1,i+1) 173: * 174: ALPHA = A( I, I+1 ) 175: CALL ZLARFG( I, ALPHA, A( 1, I+1 ), 1, TAUI ) 176: E( I ) = ALPHA 177: * 178: IF( TAUI.NE.ZERO ) THEN 179: * 180: * Apply H(i) from both sides to A(1:i,1:i) 181: * 182: A( I, I+1 ) = ONE 183: * 184: * Compute x := tau * A * v storing x in TAU(1:i) 185: * 186: CALL ZHEMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, 187: $ TAU, 1 ) 188: * 189: * Compute w := x - 1/2 * tau * (x'*v) * v 190: * 191: ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, A( 1, I+1 ), 1 ) 192: CALL ZAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 ) 193: * 194: * Apply the transformation as a rank-2 update: 195: * A := A - v * w' - w * v' 196: * 197: CALL ZHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, 198: $ LDA ) 199: * 200: ELSE 201: A( I, I ) = DBLE( A( I, I ) ) 202: END IF 203: A( I, I+1 ) = E( I ) 204: D( I+1 ) = A( I+1, I+1 ) 205: TAU( I ) = TAUI 206: 10 CONTINUE 207: D( 1 ) = A( 1, 1 ) 208: ELSE 209: * 210: * Reduce the lower triangle of A 211: * 212: A( 1, 1 ) = DBLE( A( 1, 1 ) ) 213: DO 20 I = 1, N - 1 214: * 215: * Generate elementary reflector H(i) = I - tau * v * v' 216: * to annihilate A(i+2:n,i) 217: * 218: ALPHA = A( I+1, I ) 219: CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAUI ) 220: E( I ) = ALPHA 221: * 222: IF( TAUI.NE.ZERO ) THEN 223: * 224: * Apply H(i) from both sides to A(i+1:n,i+1:n) 225: * 226: A( I+1, I ) = ONE 227: * 228: * Compute x := tau * A * v storing y in TAU(i:n-1) 229: * 230: CALL ZHEMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, 231: $ A( I+1, I ), 1, ZERO, TAU( I ), 1 ) 232: * 233: * Compute w := x - 1/2 * tau * (x'*v) * v 234: * 235: ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, A( I+1, I ), 236: $ 1 ) 237: CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 ) 238: * 239: * Apply the transformation as a rank-2 update: 240: * A := A - v * w' - w * v' 241: * 242: CALL ZHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, 243: $ A( I+1, I+1 ), LDA ) 244: * 245: ELSE 246: A( I+1, I+1 ) = DBLE( A( I+1, I+1 ) ) 247: END IF 248: A( I+1, I ) = E( I ) 249: D( I ) = A( I, I ) 250: TAU( I ) = TAUI 251: 20 CONTINUE 252: D( N ) = A( N, N ) 253: END IF 254: * 255: RETURN 256: * 257: * End of ZHETD2 258: * 259: END