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