From baba851215b44ac3b60b9248eb02bcce7eb76247 Mon Sep 17 00:00:00 2001 From: jason Date: Tue, 28 Oct 2008 01:38:50 +0000 Subject: Move LAPACK trunk into position. --- SRC/zhpgst.f | 215 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 215 insertions(+) create mode 100644 SRC/zhpgst.f (limited to 'SRC/zhpgst.f') diff --git a/SRC/zhpgst.f b/SRC/zhpgst.f new file mode 100644 index 00000000..2a9fca87 --- /dev/null +++ b/SRC/zhpgst.f @@ -0,0 +1,215 @@ + SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, ITYPE, N +* .. +* .. Array Arguments .. + COMPLEX*16 AP( * ), BP( * ) +* .. +* +* Purpose +* ======= +* +* ZHPGST reduces a complex Hermitian-definite generalized +* eigenproblem to standard form, using packed storage. +* +* If ITYPE = 1, the problem is A*x = lambda*B*x, +* and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H) +* +* If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or +* B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L. +* +* B must have been previously factorized as U**H*U or L*L**H by ZPPTRF. +* +* Arguments +* ========= +* +* ITYPE (input) INTEGER +* = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H); +* = 2 or 3: compute U*A*U**H or L**H*A*L. +* +* UPLO (input) CHARACTER*1 +* = 'U': Upper triangle of A is stored and B is factored as +* U**H*U; +* = 'L': Lower triangle of A is stored and B is factored as +* L*L**H. +* +* N (input) INTEGER +* The order of the matrices A and B. N >= 0. +* +* AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) +* On entry, the upper or lower triangle of the Hermitian matrix +* A, packed columnwise in a linear array. The j-th column of A +* is stored in the array AP as follows: +* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; +* if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. +* +* On exit, if INFO = 0, the transformed matrix, stored in the +* same format as A. +* +* BP (input) COMPLEX*16 array, dimension (N*(N+1)/2) +* The triangular factor from the Cholesky factorization of B, +* stored in the same format as A, as returned by ZPPTRF. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, HALF + PARAMETER ( ONE = 1.0D+0, HALF = 0.5D+0 ) + COMPLEX*16 CONE + PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK + DOUBLE PRECISION AJJ, AKK, BJJ, BKK + COMPLEX*16 CT +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV, + $ ZTPSV +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE +* .. +* .. External Functions .. + LOGICAL LSAME + COMPLEX*16 ZDOTC + EXTERNAL LSAME, ZDOTC +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + UPPER = LSAME( UPLO, 'U' ) + IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN + INFO = -1 + ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -2 + ELSE IF( N.LT.0 ) THEN + INFO = -3 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZHPGST', -INFO ) + RETURN + END IF +* + IF( ITYPE.EQ.1 ) THEN + IF( UPPER ) THEN +* +* Compute inv(U')*A*inv(U) +* +* J1 and JJ are the indices of A(1,j) and A(j,j) +* + JJ = 0 + DO 10 J = 1, N + J1 = JJ + 1 + JJ = JJ + J +* +* Compute the j-th column of the upper triangle of A +* + AP( JJ ) = DBLE( AP( JJ ) ) + BJJ = BP( JJ ) + CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J, + $ BP, AP( J1 ), 1 ) + CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE, + $ AP( J1 ), 1 ) + CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 ) + AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ), + $ 1 ) ) / BJJ + 10 CONTINUE + ELSE +* +* Compute inv(L)*A*inv(L') +* +* KK and K1K1 are the indices of A(k,k) and A(k+1,k+1) +* + KK = 1 + DO 20 K = 1, N + K1K1 = KK + N - K + 1 +* +* Update the lower triangle of A(k:n,k:n) +* + AKK = AP( KK ) + BKK = BP( KK ) + AKK = AKK / BKK**2 + AP( KK ) = AKK + IF( K.LT.N ) THEN + CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 ) + CT = -HALF*AKK + CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) + CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1, + $ BP( KK+1 ), 1, AP( K1K1 ) ) + CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) + CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K, + $ BP( K1K1 ), AP( KK+1 ), 1 ) + END IF + KK = K1K1 + 20 CONTINUE + END IF + ELSE + IF( UPPER ) THEN +* +* Compute U*A*U' +* +* K1 and KK are the indices of A(1,k) and A(k,k) +* + KK = 0 + DO 30 K = 1, N + K1 = KK + 1 + KK = KK + K +* +* Update the upper triangle of A(1:k,1:k) +* + AKK = AP( KK ) + BKK = BP( KK ) + CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP, + $ AP( K1 ), 1 ) + CT = HALF*AKK + CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) + CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1, + $ AP ) + CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) + CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 ) + AP( KK ) = AKK*BKK**2 + 30 CONTINUE + ELSE +* +* Compute L'*A*L +* +* JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1) +* + JJ = 1 + DO 40 J = 1, N + J1J1 = JJ + N - J + 1 +* +* Compute the j-th column of the lower triangle of A +* + AJJ = AP( JJ ) + BJJ = BP( JJ ) + AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1, + $ BP( JJ+1 ), 1 ) + CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 ) + CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1, + $ CONE, AP( JJ+1 ), 1 ) + CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit', + $ N-J+1, BP( JJ ), AP( JJ ), 1 ) + JJ = J1J1 + 40 CONTINUE + END IF + END IF + RETURN +* +* End of ZHPGST +* + END -- cgit v1.2.3