LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zhegst()

subroutine zhegst ( integer itype,
character uplo,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
integer info )

ZHEGST

Download ZHEGST + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZHEGST reduces a complex Hermitian-definite generalized
!> eigenproblem to standard form.
!>
!> 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 ZPOTRF.
!>
Parameters
[in]ITYPE
!>          ITYPE is 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.
!>
[in]UPLO
!>          UPLO is 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.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
!>          N-by-N upper triangular part of A contains the upper
!>          triangular part of the matrix A, and the strictly lower
!>          triangular part of A is not referenced.  If UPLO = 'L', the
!>          leading N-by-N lower triangular part of A contains the lower
!>          triangular part of the matrix A, and the strictly upper
!>          triangular part of A is not referenced.
!>
!>          On exit, if INFO = 0, the transformed matrix, stored in the
!>          same format as A.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is COMPLEX*16 array, dimension (LDB,N)
!>          The triangular factor from the Cholesky factorization of B,
!>          as returned by ZPOTRF.
!>          B is modified by the routine but restored on exit.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 125 of file zhegst.f.

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*
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
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
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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
Here is the call graph for this function:
Here is the caller graph for this function: