00001 SUBROUTINE CHETD2( UPLO, N, A, LDA, D, E, TAU, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INFO, LDA, N 00011 * .. 00012 * .. Array Arguments .. 00013 REAL D( * ), E( * ) 00014 COMPLEX A( LDA, * ), TAU( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CHETD2 reduces a complex Hermitian matrix A to real symmetric 00021 * tridiagonal form T by a unitary similarity transformation: 00022 * Q' * A * Q = T. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * UPLO (input) CHARACTER*1 00028 * Specifies whether the upper or lower triangular part of the 00029 * Hermitian matrix A is stored: 00030 * = 'U': Upper triangular 00031 * = 'L': Lower triangular 00032 * 00033 * N (input) INTEGER 00034 * The order of the matrix A. N >= 0. 00035 * 00036 * A (input/output) COMPLEX array, dimension (LDA,N) 00037 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading 00038 * n-by-n upper triangular part of A contains the upper 00039 * triangular part of the matrix A, and the strictly lower 00040 * triangular part of A is not referenced. If UPLO = 'L', the 00041 * leading n-by-n lower triangular part of A contains the lower 00042 * triangular part of the matrix A, and the strictly upper 00043 * triangular part of A is not referenced. 00044 * On exit, if UPLO = 'U', the diagonal and first superdiagonal 00045 * of A are overwritten by the corresponding elements of the 00046 * tridiagonal matrix T, and the elements above the first 00047 * superdiagonal, with the array TAU, represent the unitary 00048 * matrix Q as a product of elementary reflectors; if UPLO 00049 * = 'L', the diagonal and first subdiagonal of A are over- 00050 * written by the corresponding elements of the tridiagonal 00051 * matrix T, and the elements below the first subdiagonal, with 00052 * the array TAU, represent the unitary matrix Q as a product 00053 * of elementary reflectors. See Further Details. 00054 * 00055 * LDA (input) INTEGER 00056 * The leading dimension of the array A. LDA >= max(1,N). 00057 * 00058 * D (output) REAL array, dimension (N) 00059 * The diagonal elements of the tridiagonal matrix T: 00060 * D(i) = A(i,i). 00061 * 00062 * E (output) REAL array, dimension (N-1) 00063 * The off-diagonal elements of the tridiagonal matrix T: 00064 * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 00065 * 00066 * TAU (output) COMPLEX array, dimension (N-1) 00067 * The scalar factors of the elementary reflectors (see Further 00068 * Details). 00069 * 00070 * INFO (output) INTEGER 00071 * = 0: successful exit 00072 * < 0: if INFO = -i, the i-th argument had an illegal value. 00073 * 00074 * Further Details 00075 * =============== 00076 * 00077 * If UPLO = 'U', the matrix Q is represented as a product of elementary 00078 * reflectors 00079 * 00080 * Q = H(n-1) . . . H(2) H(1). 00081 * 00082 * Each H(i) has the form 00083 * 00084 * H(i) = I - tau * v * v' 00085 * 00086 * where tau is a complex scalar, and v is a complex vector with 00087 * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 00088 * A(1:i-1,i+1), and tau in TAU(i). 00089 * 00090 * If UPLO = 'L', the matrix Q is represented as a product of elementary 00091 * reflectors 00092 * 00093 * Q = H(1) H(2) . . . H(n-1). 00094 * 00095 * Each H(i) has the form 00096 * 00097 * H(i) = I - tau * v * v' 00098 * 00099 * where tau is a complex scalar, and v is a complex vector with 00100 * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 00101 * and tau in TAU(i). 00102 * 00103 * The contents of A on exit are illustrated by the following examples 00104 * with n = 5: 00105 * 00106 * if UPLO = 'U': if UPLO = 'L': 00107 * 00108 * ( d e v2 v3 v4 ) ( d ) 00109 * ( d e v3 v4 ) ( e d ) 00110 * ( d e v4 ) ( v1 e d ) 00111 * ( d e ) ( v1 v2 e d ) 00112 * ( d ) ( v1 v2 v3 e d ) 00113 * 00114 * where d and e denote diagonal and off-diagonal elements of T, and vi 00115 * denotes an element of the vector defining H(i). 00116 * 00117 * ===================================================================== 00118 * 00119 * .. Parameters .. 00120 COMPLEX ONE, ZERO, HALF 00121 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), 00122 $ ZERO = ( 0.0E+0, 0.0E+0 ), 00123 $ HALF = ( 0.5E+0, 0.0E+0 ) ) 00124 * .. 00125 * .. Local Scalars .. 00126 LOGICAL UPPER 00127 INTEGER I 00128 COMPLEX ALPHA, TAUI 00129 * .. 00130 * .. External Subroutines .. 00131 EXTERNAL CAXPY, CHEMV, CHER2, CLARFG, XERBLA 00132 * .. 00133 * .. External Functions .. 00134 LOGICAL LSAME 00135 COMPLEX CDOTC 00136 EXTERNAL LSAME, CDOTC 00137 * .. 00138 * .. Intrinsic Functions .. 00139 INTRINSIC MAX, MIN, REAL 00140 * .. 00141 * .. Executable Statements .. 00142 * 00143 * Test the input parameters 00144 * 00145 INFO = 0 00146 UPPER = LSAME( UPLO, 'U' ) 00147 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00148 INFO = -1 00149 ELSE IF( N.LT.0 ) THEN 00150 INFO = -2 00151 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00152 INFO = -4 00153 END IF 00154 IF( INFO.NE.0 ) THEN 00155 CALL XERBLA( 'CHETD2', -INFO ) 00156 RETURN 00157 END IF 00158 * 00159 * Quick return if possible 00160 * 00161 IF( N.LE.0 ) 00162 $ RETURN 00163 * 00164 IF( UPPER ) THEN 00165 * 00166 * Reduce the upper triangle of A 00167 * 00168 A( N, N ) = REAL( A( N, N ) ) 00169 DO 10 I = N - 1, 1, -1 00170 * 00171 * Generate elementary reflector H(i) = I - tau * v * v' 00172 * to annihilate A(1:i-1,i+1) 00173 * 00174 ALPHA = A( I, I+1 ) 00175 CALL CLARFG( I, ALPHA, A( 1, I+1 ), 1, TAUI ) 00176 E( I ) = ALPHA 00177 * 00178 IF( TAUI.NE.ZERO ) THEN 00179 * 00180 * Apply H(i) from both sides to A(1:i,1:i) 00181 * 00182 A( I, I+1 ) = ONE 00183 * 00184 * Compute x := tau * A * v storing x in TAU(1:i) 00185 * 00186 CALL CHEMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, 00187 $ TAU, 1 ) 00188 * 00189 * Compute w := x - 1/2 * tau * (x'*v) * v 00190 * 00191 ALPHA = -HALF*TAUI*CDOTC( I, TAU, 1, A( 1, I+1 ), 1 ) 00192 CALL CAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 ) 00193 * 00194 * Apply the transformation as a rank-2 update: 00195 * A := A - v * w' - w * v' 00196 * 00197 CALL CHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, 00198 $ LDA ) 00199 * 00200 ELSE 00201 A( I, I ) = REAL( A( I, I ) ) 00202 END IF 00203 A( I, I+1 ) = E( I ) 00204 D( I+1 ) = A( I+1, I+1 ) 00205 TAU( I ) = TAUI 00206 10 CONTINUE 00207 D( 1 ) = A( 1, 1 ) 00208 ELSE 00209 * 00210 * Reduce the lower triangle of A 00211 * 00212 A( 1, 1 ) = REAL( A( 1, 1 ) ) 00213 DO 20 I = 1, N - 1 00214 * 00215 * Generate elementary reflector H(i) = I - tau * v * v' 00216 * to annihilate A(i+2:n,i) 00217 * 00218 ALPHA = A( I+1, I ) 00219 CALL CLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAUI ) 00220 E( I ) = ALPHA 00221 * 00222 IF( TAUI.NE.ZERO ) THEN 00223 * 00224 * Apply H(i) from both sides to A(i+1:n,i+1:n) 00225 * 00226 A( I+1, I ) = ONE 00227 * 00228 * Compute x := tau * A * v storing y in TAU(i:n-1) 00229 * 00230 CALL CHEMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, 00231 $ A( I+1, I ), 1, ZERO, TAU( I ), 1 ) 00232 * 00233 * Compute w := x - 1/2 * tau * (x'*v) * v 00234 * 00235 ALPHA = -HALF*TAUI*CDOTC( N-I, TAU( I ), 1, A( I+1, I ), 00236 $ 1 ) 00237 CALL CAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 ) 00238 * 00239 * Apply the transformation as a rank-2 update: 00240 * A := A - v * w' - w * v' 00241 * 00242 CALL CHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, 00243 $ A( I+1, I+1 ), LDA ) 00244 * 00245 ELSE 00246 A( I+1, I+1 ) = REAL( A( I+1, I+1 ) ) 00247 END IF 00248 A( I+1, I ) = E( I ) 00249 D( I ) = A( I, I ) 00250 TAU( I ) = TAUI 00251 20 CONTINUE 00252 D( N ) = A( N, N ) 00253 END IF 00254 * 00255 RETURN 00256 * 00257 * End of CHETD2 00258 * 00259 END
1.7.2