1*> \brief \b ZHPGST 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZHPGST + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpgst.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpgst.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpgst.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO ) 22* 23* .. Scalar Arguments .. 24* CHARACTER UPLO 25* INTEGER INFO, ITYPE, N 26* .. 27* .. Array Arguments .. 28* COMPLEX*16 AP( * ), BP( * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> ZHPGST reduces a complex Hermitian-definite generalized 38*> eigenproblem to standard form, using packed storage. 39*> 40*> If ITYPE = 1, the problem is A*x = lambda*B*x, 41*> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H) 42*> 43*> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or 44*> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L. 45*> 46*> B must have been previously factorized as U**H*U or L*L**H by ZPPTRF. 47*> \endverbatim 48* 49* Arguments: 50* ========== 51* 52*> \param[in] ITYPE 53*> \verbatim 54*> ITYPE is INTEGER 55*> = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H); 56*> = 2 or 3: compute U*A*U**H or L**H*A*L. 57*> \endverbatim 58*> 59*> \param[in] UPLO 60*> \verbatim 61*> UPLO is CHARACTER*1 62*> = 'U': Upper triangle of A is stored and B is factored as 63*> U**H*U; 64*> = 'L': Lower triangle of A is stored and B is factored as 65*> L*L**H. 66*> \endverbatim 67*> 68*> \param[in] N 69*> \verbatim 70*> N is INTEGER 71*> The order of the matrices A and B. N >= 0. 72*> \endverbatim 73*> 74*> \param[in,out] AP 75*> \verbatim 76*> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 77*> On entry, the upper or lower triangle of the Hermitian matrix 78*> A, packed columnwise in a linear array. The j-th column of A 79*> is stored in the array AP as follows: 80*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 81*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 82*> 83*> On exit, if INFO = 0, the transformed matrix, stored in the 84*> same format as A. 85*> \endverbatim 86*> 87*> \param[in] BP 88*> \verbatim 89*> BP is COMPLEX*16 array, dimension (N*(N+1)/2) 90*> The triangular factor from the Cholesky factorization of B, 91*> stored in the same format as A, as returned by ZPPTRF. 92*> \endverbatim 93*> 94*> \param[out] INFO 95*> \verbatim 96*> INFO is INTEGER 97*> = 0: successful exit 98*> < 0: if INFO = -i, the i-th argument had an illegal value 99*> \endverbatim 100* 101* Authors: 102* ======== 103* 104*> \author Univ. of Tennessee 105*> \author Univ. of California Berkeley 106*> \author Univ. of Colorado Denver 107*> \author NAG Ltd. 108* 109*> \ingroup complex16OTHERcomputational 110* 111* ===================================================================== 112 SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO ) 113* 114* -- LAPACK computational routine -- 115* -- LAPACK is a software package provided by Univ. of Tennessee, -- 116* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 117* 118* .. Scalar Arguments .. 119 CHARACTER UPLO 120 INTEGER INFO, ITYPE, N 121* .. 122* .. Array Arguments .. 123 COMPLEX*16 AP( * ), BP( * ) 124* .. 125* 126* ===================================================================== 127* 128* .. Parameters .. 129 DOUBLE PRECISION ONE, HALF 130 PARAMETER ( ONE = 1.0D+0, HALF = 0.5D+0 ) 131 COMPLEX*16 CONE 132 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 133* .. 134* .. Local Scalars .. 135 LOGICAL UPPER 136 INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK 137 DOUBLE PRECISION AJJ, AKK, BJJ, BKK 138 COMPLEX*16 CT 139* .. 140* .. External Subroutines .. 141 EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV, 142 $ ZTPSV 143* .. 144* .. Intrinsic Functions .. 145 INTRINSIC DBLE 146* .. 147* .. External Functions .. 148 LOGICAL LSAME 149 COMPLEX*16 ZDOTC 150 EXTERNAL LSAME, ZDOTC 151* .. 152* .. Executable Statements .. 153* 154* Test the input parameters. 155* 156 INFO = 0 157 UPPER = LSAME( UPLO, 'U' ) 158 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 159 INFO = -1 160 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 161 INFO = -2 162 ELSE IF( N.LT.0 ) THEN 163 INFO = -3 164 END IF 165 IF( INFO.NE.0 ) THEN 166 CALL XERBLA( 'ZHPGST', -INFO ) 167 RETURN 168 END IF 169* 170 IF( ITYPE.EQ.1 ) THEN 171 IF( UPPER ) THEN 172* 173* Compute inv(U**H)*A*inv(U) 174* 175* J1 and JJ are the indices of A(1,j) and A(j,j) 176* 177 JJ = 0 178 DO 10 J = 1, N 179 J1 = JJ + 1 180 JJ = JJ + J 181* 182* Compute the j-th column of the upper triangle of A 183* 184 AP( JJ ) = DBLE( AP( JJ ) ) 185 BJJ = DBLE( BP( JJ ) ) 186 CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J, 187 $ BP, AP( J1 ), 1 ) 188 CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE, 189 $ AP( J1 ), 1 ) 190 CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 ) 191 AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ), 192 $ 1 ) ) / BJJ 193 10 CONTINUE 194 ELSE 195* 196* Compute inv(L)*A*inv(L**H) 197* 198* KK and K1K1 are the indices of A(k,k) and A(k+1,k+1) 199* 200 KK = 1 201 DO 20 K = 1, N 202 K1K1 = KK + N - K + 1 203* 204* Update the lower triangle of A(k:n,k:n) 205* 206 AKK = DBLE( AP( KK ) ) 207 BKK = DBLE( BP( KK ) ) 208 AKK = AKK / BKK**2 209 AP( KK ) = AKK 210 IF( K.LT.N ) THEN 211 CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 ) 212 CT = -HALF*AKK 213 CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 214 CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1, 215 $ BP( KK+1 ), 1, AP( K1K1 ) ) 216 CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 217 CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K, 218 $ BP( K1K1 ), AP( KK+1 ), 1 ) 219 END IF 220 KK = K1K1 221 20 CONTINUE 222 END IF 223 ELSE 224 IF( UPPER ) THEN 225* 226* Compute U*A*U**H 227* 228* K1 and KK are the indices of A(1,k) and A(k,k) 229* 230 KK = 0 231 DO 30 K = 1, N 232 K1 = KK + 1 233 KK = KK + K 234* 235* Update the upper triangle of A(1:k,1:k) 236* 237 AKK = DBLE( AP( KK ) ) 238 BKK = DBLE( BP( KK ) ) 239 CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP, 240 $ AP( K1 ), 1 ) 241 CT = HALF*AKK 242 CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 243 CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1, 244 $ AP ) 245 CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 246 CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 ) 247 AP( KK ) = AKK*BKK**2 248 30 CONTINUE 249 ELSE 250* 251* Compute L**H *A*L 252* 253* JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1) 254* 255 JJ = 1 256 DO 40 J = 1, N 257 J1J1 = JJ + N - J + 1 258* 259* Compute the j-th column of the lower triangle of A 260* 261 AJJ = DBLE( AP( JJ ) ) 262 BJJ = DBLE( BP( JJ ) ) 263 AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1, 264 $ BP( JJ+1 ), 1 ) 265 CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 ) 266 CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1, 267 $ CONE, AP( JJ+1 ), 1 ) 268 CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit', 269 $ N-J+1, BP( JJ ), AP( JJ ), 1 ) 270 JJ = J1J1 271 40 CONTINUE 272 END IF 273 END IF 274 RETURN 275* 276* End of ZHPGST 277* 278 END 279