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*> \date December 2016 110* 111*> \ingroup complex16OTHERcomputational 112* 113* ===================================================================== 114 SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO ) 115* 116* -- LAPACK computational routine (version 3.7.0) -- 117* -- LAPACK is a software package provided by Univ. of Tennessee, -- 118* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 119* December 2016 120* 121* .. Scalar Arguments .. 122 CHARACTER UPLO 123 INTEGER INFO, ITYPE, N 124* .. 125* .. Array Arguments .. 126 COMPLEX*16 AP( * ), BP( * ) 127* .. 128* 129* ===================================================================== 130* 131* .. Parameters .. 132 DOUBLE PRECISION ONE, HALF 133 PARAMETER ( ONE = 1.0D+0, HALF = 0.5D+0 ) 134 COMPLEX*16 CONE 135 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 136* .. 137* .. Local Scalars .. 138 LOGICAL UPPER 139 INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK 140 DOUBLE PRECISION AJJ, AKK, BJJ, BKK 141 COMPLEX*16 CT 142* .. 143* .. External Subroutines .. 144 EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV, 145 $ ZTPSV 146* .. 147* .. Intrinsic Functions .. 148 INTRINSIC DBLE 149* .. 150* .. External Functions .. 151 LOGICAL LSAME 152 COMPLEX*16 ZDOTC 153 EXTERNAL LSAME, ZDOTC 154* .. 155* .. Executable Statements .. 156* 157* Test the input parameters. 158* 159 INFO = 0 160 UPPER = LSAME( UPLO, 'U' ) 161 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 162 INFO = -1 163 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 164 INFO = -2 165 ELSE IF( N.LT.0 ) THEN 166 INFO = -3 167 END IF 168 IF( INFO.NE.0 ) THEN 169 CALL XERBLA( 'ZHPGST', -INFO ) 170 RETURN 171 END IF 172* 173 IF( ITYPE.EQ.1 ) THEN 174 IF( UPPER ) THEN 175* 176* Compute inv(U**H)*A*inv(U) 177* 178* J1 and JJ are the indices of A(1,j) and A(j,j) 179* 180 JJ = 0 181 DO 10 J = 1, N 182 J1 = JJ + 1 183 JJ = JJ + J 184* 185* Compute the j-th column of the upper triangle of A 186* 187 AP( JJ ) = DBLE( AP( JJ ) ) 188 BJJ = BP( JJ ) 189 CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J, 190 $ BP, AP( J1 ), 1 ) 191 CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE, 192 $ AP( J1 ), 1 ) 193 CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 ) 194 AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ), 195 $ 1 ) ) / BJJ 196 10 CONTINUE 197 ELSE 198* 199* Compute inv(L)*A*inv(L**H) 200* 201* KK and K1K1 are the indices of A(k,k) and A(k+1,k+1) 202* 203 KK = 1 204 DO 20 K = 1, N 205 K1K1 = KK + N - K + 1 206* 207* Update the lower triangle of A(k:n,k:n) 208* 209 AKK = AP( KK ) 210 BKK = BP( KK ) 211 AKK = AKK / BKK**2 212 AP( KK ) = AKK 213 IF( K.LT.N ) THEN 214 CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 ) 215 CT = -HALF*AKK 216 CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 217 CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1, 218 $ BP( KK+1 ), 1, AP( K1K1 ) ) 219 CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 220 CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K, 221 $ BP( K1K1 ), AP( KK+1 ), 1 ) 222 END IF 223 KK = K1K1 224 20 CONTINUE 225 END IF 226 ELSE 227 IF( UPPER ) THEN 228* 229* Compute U*A*U**H 230* 231* K1 and KK are the indices of A(1,k) and A(k,k) 232* 233 KK = 0 234 DO 30 K = 1, N 235 K1 = KK + 1 236 KK = KK + K 237* 238* Update the upper triangle of A(1:k,1:k) 239* 240 AKK = AP( KK ) 241 BKK = BP( KK ) 242 CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP, 243 $ AP( K1 ), 1 ) 244 CT = HALF*AKK 245 CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 246 CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1, 247 $ AP ) 248 CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 249 CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 ) 250 AP( KK ) = AKK*BKK**2 251 30 CONTINUE 252 ELSE 253* 254* Compute L**H *A*L 255* 256* JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1) 257* 258 JJ = 1 259 DO 40 J = 1, N 260 J1J1 = JJ + N - J + 1 261* 262* Compute the j-th column of the lower triangle of A 263* 264 AJJ = AP( JJ ) 265 BJJ = BP( JJ ) 266 AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1, 267 $ BP( JJ+1 ), 1 ) 268 CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 ) 269 CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1, 270 $ CONE, AP( JJ+1 ), 1 ) 271 CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit', 272 $ N-J+1, BP( JJ ), AP( JJ ), 1 ) 273 JJ = J1J1 274 40 CONTINUE 275 END IF 276 END IF 277 RETURN 278* 279* End of ZHPGST 280* 281 END 282