![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHPTRD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHPTRD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptrd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptrd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptrd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, N 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION D( * ), E( * ) 00029 * COMPLEX*16 AP( * ), TAU( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> ZHPTRD reduces a complex Hermitian matrix A stored in packed form to 00039 *> real symmetric tridiagonal form T by a unitary similarity 00040 *> transformation: Q**H * A * Q = T. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in] UPLO 00047 *> \verbatim 00048 *> UPLO is CHARACTER*1 00049 *> = 'U': Upper triangle of A is stored; 00050 *> = 'L': Lower triangle of A is stored. 00051 *> \endverbatim 00052 *> 00053 *> \param[in] N 00054 *> \verbatim 00055 *> N is INTEGER 00056 *> The order of the matrix A. N >= 0. 00057 *> \endverbatim 00058 *> 00059 *> \param[in,out] AP 00060 *> \verbatim 00061 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 00062 *> On entry, the upper or lower triangle of the Hermitian matrix 00063 *> A, packed columnwise in a linear array. The j-th column of A 00064 *> is stored in the array AP as follows: 00065 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00066 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 00067 *> On exit, if UPLO = 'U', the diagonal and first superdiagonal 00068 *> of A are overwritten by the corresponding elements of the 00069 *> tridiagonal matrix T, and the elements above the first 00070 *> superdiagonal, with the array TAU, represent the unitary 00071 *> matrix Q as a product of elementary reflectors; if UPLO 00072 *> = 'L', the diagonal and first subdiagonal of A are over- 00073 *> written by the corresponding elements of the tridiagonal 00074 *> matrix T, and the elements below the first subdiagonal, with 00075 *> the array TAU, represent the unitary matrix Q as a product 00076 *> of elementary reflectors. See Further Details. 00077 *> \endverbatim 00078 *> 00079 *> \param[out] D 00080 *> \verbatim 00081 *> D is DOUBLE PRECISION array, dimension (N) 00082 *> The diagonal elements of the tridiagonal matrix T: 00083 *> D(i) = A(i,i). 00084 *> \endverbatim 00085 *> 00086 *> \param[out] E 00087 *> \verbatim 00088 *> E is DOUBLE PRECISION array, dimension (N-1) 00089 *> The off-diagonal elements of the tridiagonal matrix T: 00090 *> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 00091 *> \endverbatim 00092 *> 00093 *> \param[out] TAU 00094 *> \verbatim 00095 *> TAU is COMPLEX*16 array, dimension (N-1) 00096 *> The scalar factors of the elementary reflectors (see Further 00097 *> Details). 00098 *> \endverbatim 00099 *> 00100 *> \param[out] INFO 00101 *> \verbatim 00102 *> INFO is INTEGER 00103 *> = 0: successful exit 00104 *> < 0: if INFO = -i, the i-th argument had an illegal value 00105 *> \endverbatim 00106 * 00107 * Authors: 00108 * ======== 00109 * 00110 *> \author Univ. of Tennessee 00111 *> \author Univ. of California Berkeley 00112 *> \author Univ. of Colorado Denver 00113 *> \author NAG Ltd. 00114 * 00115 *> \date November 2011 00116 * 00117 *> \ingroup complex16OTHERcomputational 00118 * 00119 *> \par Further Details: 00120 * ===================== 00121 *> 00122 *> \verbatim 00123 *> 00124 *> If UPLO = 'U', the matrix Q is represented as a product of elementary 00125 *> reflectors 00126 *> 00127 *> Q = H(n-1) . . . H(2) H(1). 00128 *> 00129 *> Each H(i) has the form 00130 *> 00131 *> H(i) = I - tau * v * v**H 00132 *> 00133 *> where tau is a complex scalar, and v is a complex vector with 00134 *> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP, 00135 *> overwriting A(1:i-1,i+1), and tau is stored in TAU(i). 00136 *> 00137 *> If UPLO = 'L', the matrix Q is represented as a product of elementary 00138 *> reflectors 00139 *> 00140 *> Q = H(1) H(2) . . . H(n-1). 00141 *> 00142 *> Each H(i) has the form 00143 *> 00144 *> H(i) = I - tau * v * v**H 00145 *> 00146 *> where tau is a complex scalar, and v is a complex vector with 00147 *> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP, 00148 *> overwriting A(i+2:n,i), and tau is stored in TAU(i). 00149 *> \endverbatim 00150 *> 00151 * ===================================================================== 00152 SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO ) 00153 * 00154 * -- LAPACK computational routine (version 3.4.0) -- 00155 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00157 * November 2011 00158 * 00159 * .. Scalar Arguments .. 00160 CHARACTER UPLO 00161 INTEGER INFO, N 00162 * .. 00163 * .. Array Arguments .. 00164 DOUBLE PRECISION D( * ), E( * ) 00165 COMPLEX*16 AP( * ), TAU( * ) 00166 * .. 00167 * 00168 * ===================================================================== 00169 * 00170 * .. Parameters .. 00171 COMPLEX*16 ONE, ZERO, HALF 00172 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 00173 $ ZERO = ( 0.0D+0, 0.0D+0 ), 00174 $ HALF = ( 0.5D+0, 0.0D+0 ) ) 00175 * .. 00176 * .. Local Scalars .. 00177 LOGICAL UPPER 00178 INTEGER I, I1, I1I1, II 00179 COMPLEX*16 ALPHA, TAUI 00180 * .. 00181 * .. External Subroutines .. 00182 EXTERNAL XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG 00183 * .. 00184 * .. External Functions .. 00185 LOGICAL LSAME 00186 COMPLEX*16 ZDOTC 00187 EXTERNAL LSAME, ZDOTC 00188 * .. 00189 * .. Intrinsic Functions .. 00190 INTRINSIC DBLE 00191 * .. 00192 * .. Executable Statements .. 00193 * 00194 * Test the input parameters 00195 * 00196 INFO = 0 00197 UPPER = LSAME( UPLO, 'U' ) 00198 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00199 INFO = -1 00200 ELSE IF( N.LT.0 ) THEN 00201 INFO = -2 00202 END IF 00203 IF( INFO.NE.0 ) THEN 00204 CALL XERBLA( 'ZHPTRD', -INFO ) 00205 RETURN 00206 END IF 00207 * 00208 * Quick return if possible 00209 * 00210 IF( N.LE.0 ) 00211 $ RETURN 00212 * 00213 IF( UPPER ) THEN 00214 * 00215 * Reduce the upper triangle of A. 00216 * I1 is the index in AP of A(1,I+1). 00217 * 00218 I1 = N*( N-1 ) / 2 + 1 00219 AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) ) 00220 DO 10 I = N - 1, 1, -1 00221 * 00222 * Generate elementary reflector H(i) = I - tau * v * v**H 00223 * to annihilate A(1:i-1,i+1) 00224 * 00225 ALPHA = AP( I1+I-1 ) 00226 CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI ) 00227 E( I ) = ALPHA 00228 * 00229 IF( TAUI.NE.ZERO ) THEN 00230 * 00231 * Apply H(i) from both sides to A(1:i,1:i) 00232 * 00233 AP( I1+I-1 ) = ONE 00234 * 00235 * Compute y := tau * A * v storing y in TAU(1:i) 00236 * 00237 CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU, 00238 $ 1 ) 00239 * 00240 * Compute w := y - 1/2 * tau * (y**H *v) * v 00241 * 00242 ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 ) 00243 CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 ) 00244 * 00245 * Apply the transformation as a rank-2 update: 00246 * A := A - v * w**H - w * v**H 00247 * 00248 CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP ) 00249 * 00250 END IF 00251 AP( I1+I-1 ) = E( I ) 00252 D( I+1 ) = AP( I1+I ) 00253 TAU( I ) = TAUI 00254 I1 = I1 - I 00255 10 CONTINUE 00256 D( 1 ) = AP( 1 ) 00257 ELSE 00258 * 00259 * Reduce the lower triangle of A. II is the index in AP of 00260 * A(i,i) and I1I1 is the index of A(i+1,i+1). 00261 * 00262 II = 1 00263 AP( 1 ) = DBLE( AP( 1 ) ) 00264 DO 20 I = 1, N - 1 00265 I1I1 = II + N - I + 1 00266 * 00267 * Generate elementary reflector H(i) = I - tau * v * v**H 00268 * to annihilate A(i+2:n,i) 00269 * 00270 ALPHA = AP( II+1 ) 00271 CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI ) 00272 E( I ) = ALPHA 00273 * 00274 IF( TAUI.NE.ZERO ) THEN 00275 * 00276 * Apply H(i) from both sides to A(i+1:n,i+1:n) 00277 * 00278 AP( II+1 ) = ONE 00279 * 00280 * Compute y := tau * A * v storing y in TAU(i:n-1) 00281 * 00282 CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1, 00283 $ ZERO, TAU( I ), 1 ) 00284 * 00285 * Compute w := y - 1/2 * tau * (y**H *v) * v 00286 * 00287 ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ), 00288 $ 1 ) 00289 CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 ) 00290 * 00291 * Apply the transformation as a rank-2 update: 00292 * A := A - v * w**H - w * v**H 00293 * 00294 CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1, 00295 $ AP( I1I1 ) ) 00296 * 00297 END IF 00298 AP( II+1 ) = E( I ) 00299 D( I ) = AP( II ) 00300 TAU( I ) = TAUI 00301 II = I1I1 00302 20 CONTINUE 00303 D( N ) = AP( II ) 00304 END IF 00305 * 00306 RETURN 00307 * 00308 * End of ZHPTRD 00309 * 00310 END