1 SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, 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, LDA, LDB, N
11 * ..
12 * .. Array Arguments ..
13 COMPLEX*16 A( LDA, * ), B( LDB, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZHEGST reduces a complex Hermitian-definite generalized
20 * eigenproblem to standard form.
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 ZPOTRF.
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 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
47 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
48 * N-by-N upper triangular part of A contains the upper
49 * triangular part of the matrix A, and the strictly lower
50 * triangular part of A is not referenced. If UPLO = 'L', the
51 * leading N-by-N lower triangular part of A contains the lower
52 * triangular part of the matrix A, and the strictly upper
53 * triangular part of A is not referenced.
54 *
55 * On exit, if INFO = 0, the transformed matrix, stored in the
56 * same format as A.
57 *
58 * LDA (input) INTEGER
59 * The leading dimension of the array A. LDA >= max(1,N).
60 *
61 * B (input) COMPLEX*16 array, dimension (LDB,N)
62 * The triangular factor from the Cholesky factorization of B,
63 * as returned by ZPOTRF.
64 *
65 * LDB (input) INTEGER
66 * The leading dimension of the array B. LDB >= max(1,N).
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * < 0: if INFO = -i, the i-th argument had an illegal value
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ONE
76 PARAMETER ( ONE = 1.0D+0 )
77 COMPLEX*16 CONE, HALF
78 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
79 $ HALF = ( 0.5D+0, 0.0D+0 ) )
80 * ..
81 * .. Local Scalars ..
82 LOGICAL UPPER
83 INTEGER K, KB, NB
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL XERBLA, ZHEGS2, ZHEMM, ZHER2K, ZTRMM, ZTRSM
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC MAX, MIN
90 * ..
91 * .. External Functions ..
92 LOGICAL LSAME
93 INTEGER ILAENV
94 EXTERNAL LSAME, ILAENV
95 * ..
96 * .. Executable Statements ..
97 *
98 * Test the input parameters.
99 *
100 INFO = 0
101 UPPER = LSAME( UPLO, 'U' )
102 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
103 INFO = -1
104 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
105 INFO = -2
106 ELSE IF( N.LT.0 ) THEN
107 INFO = -3
108 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
109 INFO = -5
110 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
111 INFO = -7
112 END IF
113 IF( INFO.NE.0 ) THEN
114 CALL XERBLA( 'ZHEGST', -INFO )
115 RETURN
116 END IF
117 *
118 * Quick return if possible
119 *
120 IF( N.EQ.0 )
121 $ RETURN
122 *
123 * Determine the block size for this environment.
124 *
125 NB = ILAENV( 1, 'ZHEGST', UPLO, N, -1, -1, -1 )
126 *
127 IF( NB.LE.1 .OR. NB.GE.N ) THEN
128 *
129 * Use unblocked code
130 *
131 CALL ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
132 ELSE
133 *
134 * Use blocked code
135 *
136 IF( ITYPE.EQ.1 ) THEN
137 IF( UPPER ) THEN
138 *
139 * Compute inv(U**H)*A*inv(U)
140 *
141 DO 10 K = 1, N, NB
142 KB = MIN( N-K+1, NB )
143 *
144 * Update the upper triangle of A(k:n,k:n)
145 *
146 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
147 $ B( K, K ), LDB, INFO )
148 IF( K+KB.LE.N ) THEN
149 CALL ZTRSM( 'Left', UPLO, 'Conjugate transpose',
150 $ 'Non-unit', KB, N-K-KB+1, CONE,
151 $ B( K, K ), LDB, A( K, K+KB ), LDA )
152 CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
153 $ A( K, K ), LDA, B( K, K+KB ), LDB,
154 $ CONE, A( K, K+KB ), LDA )
155 CALL ZHER2K( UPLO, 'Conjugate transpose', N-K-KB+1,
156 $ KB, -CONE, A( K, K+KB ), LDA,
157 $ B( K, K+KB ), LDB, ONE,
158 $ A( K+KB, K+KB ), LDA )
159 CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
160 $ A( K, K ), LDA, B( K, K+KB ), LDB,
161 $ CONE, A( K, K+KB ), LDA )
162 CALL ZTRSM( 'Right', UPLO, 'No transpose',
163 $ 'Non-unit', KB, N-K-KB+1, CONE,
164 $ B( K+KB, K+KB ), LDB, A( K, K+KB ),
165 $ LDA )
166 END IF
167 10 CONTINUE
168 ELSE
169 *
170 * Compute inv(L)*A*inv(L**H)
171 *
172 DO 20 K = 1, N, NB
173 KB = MIN( N-K+1, NB )
174 *
175 * Update the lower triangle of A(k:n,k:n)
176 *
177 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
178 $ B( K, K ), LDB, INFO )
179 IF( K+KB.LE.N ) THEN
180 CALL ZTRSM( 'Right', UPLO, 'Conjugate transpose',
181 $ 'Non-unit', N-K-KB+1, KB, CONE,
182 $ B( K, K ), LDB, A( K+KB, K ), LDA )
183 CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
184 $ A( K, K ), LDA, B( K+KB, K ), LDB,
185 $ CONE, A( K+KB, K ), LDA )
186 CALL ZHER2K( UPLO, 'No transpose', N-K-KB+1, KB,
187 $ -CONE, A( K+KB, K ), LDA,
188 $ B( K+KB, K ), LDB, ONE,
189 $ A( K+KB, K+KB ), LDA )
190 CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
191 $ A( K, K ), LDA, B( K+KB, K ), LDB,
192 $ CONE, A( K+KB, K ), LDA )
193 CALL ZTRSM( 'Left', UPLO, 'No transpose',
194 $ 'Non-unit', N-K-KB+1, KB, CONE,
195 $ B( K+KB, K+KB ), LDB, A( K+KB, K ),
196 $ LDA )
197 END IF
198 20 CONTINUE
199 END IF
200 ELSE
201 IF( UPPER ) THEN
202 *
203 * Compute U*A*U**H
204 *
205 DO 30 K = 1, N, NB
206 KB = MIN( N-K+1, NB )
207 *
208 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
209 *
210 CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
211 $ K-1, KB, CONE, B, LDB, A( 1, K ), LDA )
212 CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
213 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ),
214 $ LDA )
215 CALL ZHER2K( UPLO, 'No transpose', K-1, KB, CONE,
216 $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
217 $ LDA )
218 CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
219 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ),
220 $ LDA )
221 CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose',
222 $ 'Non-unit', K-1, KB, CONE, B( K, K ), LDB,
223 $ A( 1, K ), LDA )
224 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
225 $ B( K, K ), LDB, INFO )
226 30 CONTINUE
227 ELSE
228 *
229 * Compute L**H*A*L
230 *
231 DO 40 K = 1, N, NB
232 KB = MIN( N-K+1, NB )
233 *
234 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
235 *
236 CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
237 $ KB, K-1, CONE, B, LDB, A( K, 1 ), LDA )
238 CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
239 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
240 $ LDA )
241 CALL ZHER2K( UPLO, 'Conjugate transpose', K-1, KB,
242 $ CONE, A( K, 1 ), LDA, B( K, 1 ), LDB,
243 $ ONE, A, LDA )
244 CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
245 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
246 $ LDA )
247 CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose',
248 $ 'Non-unit', KB, K-1, CONE, B( K, K ), LDB,
249 $ A( K, 1 ), LDA )
250 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
251 $ B( K, K ), LDB, INFO )
252 40 CONTINUE
253 END IF
254 END IF
255 END IF
256 RETURN
257 *
258 * End of ZHEGST
259 *
260 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, LDA, LDB, N
11 * ..
12 * .. Array Arguments ..
13 COMPLEX*16 A( LDA, * ), B( LDB, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZHEGST reduces a complex Hermitian-definite generalized
20 * eigenproblem to standard form.
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 ZPOTRF.
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 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
47 * On entry, the Hermitian matrix A. If UPLO = 'U', the leading
48 * N-by-N upper triangular part of A contains the upper
49 * triangular part of the matrix A, and the strictly lower
50 * triangular part of A is not referenced. If UPLO = 'L', the
51 * leading N-by-N lower triangular part of A contains the lower
52 * triangular part of the matrix A, and the strictly upper
53 * triangular part of A is not referenced.
54 *
55 * On exit, if INFO = 0, the transformed matrix, stored in the
56 * same format as A.
57 *
58 * LDA (input) INTEGER
59 * The leading dimension of the array A. LDA >= max(1,N).
60 *
61 * B (input) COMPLEX*16 array, dimension (LDB,N)
62 * The triangular factor from the Cholesky factorization of B,
63 * as returned by ZPOTRF.
64 *
65 * LDB (input) INTEGER
66 * The leading dimension of the array B. LDB >= max(1,N).
67 *
68 * INFO (output) INTEGER
69 * = 0: successful exit
70 * < 0: if INFO = -i, the i-th argument had an illegal value
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ONE
76 PARAMETER ( ONE = 1.0D+0 )
77 COMPLEX*16 CONE, HALF
78 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
79 $ HALF = ( 0.5D+0, 0.0D+0 ) )
80 * ..
81 * .. Local Scalars ..
82 LOGICAL UPPER
83 INTEGER K, KB, NB
84 * ..
85 * .. External Subroutines ..
86 EXTERNAL XERBLA, ZHEGS2, ZHEMM, ZHER2K, ZTRMM, ZTRSM
87 * ..
88 * .. Intrinsic Functions ..
89 INTRINSIC MAX, MIN
90 * ..
91 * .. External Functions ..
92 LOGICAL LSAME
93 INTEGER ILAENV
94 EXTERNAL LSAME, ILAENV
95 * ..
96 * .. Executable Statements ..
97 *
98 * Test the input parameters.
99 *
100 INFO = 0
101 UPPER = LSAME( UPLO, 'U' )
102 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
103 INFO = -1
104 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
105 INFO = -2
106 ELSE IF( N.LT.0 ) THEN
107 INFO = -3
108 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
109 INFO = -5
110 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
111 INFO = -7
112 END IF
113 IF( INFO.NE.0 ) THEN
114 CALL XERBLA( 'ZHEGST', -INFO )
115 RETURN
116 END IF
117 *
118 * Quick return if possible
119 *
120 IF( N.EQ.0 )
121 $ RETURN
122 *
123 * Determine the block size for this environment.
124 *
125 NB = ILAENV( 1, 'ZHEGST', UPLO, N, -1, -1, -1 )
126 *
127 IF( NB.LE.1 .OR. NB.GE.N ) THEN
128 *
129 * Use unblocked code
130 *
131 CALL ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
132 ELSE
133 *
134 * Use blocked code
135 *
136 IF( ITYPE.EQ.1 ) THEN
137 IF( UPPER ) THEN
138 *
139 * Compute inv(U**H)*A*inv(U)
140 *
141 DO 10 K = 1, N, NB
142 KB = MIN( N-K+1, NB )
143 *
144 * Update the upper triangle of A(k:n,k:n)
145 *
146 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
147 $ B( K, K ), LDB, INFO )
148 IF( K+KB.LE.N ) THEN
149 CALL ZTRSM( 'Left', UPLO, 'Conjugate transpose',
150 $ 'Non-unit', KB, N-K-KB+1, CONE,
151 $ B( K, K ), LDB, A( K, K+KB ), LDA )
152 CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
153 $ A( K, K ), LDA, B( K, K+KB ), LDB,
154 $ CONE, A( K, K+KB ), LDA )
155 CALL ZHER2K( UPLO, 'Conjugate transpose', N-K-KB+1,
156 $ KB, -CONE, A( K, K+KB ), LDA,
157 $ B( K, K+KB ), LDB, ONE,
158 $ A( K+KB, K+KB ), LDA )
159 CALL ZHEMM( 'Left', UPLO, KB, N-K-KB+1, -HALF,
160 $ A( K, K ), LDA, B( K, K+KB ), LDB,
161 $ CONE, A( K, K+KB ), LDA )
162 CALL ZTRSM( 'Right', UPLO, 'No transpose',
163 $ 'Non-unit', KB, N-K-KB+1, CONE,
164 $ B( K+KB, K+KB ), LDB, A( K, K+KB ),
165 $ LDA )
166 END IF
167 10 CONTINUE
168 ELSE
169 *
170 * Compute inv(L)*A*inv(L**H)
171 *
172 DO 20 K = 1, N, NB
173 KB = MIN( N-K+1, NB )
174 *
175 * Update the lower triangle of A(k:n,k:n)
176 *
177 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
178 $ B( K, K ), LDB, INFO )
179 IF( K+KB.LE.N ) THEN
180 CALL ZTRSM( 'Right', UPLO, 'Conjugate transpose',
181 $ 'Non-unit', N-K-KB+1, KB, CONE,
182 $ B( K, K ), LDB, A( K+KB, K ), LDA )
183 CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
184 $ A( K, K ), LDA, B( K+KB, K ), LDB,
185 $ CONE, A( K+KB, K ), LDA )
186 CALL ZHER2K( UPLO, 'No transpose', N-K-KB+1, KB,
187 $ -CONE, A( K+KB, K ), LDA,
188 $ B( K+KB, K ), LDB, ONE,
189 $ A( K+KB, K+KB ), LDA )
190 CALL ZHEMM( 'Right', UPLO, N-K-KB+1, KB, -HALF,
191 $ A( K, K ), LDA, B( K+KB, K ), LDB,
192 $ CONE, A( K+KB, K ), LDA )
193 CALL ZTRSM( 'Left', UPLO, 'No transpose',
194 $ 'Non-unit', N-K-KB+1, KB, CONE,
195 $ B( K+KB, K+KB ), LDB, A( K+KB, K ),
196 $ LDA )
197 END IF
198 20 CONTINUE
199 END IF
200 ELSE
201 IF( UPPER ) THEN
202 *
203 * Compute U*A*U**H
204 *
205 DO 30 K = 1, N, NB
206 KB = MIN( N-K+1, NB )
207 *
208 * Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
209 *
210 CALL ZTRMM( 'Left', UPLO, 'No transpose', 'Non-unit',
211 $ K-1, KB, CONE, B, LDB, A( 1, K ), LDA )
212 CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
213 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ),
214 $ LDA )
215 CALL ZHER2K( UPLO, 'No transpose', K-1, KB, CONE,
216 $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
217 $ LDA )
218 CALL ZHEMM( 'Right', UPLO, K-1, KB, HALF, A( K, K ),
219 $ LDA, B( 1, K ), LDB, CONE, A( 1, K ),
220 $ LDA )
221 CALL ZTRMM( 'Right', UPLO, 'Conjugate transpose',
222 $ 'Non-unit', K-1, KB, CONE, B( K, K ), LDB,
223 $ A( 1, K ), LDA )
224 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
225 $ B( K, K ), LDB, INFO )
226 30 CONTINUE
227 ELSE
228 *
229 * Compute L**H*A*L
230 *
231 DO 40 K = 1, N, NB
232 KB = MIN( N-K+1, NB )
233 *
234 * Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
235 *
236 CALL ZTRMM( 'Right', UPLO, 'No transpose', 'Non-unit',
237 $ KB, K-1, CONE, B, LDB, A( K, 1 ), LDA )
238 CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
239 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
240 $ LDA )
241 CALL ZHER2K( UPLO, 'Conjugate transpose', K-1, KB,
242 $ CONE, A( K, 1 ), LDA, B( K, 1 ), LDB,
243 $ ONE, A, LDA )
244 CALL ZHEMM( 'Left', UPLO, KB, K-1, HALF, A( K, K ),
245 $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
246 $ LDA )
247 CALL ZTRMM( 'Left', UPLO, 'Conjugate transpose',
248 $ 'Non-unit', KB, K-1, CONE, B( K, K ), LDB,
249 $ A( K, 1 ), LDA )
250 CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
251 $ B( K, K ), LDB, INFO )
252 40 CONTINUE
253 END IF
254 END IF
255 END IF
256 RETURN
257 *
258 * End of ZHEGST
259 *
260 END