![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHPGST 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHPGST + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpgst.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpgst.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpgst.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, ITYPE, N 00026 * .. 00027 * .. Array Arguments .. 00028 * COMPLEX*16 AP( * ), BP( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> ZHPGST reduces a complex Hermitian-definite generalized 00038 *> eigenproblem to standard form, using packed storage. 00039 *> 00040 *> If ITYPE = 1, the problem is A*x = lambda*B*x, 00041 *> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H) 00042 *> 00043 *> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or 00044 *> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L. 00045 *> 00046 *> B must have been previously factorized as U**H*U or L*L**H by ZPPTRF. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] ITYPE 00053 *> \verbatim 00054 *> ITYPE is INTEGER 00055 *> = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H); 00056 *> = 2 or 3: compute U*A*U**H or L**H*A*L. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] UPLO 00060 *> \verbatim 00061 *> UPLO is CHARACTER*1 00062 *> = 'U': Upper triangle of A is stored and B is factored as 00063 *> U**H*U; 00064 *> = 'L': Lower triangle of A is stored and B is factored as 00065 *> L*L**H. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] N 00069 *> \verbatim 00070 *> N is INTEGER 00071 *> The order of the matrices A and B. N >= 0. 00072 *> \endverbatim 00073 *> 00074 *> \param[in,out] AP 00075 *> \verbatim 00076 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 00077 *> On entry, the upper or lower triangle of the Hermitian matrix 00078 *> A, packed columnwise in a linear array. The j-th column of A 00079 *> is stored in the array AP as follows: 00080 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00081 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 00082 *> 00083 *> On exit, if INFO = 0, the transformed matrix, stored in the 00084 *> same format as A. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] BP 00088 *> \verbatim 00089 *> BP is COMPLEX*16 array, dimension (N*(N+1)/2) 00090 *> The triangular factor from the Cholesky factorization of B, 00091 *> stored in the same format as A, as returned by ZPPTRF. 00092 *> \endverbatim 00093 *> 00094 *> \param[out] INFO 00095 *> \verbatim 00096 *> INFO is INTEGER 00097 *> = 0: successful exit 00098 *> < 0: if INFO = -i, the i-th argument had an illegal value 00099 *> \endverbatim 00100 * 00101 * Authors: 00102 * ======== 00103 * 00104 *> \author Univ. of Tennessee 00105 *> \author Univ. of California Berkeley 00106 *> \author Univ. of Colorado Denver 00107 *> \author NAG Ltd. 00108 * 00109 *> \date November 2011 00110 * 00111 *> \ingroup complex16OTHERcomputational 00112 * 00113 * ===================================================================== 00114 SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO ) 00115 * 00116 * -- LAPACK computational routine (version 3.4.0) -- 00117 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00118 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00119 * November 2011 00120 * 00121 * .. Scalar Arguments .. 00122 CHARACTER UPLO 00123 INTEGER INFO, ITYPE, N 00124 * .. 00125 * .. Array Arguments .. 00126 COMPLEX*16 AP( * ), BP( * ) 00127 * .. 00128 * 00129 * ===================================================================== 00130 * 00131 * .. Parameters .. 00132 DOUBLE PRECISION ONE, HALF 00133 PARAMETER ( ONE = 1.0D+0, HALF = 0.5D+0 ) 00134 COMPLEX*16 CONE 00135 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 00136 * .. 00137 * .. Local Scalars .. 00138 LOGICAL UPPER 00139 INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK 00140 DOUBLE PRECISION AJJ, AKK, BJJ, BKK 00141 COMPLEX*16 CT 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV, 00145 $ ZTPSV 00146 * .. 00147 * .. Intrinsic Functions .. 00148 INTRINSIC DBLE 00149 * .. 00150 * .. External Functions .. 00151 LOGICAL LSAME 00152 COMPLEX*16 ZDOTC 00153 EXTERNAL LSAME, ZDOTC 00154 * .. 00155 * .. Executable Statements .. 00156 * 00157 * Test the input parameters. 00158 * 00159 INFO = 0 00160 UPPER = LSAME( UPLO, 'U' ) 00161 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00162 INFO = -1 00163 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00164 INFO = -2 00165 ELSE IF( N.LT.0 ) THEN 00166 INFO = -3 00167 END IF 00168 IF( INFO.NE.0 ) THEN 00169 CALL XERBLA( 'ZHPGST', -INFO ) 00170 RETURN 00171 END IF 00172 * 00173 IF( ITYPE.EQ.1 ) THEN 00174 IF( UPPER ) THEN 00175 * 00176 * Compute inv(U**H)*A*inv(U) 00177 * 00178 * J1 and JJ are the indices of A(1,j) and A(j,j) 00179 * 00180 JJ = 0 00181 DO 10 J = 1, N 00182 J1 = JJ + 1 00183 JJ = JJ + J 00184 * 00185 * Compute the j-th column of the upper triangle of A 00186 * 00187 AP( JJ ) = DBLE( AP( JJ ) ) 00188 BJJ = BP( JJ ) 00189 CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J, 00190 $ BP, AP( J1 ), 1 ) 00191 CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE, 00192 $ AP( J1 ), 1 ) 00193 CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 ) 00194 AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ), 00195 $ 1 ) ) / BJJ 00196 10 CONTINUE 00197 ELSE 00198 * 00199 * Compute inv(L)*A*inv(L**H) 00200 * 00201 * KK and K1K1 are the indices of A(k,k) and A(k+1,k+1) 00202 * 00203 KK = 1 00204 DO 20 K = 1, N 00205 K1K1 = KK + N - K + 1 00206 * 00207 * Update the lower triangle of A(k:n,k:n) 00208 * 00209 AKK = AP( KK ) 00210 BKK = BP( KK ) 00211 AKK = AKK / BKK**2 00212 AP( KK ) = AKK 00213 IF( K.LT.N ) THEN 00214 CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 ) 00215 CT = -HALF*AKK 00216 CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 00217 CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1, 00218 $ BP( KK+1 ), 1, AP( K1K1 ) ) 00219 CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 00220 CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K, 00221 $ BP( K1K1 ), AP( KK+1 ), 1 ) 00222 END IF 00223 KK = K1K1 00224 20 CONTINUE 00225 END IF 00226 ELSE 00227 IF( UPPER ) THEN 00228 * 00229 * Compute U*A*U**H 00230 * 00231 * K1 and KK are the indices of A(1,k) and A(k,k) 00232 * 00233 KK = 0 00234 DO 30 K = 1, N 00235 K1 = KK + 1 00236 KK = KK + K 00237 * 00238 * Update the upper triangle of A(1:k,1:k) 00239 * 00240 AKK = AP( KK ) 00241 BKK = BP( KK ) 00242 CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP, 00243 $ AP( K1 ), 1 ) 00244 CT = HALF*AKK 00245 CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 00246 CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1, 00247 $ AP ) 00248 CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 00249 CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 ) 00250 AP( KK ) = AKK*BKK**2 00251 30 CONTINUE 00252 ELSE 00253 * 00254 * Compute L**H *A*L 00255 * 00256 * JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1) 00257 * 00258 JJ = 1 00259 DO 40 J = 1, N 00260 J1J1 = JJ + N - J + 1 00261 * 00262 * Compute the j-th column of the lower triangle of A 00263 * 00264 AJJ = AP( JJ ) 00265 BJJ = BP( JJ ) 00266 AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1, 00267 $ BP( JJ+1 ), 1 ) 00268 CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 ) 00269 CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1, 00270 $ CONE, AP( JJ+1 ), 1 ) 00271 CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit', 00272 $ N-J+1, BP( JJ ), AP( JJ ), 1 ) 00273 JJ = J1J1 00274 40 CONTINUE 00275 END IF 00276 END IF 00277 RETURN 00278 * 00279 * End of ZHPGST 00280 * 00281 END