1 SUBROUTINE ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**H)*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**H)
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**H
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**H *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
2 *
3 * -- LAPACK routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**H)*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**H)
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**H
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**H *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