LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zhegst.f
Go to the documentation of this file.
1*> \brief \b ZHEGST
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZHEGST + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegst.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegst.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegst.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
20*
21* .. Scalar Arguments ..
22* CHARACTER UPLO
23* INTEGER INFO, ITYPE, LDA, LDB, N
24* ..
25* .. Array Arguments ..
26* COMPLEX*16 A( LDA, * ), B( LDB, * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> ZHEGST reduces a complex Hermitian-definite generalized
36*> eigenproblem to standard form.
37*>
38*> If ITYPE = 1, the problem is A*x = lambda*B*x,
39*> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
40*>
41*> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
42*> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
43*>
44*> B must have been previously factorized as U**H*U or L*L**H by ZPOTRF.
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] ITYPE
51*> \verbatim
52*> ITYPE is INTEGER
53*> = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
54*> = 2 or 3: compute U*A*U**H or L**H*A*L.
55*> \endverbatim
56*>
57*> \param[in] UPLO
58*> \verbatim
59*> UPLO is CHARACTER*1
60*> = 'U': Upper triangle of A is stored and B is factored as
61*> U**H*U;
62*> = 'L': Lower triangle of A is stored and B is factored as
63*> L*L**H.
64*> \endverbatim
65*>
66*> \param[in] N
67*> \verbatim
68*> N is INTEGER
69*> The order of the matrices A and B. N >= 0.
70*> \endverbatim
71*>
72*> \param[in,out] A
73*> \verbatim
74*> A is COMPLEX*16 array, dimension (LDA,N)
75*> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
76*> N-by-N upper triangular part of A contains the upper
77*> triangular part of the matrix A, and the strictly lower
78*> triangular part of A is not referenced. If UPLO = 'L', the
79*> leading N-by-N lower triangular part of A contains the lower
80*> triangular part of the matrix A, and the strictly upper
81*> triangular part of A is not referenced.
82*>
83*> On exit, if INFO = 0, the transformed matrix, stored in the
84*> same format as A.
85*> \endverbatim
86*>
87*> \param[in] LDA
88*> \verbatim
89*> LDA is INTEGER
90*> The leading dimension of the array A. LDA >= max(1,N).
91*> \endverbatim
92*>
93*> \param[in,out] B
94*> \verbatim
95*> B is COMPLEX*16 array, dimension (LDB,N)
96*> The triangular factor from the Cholesky factorization of B,
97*> as returned by ZPOTRF.
98*> B is modified by the routine but restored on exit.
99*> \endverbatim
100*>
101*> \param[in] LDB
102*> \verbatim
103*> LDB is INTEGER
104*> The leading dimension of the array B. LDB >= max(1,N).
105*> \endverbatim
106*>
107*> \param[out] INFO
108*> \verbatim
109*> INFO is INTEGER
110*> = 0: successful exit
111*> < 0: if INFO = -i, the i-th argument had an illegal value
112*> \endverbatim
113*
114* Authors:
115* ========
116*
117*> \author Univ. of Tennessee
118*> \author Univ. of California Berkeley
119*> \author Univ. of Colorado Denver
120*> \author NAG Ltd.
121*
122*> \ingroup hegst
123*
124* =====================================================================
125 SUBROUTINE zhegst( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
126*
127* -- LAPACK computational routine --
128* -- LAPACK is a software package provided by Univ. of Tennessee, --
129* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130*
131* .. Scalar Arguments ..
132 CHARACTER UPLO
133 INTEGER INFO, ITYPE, LDA, LDB, N
134* ..
135* .. Array Arguments ..
136 COMPLEX*16 A( LDA, * ), B( LDB, * )
137* ..
138*
139* =====================================================================
140*
141* .. Parameters ..
142 DOUBLE PRECISION ONE
143 parameter( one = 1.0d+0 )
144 COMPLEX*16 CONE, HALF
145 parameter( cone = ( 1.0d+0, 0.0d+0 ),
146 $ half = ( 0.5d+0, 0.0d+0 ) )
147* ..
148* .. Local Scalars ..
149 LOGICAL UPPER
150 INTEGER K, KB, NB
151* ..
152* .. External Subroutines ..
153 EXTERNAL xerbla, zhegs2, zhemm, zher2k, ztrmm,
154 $ ztrsm
155* ..
156* .. Intrinsic Functions ..
157 INTRINSIC max, min
158* ..
159* .. External Functions ..
160 LOGICAL LSAME
161 INTEGER ILAENV
162 EXTERNAL lsame, ilaenv
163* ..
164* .. Executable Statements ..
165*
166* Test the input parameters.
167*
168 info = 0
169 upper = lsame( uplo, 'U' )
170 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
171 info = -1
172 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
173 info = -2
174 ELSE IF( n.LT.0 ) THEN
175 info = -3
176 ELSE IF( lda.LT.max( 1, n ) ) THEN
177 info = -5
178 ELSE IF( ldb.LT.max( 1, n ) ) THEN
179 info = -7
180 END IF
181 IF( info.NE.0 ) THEN
182 CALL xerbla( 'ZHEGST', -info )
183 RETURN
184 END IF
185*
186* Quick return if possible
187*
188 IF( n.EQ.0 )
189 $ RETURN
190*
191* Determine the block size for this environment.
192*
193 nb = ilaenv( 1, 'ZHEGST', uplo, n, -1, -1, -1 )
194*
195 IF( nb.LE.1 .OR. nb.GE.n ) THEN
196*
197* Use unblocked code
198*
199 CALL zhegs2( itype, uplo, n, a, lda, b, ldb, info )
200 ELSE
201*
202* Use blocked code
203*
204 IF( itype.EQ.1 ) THEN
205 IF( upper ) THEN
206*
207* Compute inv(U**H)*A*inv(U)
208*
209 DO 10 k = 1, n, nb
210 kb = min( n-k+1, nb )
211*
212* Update the upper triangle of A(k:n,k:n)
213*
214 CALL zhegs2( itype, uplo, kb, a( k, k ), lda,
215 $ b( k, k ), ldb, info )
216 IF( k+kb.LE.n ) THEN
217 CALL ztrsm( 'Left', uplo, 'Conjugate transpose',
218 $ 'Non-unit', kb, n-k-kb+1, cone,
219 $ b( k, k ), ldb, a( k, k+kb ), lda )
220 CALL zhemm( 'Left', uplo, kb, n-k-kb+1, -half,
221 $ a( k, k ), lda, b( k, k+kb ), ldb,
222 $ cone, a( k, k+kb ), lda )
223 CALL zher2k( uplo, 'Conjugate transpose',
224 $ n-k-kb+1,
225 $ kb, -cone, a( k, k+kb ), lda,
226 $ b( k, k+kb ), ldb, one,
227 $ a( k+kb, k+kb ), lda )
228 CALL zhemm( 'Left', uplo, kb, n-k-kb+1, -half,
229 $ a( k, k ), lda, b( k, k+kb ), ldb,
230 $ cone, a( k, k+kb ), lda )
231 CALL ztrsm( 'Right', uplo, 'No transpose',
232 $ 'Non-unit', kb, n-k-kb+1, cone,
233 $ b( k+kb, k+kb ), ldb, a( k, k+kb ),
234 $ lda )
235 END IF
236 10 CONTINUE
237 ELSE
238*
239* Compute inv(L)*A*inv(L**H)
240*
241 DO 20 k = 1, n, nb
242 kb = min( n-k+1, nb )
243*
244* Update the lower triangle of A(k:n,k:n)
245*
246 CALL zhegs2( itype, uplo, kb, a( k, k ), lda,
247 $ b( k, k ), ldb, info )
248 IF( k+kb.LE.n ) THEN
249 CALL ztrsm( 'Right', uplo,
250 $ 'Conjugate transpose',
251 $ 'Non-unit', n-k-kb+1, kb, cone,
252 $ b( k, k ), ldb, a( k+kb, k ), lda )
253 CALL zhemm( 'Right', uplo, n-k-kb+1, kb, -half,
254 $ a( k, k ), lda, b( k+kb, k ), ldb,
255 $ cone, a( k+kb, k ), lda )
256 CALL zher2k( uplo, 'No transpose', n-k-kb+1, kb,
257 $ -cone, a( k+kb, k ), lda,
258 $ b( k+kb, k ), ldb, one,
259 $ a( k+kb, k+kb ), lda )
260 CALL zhemm( 'Right', uplo, n-k-kb+1, kb, -half,
261 $ a( k, k ), lda, b( k+kb, k ), ldb,
262 $ cone, a( k+kb, k ), lda )
263 CALL ztrsm( 'Left', uplo, 'No transpose',
264 $ 'Non-unit', n-k-kb+1, kb, cone,
265 $ b( k+kb, k+kb ), ldb, a( k+kb, k ),
266 $ lda )
267 END IF
268 20 CONTINUE
269 END IF
270 ELSE
271 IF( upper ) THEN
272*
273* Compute U*A*U**H
274*
275 DO 30 k = 1, n, nb
276 kb = min( n-k+1, nb )
277*
278* Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
279*
280 CALL ztrmm( 'Left', uplo, 'No transpose',
281 $ 'Non-unit',
282 $ k-1, kb, cone, b, ldb, a( 1, k ), lda )
283 CALL zhemm( 'Right', uplo, k-1, kb, half, a( k,
284 $ k ),
285 $ lda, b( 1, k ), ldb, cone, a( 1, k ),
286 $ lda )
287 CALL zher2k( uplo, 'No transpose', k-1, kb, cone,
288 $ a( 1, k ), lda, b( 1, k ), ldb, one, a,
289 $ lda )
290 CALL zhemm( 'Right', uplo, k-1, kb, half, a( k,
291 $ k ),
292 $ lda, b( 1, k ), ldb, cone, a( 1, k ),
293 $ lda )
294 CALL ztrmm( 'Right', uplo, 'Conjugate transpose',
295 $ 'Non-unit', k-1, kb, cone, b( k, k ), ldb,
296 $ a( 1, k ), lda )
297 CALL zhegs2( itype, uplo, kb, a( k, k ), lda,
298 $ b( k, k ), ldb, info )
299 30 CONTINUE
300 ELSE
301*
302* Compute L**H*A*L
303*
304 DO 40 k = 1, n, nb
305 kb = min( n-k+1, nb )
306*
307* Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
308*
309 CALL ztrmm( 'Right', uplo, 'No transpose',
310 $ 'Non-unit',
311 $ kb, k-1, cone, b, ldb, a( k, 1 ), lda )
312 CALL zhemm( 'Left', uplo, kb, k-1, half, a( k, k ),
313 $ lda, b( k, 1 ), ldb, cone, a( k, 1 ),
314 $ lda )
315 CALL zher2k( uplo, 'Conjugate transpose', k-1, kb,
316 $ cone, a( k, 1 ), lda, b( k, 1 ), ldb,
317 $ one, a, lda )
318 CALL zhemm( 'Left', uplo, kb, k-1, half, a( k, k ),
319 $ lda, b( k, 1 ), ldb, cone, a( k, 1 ),
320 $ lda )
321 CALL ztrmm( 'Left', uplo, 'Conjugate transpose',
322 $ 'Non-unit', kb, k-1, cone, b( k, k ), ldb,
323 $ a( k, 1 ), lda )
324 CALL zhegs2( itype, uplo, kb, a( k, k ), lda,
325 $ b( k, k ), ldb, info )
326 40 CONTINUE
327 END IF
328 END IF
329 END IF
330 RETURN
331*
332* End of ZHEGST
333*
334 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zhegs2
subroutine zhegs2(itype, uplo, n, a, lda, b, ldb, info)
ZHEGS2 reduces a Hermitian definite generalized eigenproblem to standard form, using the factorizatio...
Definition zhegs2.f:126
zhegst
subroutine zhegst(itype, uplo, n, a, lda, b, ldb, info)
ZHEGST
Definition zhegst.f:126
zhemm
subroutine zhemm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
ZHEMM
Definition zhemm.f:191
zher2k
subroutine zher2k(uplo, trans, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZHER2K
Definition zher2k.f:198
ztrmm
subroutine ztrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRMM
Definition ztrmm.f:177
ztrsm
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180