Annotation of rpl/lapack/lapack/zhpgst.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: CHARACTER UPLO
! 10: INTEGER INFO, ITYPE, N
! 11: * ..
! 12: * .. Array Arguments ..
! 13: COMPLEX*16 AP( * ), BP( * )
! 14: * ..
! 15: *
! 16: * Purpose
! 17: * =======
! 18: *
! 19: * ZHPGST reduces a complex Hermitian-definite generalized
! 20: * eigenproblem to standard form, using packed storage.
! 21: *
! 22: * If ITYPE = 1, the problem is A*x = lambda*B*x,
! 23: * and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
! 24: *
! 25: * If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
! 26: * B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
! 27: *
! 28: * B must have been previously factorized as U**H*U or L*L**H by ZPPTRF.
! 29: *
! 30: * Arguments
! 31: * =========
! 32: *
! 33: * ITYPE (input) INTEGER
! 34: * = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
! 35: * = 2 or 3: compute U*A*U**H or L**H*A*L.
! 36: *
! 37: * UPLO (input) CHARACTER*1
! 38: * = 'U': Upper triangle of A is stored and B is factored as
! 39: * U**H*U;
! 40: * = 'L': Lower triangle of A is stored and B is factored as
! 41: * L*L**H.
! 42: *
! 43: * N (input) INTEGER
! 44: * The order of the matrices A and B. N >= 0.
! 45: *
! 46: * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2)
! 47: * On entry, the upper or lower triangle of the Hermitian matrix
! 48: * A, packed columnwise in a linear array. The j-th column of A
! 49: * is stored in the array AP as follows:
! 50: * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
! 51: * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
! 52: *
! 53: * On exit, if INFO = 0, the transformed matrix, stored in the
! 54: * same format as A.
! 55: *
! 56: * BP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
! 57: * The triangular factor from the Cholesky factorization of B,
! 58: * stored in the same format as A, as returned by ZPPTRF.
! 59: *
! 60: * INFO (output) INTEGER
! 61: * = 0: successful exit
! 62: * < 0: if INFO = -i, the i-th argument had an illegal value
! 63: *
! 64: * =====================================================================
! 65: *
! 66: * .. Parameters ..
! 67: DOUBLE PRECISION ONE, HALF
! 68: PARAMETER ( ONE = 1.0D+0, HALF = 0.5D+0 )
! 69: COMPLEX*16 CONE
! 70: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
! 71: * ..
! 72: * .. Local Scalars ..
! 73: LOGICAL UPPER
! 74: INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK
! 75: DOUBLE PRECISION AJJ, AKK, BJJ, BKK
! 76: COMPLEX*16 CT
! 77: * ..
! 78: * .. External Subroutines ..
! 79: EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHPMV, ZHPR2, ZTPMV,
! 80: $ ZTPSV
! 81: * ..
! 82: * .. Intrinsic Functions ..
! 83: INTRINSIC DBLE
! 84: * ..
! 85: * .. External Functions ..
! 86: LOGICAL LSAME
! 87: COMPLEX*16 ZDOTC
! 88: EXTERNAL LSAME, ZDOTC
! 89: * ..
! 90: * .. Executable Statements ..
! 91: *
! 92: * Test the input parameters.
! 93: *
! 94: INFO = 0
! 95: UPPER = LSAME( UPLO, 'U' )
! 96: IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
! 97: INFO = -1
! 98: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 99: INFO = -2
! 100: ELSE IF( N.LT.0 ) THEN
! 101: INFO = -3
! 102: END IF
! 103: IF( INFO.NE.0 ) THEN
! 104: CALL XERBLA( 'ZHPGST', -INFO )
! 105: RETURN
! 106: END IF
! 107: *
! 108: IF( ITYPE.EQ.1 ) THEN
! 109: IF( UPPER ) THEN
! 110: *
! 111: * Compute inv(U')*A*inv(U)
! 112: *
! 113: * J1 and JJ are the indices of A(1,j) and A(j,j)
! 114: *
! 115: JJ = 0
! 116: DO 10 J = 1, N
! 117: J1 = JJ + 1
! 118: JJ = JJ + J
! 119: *
! 120: * Compute the j-th column of the upper triangle of A
! 121: *
! 122: AP( JJ ) = DBLE( AP( JJ ) )
! 123: BJJ = BP( JJ )
! 124: CALL ZTPSV( UPLO, 'Conjugate transpose', 'Non-unit', J,
! 125: $ BP, AP( J1 ), 1 )
! 126: CALL ZHPMV( UPLO, J-1, -CONE, AP, BP( J1 ), 1, CONE,
! 127: $ AP( J1 ), 1 )
! 128: CALL ZDSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
! 129: AP( JJ ) = ( AP( JJ )-ZDOTC( J-1, AP( J1 ), 1, BP( J1 ),
! 130: $ 1 ) ) / BJJ
! 131: 10 CONTINUE
! 132: ELSE
! 133: *
! 134: * Compute inv(L)*A*inv(L')
! 135: *
! 136: * KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
! 137: *
! 138: KK = 1
! 139: DO 20 K = 1, N
! 140: K1K1 = KK + N - K + 1
! 141: *
! 142: * Update the lower triangle of A(k:n,k:n)
! 143: *
! 144: AKK = AP( KK )
! 145: BKK = BP( KK )
! 146: AKK = AKK / BKK**2
! 147: AP( KK ) = AKK
! 148: IF( K.LT.N ) THEN
! 149: CALL ZDSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
! 150: CT = -HALF*AKK
! 151: CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
! 152: CALL ZHPR2( UPLO, N-K, -CONE, AP( KK+1 ), 1,
! 153: $ BP( KK+1 ), 1, AP( K1K1 ) )
! 154: CALL ZAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
! 155: CALL ZTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
! 156: $ BP( K1K1 ), AP( KK+1 ), 1 )
! 157: END IF
! 158: KK = K1K1
! 159: 20 CONTINUE
! 160: END IF
! 161: ELSE
! 162: IF( UPPER ) THEN
! 163: *
! 164: * Compute U*A*U'
! 165: *
! 166: * K1 and KK are the indices of A(1,k) and A(k,k)
! 167: *
! 168: KK = 0
! 169: DO 30 K = 1, N
! 170: K1 = KK + 1
! 171: KK = KK + K
! 172: *
! 173: * Update the upper triangle of A(1:k,1:k)
! 174: *
! 175: AKK = AP( KK )
! 176: BKK = BP( KK )
! 177: CALL ZTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
! 178: $ AP( K1 ), 1 )
! 179: CT = HALF*AKK
! 180: CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
! 181: CALL ZHPR2( UPLO, K-1, CONE, AP( K1 ), 1, BP( K1 ), 1,
! 182: $ AP )
! 183: CALL ZAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
! 184: CALL ZDSCAL( K-1, BKK, AP( K1 ), 1 )
! 185: AP( KK ) = AKK*BKK**2
! 186: 30 CONTINUE
! 187: ELSE
! 188: *
! 189: * Compute L'*A*L
! 190: *
! 191: * JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
! 192: *
! 193: JJ = 1
! 194: DO 40 J = 1, N
! 195: J1J1 = JJ + N - J + 1
! 196: *
! 197: * Compute the j-th column of the lower triangle of A
! 198: *
! 199: AJJ = AP( JJ )
! 200: BJJ = BP( JJ )
! 201: AP( JJ ) = AJJ*BJJ + ZDOTC( N-J, AP( JJ+1 ), 1,
! 202: $ BP( JJ+1 ), 1 )
! 203: CALL ZDSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
! 204: CALL ZHPMV( UPLO, N-J, CONE, AP( J1J1 ), BP( JJ+1 ), 1,
! 205: $ CONE, AP( JJ+1 ), 1 )
! 206: CALL ZTPMV( UPLO, 'Conjugate transpose', 'Non-unit',
! 207: $ N-J+1, BP( JJ ), AP( JJ ), 1 )
! 208: JJ = J1J1
! 209: 40 CONTINUE
! 210: END IF
! 211: END IF
! 212: RETURN
! 213: *
! 214: * End of ZHPGST
! 215: *
! 216: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhpgst.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack