00001 SUBROUTINE CHPTRD( UPLO, N, AP, 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, N 00011 * .. 00012 * .. Array Arguments .. 00013 REAL D( * ), E( * ) 00014 COMPLEX AP( * ), TAU( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CHPTRD reduces a complex Hermitian matrix A stored in packed form to 00021 * real symmetric tridiagonal form T by a unitary similarity 00022 * transformation: Q**H * A * Q = T. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * UPLO (input) CHARACTER*1 00028 * = 'U': Upper triangle of A is stored; 00029 * = 'L': Lower triangle of A is stored. 00030 * 00031 * N (input) INTEGER 00032 * The order of the matrix A. N >= 0. 00033 * 00034 * AP (input/output) COMPLEX array, dimension (N*(N+1)/2) 00035 * On entry, the upper or lower triangle of the Hermitian matrix 00036 * A, packed columnwise in a linear array. The j-th column of A 00037 * is stored in the array AP as follows: 00038 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00039 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 00040 * On exit, if UPLO = 'U', the diagonal and first superdiagonal 00041 * of A are overwritten by the corresponding elements of the 00042 * tridiagonal matrix T, and the elements above the first 00043 * superdiagonal, with the array TAU, represent the unitary 00044 * matrix Q as a product of elementary reflectors; if UPLO 00045 * = 'L', the diagonal and first subdiagonal of A are over- 00046 * written by the corresponding elements of the tridiagonal 00047 * matrix T, and the elements below the first subdiagonal, with 00048 * the array TAU, represent the unitary matrix Q as a product 00049 * of elementary reflectors. See Further Details. 00050 * 00051 * D (output) REAL array, dimension (N) 00052 * The diagonal elements of the tridiagonal matrix T: 00053 * D(i) = A(i,i). 00054 * 00055 * E (output) REAL array, dimension (N-1) 00056 * The off-diagonal elements of the tridiagonal matrix T: 00057 * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 00058 * 00059 * TAU (output) COMPLEX array, dimension (N-1) 00060 * The scalar factors of the elementary reflectors (see Further 00061 * Details). 00062 * 00063 * INFO (output) INTEGER 00064 * = 0: successful exit 00065 * < 0: if INFO = -i, the i-th argument had an illegal value 00066 * 00067 * Further Details 00068 * =============== 00069 * 00070 * If UPLO = 'U', the matrix Q is represented as a product of elementary 00071 * reflectors 00072 * 00073 * Q = H(n-1) . . . H(2) H(1). 00074 * 00075 * Each H(i) has the form 00076 * 00077 * H(i) = I - tau * v * v' 00078 * 00079 * where tau is a complex scalar, and v is a complex vector with 00080 * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP, 00081 * overwriting A(1:i-1,i+1), and tau is stored in TAU(i). 00082 * 00083 * If UPLO = 'L', the matrix Q is represented as a product of elementary 00084 * reflectors 00085 * 00086 * Q = H(1) H(2) . . . H(n-1). 00087 * 00088 * Each H(i) has the form 00089 * 00090 * H(i) = I - tau * v * v' 00091 * 00092 * where tau is a complex scalar, and v is a complex vector with 00093 * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP, 00094 * overwriting A(i+2:n,i), and tau is stored in TAU(i). 00095 * 00096 * ===================================================================== 00097 * 00098 * .. Parameters .. 00099 COMPLEX ONE, ZERO, HALF 00100 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), 00101 $ ZERO = ( 0.0E+0, 0.0E+0 ), 00102 $ HALF = ( 0.5E+0, 0.0E+0 ) ) 00103 * .. 00104 * .. Local Scalars .. 00105 LOGICAL UPPER 00106 INTEGER I, I1, I1I1, II 00107 COMPLEX ALPHA, TAUI 00108 * .. 00109 * .. External Subroutines .. 00110 EXTERNAL CAXPY, CHPMV, CHPR2, CLARFG, XERBLA 00111 * .. 00112 * .. External Functions .. 00113 LOGICAL LSAME 00114 COMPLEX CDOTC 00115 EXTERNAL LSAME, CDOTC 00116 * .. 00117 * .. Intrinsic Functions .. 00118 INTRINSIC REAL 00119 * .. 00120 * .. Executable Statements .. 00121 * 00122 * Test the input parameters 00123 * 00124 INFO = 0 00125 UPPER = LSAME( UPLO, 'U' ) 00126 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00127 INFO = -1 00128 ELSE IF( N.LT.0 ) THEN 00129 INFO = -2 00130 END IF 00131 IF( INFO.NE.0 ) THEN 00132 CALL XERBLA( 'CHPTRD', -INFO ) 00133 RETURN 00134 END IF 00135 * 00136 * Quick return if possible 00137 * 00138 IF( N.LE.0 ) 00139 $ RETURN 00140 * 00141 IF( UPPER ) THEN 00142 * 00143 * Reduce the upper triangle of A. 00144 * I1 is the index in AP of A(1,I+1). 00145 * 00146 I1 = N*( N-1 ) / 2 + 1 00147 AP( I1+N-1 ) = REAL( AP( I1+N-1 ) ) 00148 DO 10 I = N - 1, 1, -1 00149 * 00150 * Generate elementary reflector H(i) = I - tau * v * v' 00151 * to annihilate A(1:i-1,i+1) 00152 * 00153 ALPHA = AP( I1+I-1 ) 00154 CALL CLARFG( I, ALPHA, AP( I1 ), 1, TAUI ) 00155 E( I ) = ALPHA 00156 * 00157 IF( TAUI.NE.ZERO ) THEN 00158 * 00159 * Apply H(i) from both sides to A(1:i,1:i) 00160 * 00161 AP( I1+I-1 ) = ONE 00162 * 00163 * Compute y := tau * A * v storing y in TAU(1:i) 00164 * 00165 CALL CHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU, 00166 $ 1 ) 00167 * 00168 * Compute w := y - 1/2 * tau * (y'*v) * v 00169 * 00170 ALPHA = -HALF*TAUI*CDOTC( I, TAU, 1, AP( I1 ), 1 ) 00171 CALL CAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 ) 00172 * 00173 * Apply the transformation as a rank-2 update: 00174 * A := A - v * w' - w * v' 00175 * 00176 CALL CHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP ) 00177 * 00178 END IF 00179 AP( I1+I-1 ) = E( I ) 00180 D( I+1 ) = AP( I1+I ) 00181 TAU( I ) = TAUI 00182 I1 = I1 - I 00183 10 CONTINUE 00184 D( 1 ) = AP( 1 ) 00185 ELSE 00186 * 00187 * Reduce the lower triangle of A. II is the index in AP of 00188 * A(i,i) and I1I1 is the index of A(i+1,i+1). 00189 * 00190 II = 1 00191 AP( 1 ) = REAL( AP( 1 ) ) 00192 DO 20 I = 1, N - 1 00193 I1I1 = II + N - I + 1 00194 * 00195 * Generate elementary reflector H(i) = I - tau * v * v' 00196 * to annihilate A(i+2:n,i) 00197 * 00198 ALPHA = AP( II+1 ) 00199 CALL CLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI ) 00200 E( I ) = ALPHA 00201 * 00202 IF( TAUI.NE.ZERO ) THEN 00203 * 00204 * Apply H(i) from both sides to A(i+1:n,i+1:n) 00205 * 00206 AP( II+1 ) = ONE 00207 * 00208 * Compute y := tau * A * v storing y in TAU(i:n-1) 00209 * 00210 CALL CHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1, 00211 $ ZERO, TAU( I ), 1 ) 00212 * 00213 * Compute w := y - 1/2 * tau * (y'*v) * v 00214 * 00215 ALPHA = -HALF*TAUI*CDOTC( N-I, TAU( I ), 1, AP( II+1 ), 00216 $ 1 ) 00217 CALL CAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 ) 00218 * 00219 * Apply the transformation as a rank-2 update: 00220 * A := A - v * w' - w * v' 00221 * 00222 CALL CHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1, 00223 $ AP( I1I1 ) ) 00224 * 00225 END IF 00226 AP( II+1 ) = E( I ) 00227 D( I ) = AP( II ) 00228 TAU( I ) = TAUI 00229 II = I1I1 00230 20 CONTINUE 00231 D( N ) = AP( II ) 00232 END IF 00233 * 00234 RETURN 00235 * 00236 * End of CHPTRD 00237 * 00238 END
1.7.2