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*> \ingroup complex16OTHERcomputational 116* 117*> \par Further Details: 118* ===================== 119*> 120*> \verbatim 121*> 122*> If UPLO = 'U', the matrix Q is represented as a product of elementary 123*> reflectors 124*> 125*> Q = H(n-1) . . . H(2) H(1). 126*> 127*> Each H(i) has the form 128*> 129*> H(i) = I - tau * v * v**H 130*> 131*> where tau is a complex scalar, and v is a complex vector with 132*> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP, 133*> overwriting A(1:i-1,i+1), and tau is stored in TAU(i). 134*> 135*> If UPLO = 'L', the matrix Q is represented as a product of elementary 136*> reflectors 137*> 138*> Q = H(1) H(2) . . . H(n-1). 139*> 140*> Each H(i) has the form 141*> 142*> H(i) = I - tau * v * v**H 143*> 144*> where tau is a complex scalar, and v is a complex vector with 145*> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP, 146*> overwriting A(i+2:n,i), and tau is stored in TAU(i). 147*> \endverbatim 148*> 149* ===================================================================== 150 SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO ) 151* 152* -- LAPACK computational routine -- 153* -- LAPACK is a software package provided by Univ. of Tennessee, -- 154* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 155* 156* .. Scalar Arguments .. 157 CHARACTER UPLO 158 INTEGER INFO, N 159* .. 160* .. Array Arguments .. 161 DOUBLE PRECISION D( * ), E( * ) 162 COMPLEX*16 AP( * ), TAU( * ) 163* .. 164* 165* ===================================================================== 166* 167* .. Parameters .. 168 COMPLEX*16 ONE, ZERO, HALF 169 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 170 $ ZERO = ( 0.0D+0, 0.0D+0 ), 171 $ HALF = ( 0.5D+0, 0.0D+0 ) ) 172* .. 173* .. Local Scalars .. 174 LOGICAL UPPER 175 INTEGER I, I1, I1I1, II 176 COMPLEX*16 ALPHA, TAUI 177* .. 178* .. External Subroutines .. 179 EXTERNAL XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG 180* .. 181* .. External Functions .. 182 LOGICAL LSAME 183 COMPLEX*16 ZDOTC 184 EXTERNAL LSAME, ZDOTC 185* .. 186* .. Intrinsic Functions .. 187 INTRINSIC DBLE 188* .. 189* .. Executable Statements .. 190* 191* Test the input parameters 192* 193 INFO = 0 194 UPPER = LSAME( UPLO, 'U' ) 195 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 196 INFO = -1 197 ELSE IF( N.LT.0 ) THEN 198 INFO = -2 199 END IF 200 IF( INFO.NE.0 ) THEN 201 CALL XERBLA( 'ZHPTRD', -INFO ) 202 RETURN 203 END IF 204* 205* Quick return if possible 206* 207 IF( N.LE.0 ) 208 $ RETURN 209* 210 IF( UPPER ) THEN 211* 212* Reduce the upper triangle of A. 213* I1 is the index in AP of A(1,I+1). 214* 215 I1 = N*( N-1 ) / 2 + 1 216 AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) ) 217 DO 10 I = N - 1, 1, -1 218* 219* Generate elementary reflector H(i) = I - tau * v * v**H 220* to annihilate A(1:i-1,i+1) 221* 222 ALPHA = AP( I1+I-1 ) 223 CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI ) 224 E( I ) = DBLE( ALPHA ) 225* 226 IF( TAUI.NE.ZERO ) THEN 227* 228* Apply H(i) from both sides to A(1:i,1:i) 229* 230 AP( I1+I-1 ) = ONE 231* 232* Compute y := tau * A * v storing y in TAU(1:i) 233* 234 CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU, 235 $ 1 ) 236* 237* Compute w := y - 1/2 * tau * (y**H *v) * v 238* 239 ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 ) 240 CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 ) 241* 242* Apply the transformation as a rank-2 update: 243* A := A - v * w**H - w * v**H 244* 245 CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP ) 246* 247 END IF 248 AP( I1+I-1 ) = E( I ) 249 D( I+1 ) = DBLE( AP( I1+I ) ) 250 TAU( I ) = TAUI 251 I1 = I1 - I 252 10 CONTINUE 253 D( 1 ) = DBLE( AP( 1 ) ) 254 ELSE 255* 256* Reduce the lower triangle of A. II is the index in AP of 257* A(i,i) and I1I1 is the index of A(i+1,i+1). 258* 259 II = 1 260 AP( 1 ) = DBLE( AP( 1 ) ) 261 DO 20 I = 1, N - 1 262 I1I1 = II + N - I + 1 263* 264* Generate elementary reflector H(i) = I - tau * v * v**H 265* to annihilate A(i+2:n,i) 266* 267 ALPHA = AP( II+1 ) 268 CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI ) 269 E( I ) = DBLE( ALPHA ) 270* 271 IF( TAUI.NE.ZERO ) THEN 272* 273* Apply H(i) from both sides to A(i+1:n,i+1:n) 274* 275 AP( II+1 ) = ONE 276* 277* Compute y := tau * A * v storing y in TAU(i:n-1) 278* 279 CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1, 280 $ ZERO, TAU( I ), 1 ) 281* 282* Compute w := y - 1/2 * tau * (y**H *v) * v 283* 284 ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ), 285 $ 1 ) 286 CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 ) 287* 288* Apply the transformation as a rank-2 update: 289* A := A - v * w**H - w * v**H 290* 291 CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1, 292 $ AP( I1I1 ) ) 293* 294 END IF 295 AP( II+1 ) = E( I ) 296 D( I ) = DBLE( AP( II ) ) 297 TAU( I ) = TAUI 298 II = I1I1 299 20 CONTINUE 300 D( N ) = DBLE( AP( II ) ) 301 END IF 302* 303 RETURN 304* 305* End of ZHPTRD 306* 307 END 308