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