1*> \brief \b ZHPTRD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZHPTRD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptrd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptrd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptrd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, N 26* .. 27* .. Array Arguments .. 28* DOUBLE PRECISION D( * ), E( * ) 29* COMPLEX*16 AP( * ), TAU( * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> ZHPTRD reduces a complex Hermitian matrix A stored in packed form to 39*> real symmetric tridiagonal form T by a unitary similarity 40*> transformation: Q**H * A * Q = T. 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] UPLO 47*> \verbatim 48*> UPLO is CHARACTER*1 49*> = 'U': Upper triangle of A is stored; 50*> = 'L': Lower triangle of A is stored. 51*> \endverbatim 52*> 53*> \param[in] N 54*> \verbatim 55*> N is INTEGER 56*> The order of the matrix A. N >= 0. 57*> \endverbatim 58*> 59*> \param[in,out] AP 60*> \verbatim 61*> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 62*> On entry, the upper or lower triangle of the Hermitian matrix 63*> A, packed columnwise in a linear array. The j-th column of A 64*> is stored in the array AP as follows: 65*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 66*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 67*> On exit, if UPLO = 'U', the diagonal and first superdiagonal 68*> of A are overwritten by the corresponding elements of the 69*> tridiagonal matrix T, and the elements above the first 70*> superdiagonal, with the array TAU, represent the unitary 71*> matrix Q as a product of elementary reflectors; if UPLO 72*> = 'L', the diagonal and first subdiagonal of A are over- 73*> written by the corresponding elements of the tridiagonal 74*> matrix T, and the elements below the first subdiagonal, with 75*> the array TAU, represent the unitary matrix Q as a product 76*> of elementary reflectors. See Further Details. 77*> \endverbatim 78*> 79*> \param[out] D 80*> \verbatim 81*> D is DOUBLE PRECISION array, dimension (N) 82*> The diagonal elements of the tridiagonal matrix T: 83*> D(i) = A(i,i). 84*> \endverbatim 85*> 86*> \param[out] E 87*> \verbatim 88*> E is DOUBLE PRECISION array, dimension (N-1) 89*> The off-diagonal elements of the tridiagonal matrix T: 90*> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 91*> \endverbatim 92*> 93*> \param[out] TAU 94*> \verbatim 95*> TAU is COMPLEX*16 array, dimension (N-1) 96*> The scalar factors of the elementary reflectors (see Further 97*> Details). 98*> \endverbatim 99*> 100*> \param[out] INFO 101*> \verbatim 102*> INFO is INTEGER 103*> = 0: successful exit 104*> < 0: if INFO = -i, the i-th argument had an illegal value 105*> \endverbatim 106* 107* Authors: 108* ======== 109* 110*> \author Univ. of Tennessee 111*> \author Univ. of California Berkeley 112*> \author Univ. of Colorado Denver 113*> \author NAG Ltd. 114* 115*> \date November 2011 116* 117*> \ingroup complex16OTHERcomputational 118* 119*> \par Further Details: 120* ===================== 121*> 122*> \verbatim 123*> 124*> If UPLO = 'U', the matrix Q is represented as a product of elementary 125*> reflectors 126*> 127*> Q = H(n-1) . . . H(2) H(1). 128*> 129*> Each H(i) has the form 130*> 131*> H(i) = I - tau * v * v**H 132*> 133*> where tau is a complex scalar, and v is a complex vector with 134*> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP, 135*> overwriting A(1:i-1,i+1), and tau is stored in TAU(i). 136*> 137*> If UPLO = 'L', the matrix Q is represented as a product of elementary 138*> reflectors 139*> 140*> Q = H(1) H(2) . . . H(n-1). 141*> 142*> Each H(i) has the form 143*> 144*> H(i) = I - tau * v * v**H 145*> 146*> where tau is a complex scalar, and v is a complex vector with 147*> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP, 148*> overwriting A(i+2:n,i), and tau is stored in TAU(i). 149*> \endverbatim 150*> 151* ===================================================================== 152 SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO ) 153* 154* -- LAPACK computational routine (version 3.4.0) -- 155* -- LAPACK is a software package provided by Univ. of Tennessee, -- 156* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 157* November 2011 158* 159* .. Scalar Arguments .. 160 CHARACTER UPLO 161 INTEGER INFO, N 162* .. 163* .. Array Arguments .. 164 DOUBLE PRECISION D( * ), E( * ) 165 COMPLEX*16 AP( * ), TAU( * ) 166* .. 167* 168* ===================================================================== 169* 170* .. Parameters .. 171 COMPLEX*16 ONE, ZERO, HALF 172 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 173 $ ZERO = ( 0.0D+0, 0.0D+0 ), 174 $ HALF = ( 0.5D+0, 0.0D+0 ) ) 175* .. 176* .. Local Scalars .. 177 LOGICAL UPPER 178 INTEGER I, I1, I1I1, II 179 COMPLEX*16 ALPHA, TAUI 180* .. 181* .. External Subroutines .. 182 EXTERNAL XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG 183* .. 184* .. External Functions .. 185 LOGICAL LSAME 186 COMPLEX*16 ZDOTC 187 EXTERNAL LSAME, ZDOTC 188* .. 189* .. Intrinsic Functions .. 190 INTRINSIC DBLE 191* .. 192* .. Executable Statements .. 193* 194* Test the input parameters 195* 196 INFO = 0 197 UPPER = LSAME( UPLO, 'U' ) 198 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 199 INFO = -1 200 ELSE IF( N.LT.0 ) THEN 201 INFO = -2 202 END IF 203 IF( INFO.NE.0 ) THEN 204 CALL XERBLA( 'ZHPTRD', -INFO ) 205 RETURN 206 END IF 207* 208* Quick return if possible 209* 210 IF( N.LE.0 ) 211 $ RETURN 212* 213 IF( UPPER ) THEN 214* 215* Reduce the upper triangle of A. 216* I1 is the index in AP of A(1,I+1). 217* 218 I1 = N*( N-1 ) / 2 + 1 219 AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) ) 220 DO 10 I = N - 1, 1, -1 221* 222* Generate elementary reflector H(i) = I - tau * v * v**H 223* to annihilate A(1:i-1,i+1) 224* 225 ALPHA = AP( I1+I-1 ) 226 CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI ) 227 E( I ) = ALPHA 228* 229 IF( TAUI.NE.ZERO ) THEN 230* 231* Apply H(i) from both sides to A(1:i,1:i) 232* 233 AP( I1+I-1 ) = ONE 234* 235* Compute y := tau * A * v storing y in TAU(1:i) 236* 237 CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU, 238 $ 1 ) 239* 240* Compute w := y - 1/2 * tau * (y**H *v) * v 241* 242 ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 ) 243 CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 ) 244* 245* Apply the transformation as a rank-2 update: 246* A := A - v * w**H - w * v**H 247* 248 CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP ) 249* 250 END IF 251 AP( I1+I-1 ) = E( I ) 252 D( I+1 ) = AP( I1+I ) 253 TAU( I ) = TAUI 254 I1 = I1 - I 255 10 CONTINUE 256 D( 1 ) = AP( 1 ) 257 ELSE 258* 259* Reduce the lower triangle of A. II is the index in AP of 260* A(i,i) and I1I1 is the index of A(i+1,i+1). 261* 262 II = 1 263 AP( 1 ) = DBLE( AP( 1 ) ) 264 DO 20 I = 1, N - 1 265 I1I1 = II + N - I + 1 266* 267* Generate elementary reflector H(i) = I - tau * v * v**H 268* to annihilate A(i+2:n,i) 269* 270 ALPHA = AP( II+1 ) 271 CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI ) 272 E( I ) = ALPHA 273* 274 IF( TAUI.NE.ZERO ) THEN 275* 276* Apply H(i) from both sides to A(i+1:n,i+1:n) 277* 278 AP( II+1 ) = ONE 279* 280* Compute y := tau * A * v storing y in TAU(i:n-1) 281* 282 CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1, 283 $ ZERO, TAU( I ), 1 ) 284* 285* Compute w := y - 1/2 * tau * (y**H *v) * v 286* 287 ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ), 288 $ 1 ) 289 CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 ) 290* 291* Apply the transformation as a rank-2 update: 292* A := A - v * w**H - w * v**H 293* 294 CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1, 295 $ AP( I1I1 ) ) 296* 297 END IF 298 AP( II+1 ) = E( I ) 299 D( I ) = AP( II ) 300 TAU( I ) = TAUI 301 II = I1I1 302 20 CONTINUE 303 D( N ) = AP( II ) 304 END IF 305* 306 RETURN 307* 308* End of ZHPTRD 309* 310 END 311