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

◆ zlaqz2()

recursive subroutine zlaqz2 ( logical, intent(in) ilschur,
logical, intent(in) ilq,
logical, intent(in) ilz,
integer, intent(in) n,
integer, intent(in) ilo,
integer, intent(in) ihi,
integer, intent(in) nw,
complex*16, dimension( lda, * ), intent(inout) a,
integer, intent(in) lda,
complex*16, dimension( ldb, * ), intent(inout) b,
integer, intent(in) ldb,
complex*16, dimension( ldq, * ), intent(inout) q,
integer, intent(in) ldq,
complex*16, dimension( ldz, * ), intent(inout) z,
integer, intent(in) ldz,
integer, intent(out) ns,
integer, intent(out) nd,
complex*16, dimension( * ), intent(inout) alpha,
complex*16, dimension( * ), intent(inout) beta,
complex*16, dimension( ldqc, * ) qc,
integer, intent(in) ldqc,
complex*16, dimension( ldzc, * ) zc,
integer, intent(in) ldzc,
complex*16, dimension( * ) work,
integer, intent(in) lwork,
double precision, dimension( * ) rwork,
integer, intent(in) rec,
integer, intent(out) info )

ZLAQZ2

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

Purpose:
!>
!> ZLAQZ2 performs AED
!>
Parameters
[in]ILSCHUR
!>          ILSCHUR is LOGICAL
!>              Determines whether or not to update the full Schur form
!>
[in]ILQ
!>          ILQ is LOGICAL
!>              Determines whether or not to update the matrix Q
!>
[in]ILZ
!>          ILZ is LOGICAL
!>              Determines whether or not to update the matrix Z
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A, B, Q, and Z.  N >= 0.
!>
[in]ILO
!>          ILO is INTEGER
!>
[in]IHI
!>          IHI is INTEGER
!>          ILO and IHI mark the rows and columns of (A,B) which
!>          are to be normalized
!>
[in]NW
!>          NW is INTEGER
!>          The desired size of the deflation window.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA, N)
!>
[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)
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max( 1, N ).
!>
[in,out]Q
!>          Q is COMPLEX*16 array, dimension (LDQ, N)
!>
[in]LDQ
!>          LDQ is INTEGER
!>
[in,out]Z
!>          Z is COMPLEX*16 array, dimension (LDZ, N)
!>
[in]LDZ
!>          LDZ is INTEGER
!>
[out]NS
!>          NS is INTEGER
!>          The number of unconverged eigenvalues available to
!>          use as shifts.
!>
[out]ND
!>          ND is INTEGER
!>          The number of converged eigenvalues found.
!>
[out]ALPHA
!>          ALPHA is COMPLEX*16 array, dimension (N)
!>          Each scalar alpha defining an eigenvalue
!>          of GNEP.
!>
[out]BETA
!>          BETA is COMPLEX*16 array, dimension (N)
!>          The scalars beta that define the eigenvalues of GNEP.
!>          Together, the quantities alpha = ALPHA(j) and
!>          beta = BETA(j) represent the j-th eigenvalue of the matrix
!>          pair (A,B), in one of the forms lambda = alpha/beta or
!>          mu = beta/alpha.  Since either lambda or mu may overflow,
!>          they should not, in general, be computed.
!>
[in,out]QC
!>          QC is COMPLEX*16 array, dimension (LDQC, NW)
!>
[in]LDQC
!>          LDQC is INTEGER
!>
[in,out]ZC
!>          ZC is COMPLEX*16 array, dimension (LDZC, NW)
!>
[in]LDZC
!>          LDZ is INTEGER
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,N).
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!>
[in]REC
!>          REC is INTEGER
!>             REC indicates the current recursion level. Should be set
!>             to 0 on first call.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -i, the i-th argument had an illegal value
!>
Author
Thijs Steel, KU Leuven
Date
May 2020

Definition at line 228 of file zlaqz2.f.

233 IMPLICIT NONE
234
235* Arguments
236 LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
237 INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
238 $ LDQC, LDZC, LWORK, REC
239
240 COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
241 $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
242 INTEGER, INTENT( OUT ) :: NS, ND, INFO
243 COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
244 DOUBLE PRECISION :: RWORK( * )
245
246* Parameters
247 COMPLEX*16 CZERO, CONE
248 parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
249 $ 0.0d+0 ) )
250 DOUBLE PRECISION :: ZERO, ONE, HALF
251 parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
252
253* Local Scalars
254 INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, ZTGEXC_INFO,
255 $ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
256 DOUBLE PRECISION ::SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
257 COMPLEX*16 :: S, S1, TEMP
258
259* External Functions
260 EXTERNAL :: xerbla, zlaqz0, zlaqz1, zlacpy, zlaset, zgemm,
261 $ ztgexc, zlartg, zrot
262 DOUBLE PRECISION, EXTERNAL :: DLAMCH
263
264 info = 0
265
266* Set up deflation window
267 jw = min( nw, ihi-ilo+1 )
268 kwtop = ihi-jw+1
269 IF ( kwtop .EQ. ilo ) THEN
270 s = czero
271 ELSE
272 s = a( kwtop, kwtop-1 )
273 END IF
274
275* Determine required workspace
276 ifst = 1
277 ilst = jw
278 CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
279 $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
280 $ ldzc, work, -1, rwork, rec+1, qz_small_info )
281 lworkreq = int( work( 1 ) )+2*jw**2
282 lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
283 IF ( lwork .EQ.-1 ) THEN
284* workspace query, quick return
285 work( 1 ) = lworkreq
286 RETURN
287 ELSE IF ( lwork .LT. lworkreq ) THEN
288 info = -26
289 END IF
290
291 IF( info.NE.0 ) THEN
292 CALL xerbla( 'ZLAQZ2', -info )
293 RETURN
294 END IF
295
296* Get machine constants
297 safmin = dlamch( 'SAFE MINIMUM' )
298 safmax = one/safmin
299 ulp = dlamch( 'PRECISION' )
300 smlnum = safmin*( dble( n )/ulp )
301
302 IF ( ihi .EQ. kwtop ) THEN
303* 1 by 1 deflation window, just try a regular deflation
304 alpha( kwtop ) = a( kwtop, kwtop )
305 beta( kwtop ) = b( kwtop, kwtop )
306 ns = 1
307 nd = 0
308 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
309 $ kwtop ) ) ) ) THEN
310 ns = 0
311 nd = 1
312 IF ( kwtop .GT. ilo ) THEN
313 a( kwtop, kwtop-1 ) = czero
314 END IF
315 END IF
316 END IF
317
318
319* Store window in case of convergence failure
320 CALL zlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
321 CALL zlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb,
322 $ work( jw**2+
323 $ 1 ), jw )
324
325* Transform window to real schur form
326 CALL zlaset( 'FULL', jw, jw, czero, cone, qc, ldqc )
327 CALL zlaset( 'FULL', jw, jw, czero, cone, zc, ldzc )
328 CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
329 $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
330 $ ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
331 $ rec+1, qz_small_info )
332
333 IF( qz_small_info .NE. 0 ) THEN
334* Convergence failure, restore the window and exit
335 nd = 0
336 ns = jw-qz_small_info
337 CALL zlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ),
338 $ lda )
339 CALL zlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
340 $ kwtop ), ldb )
341 RETURN
342 END IF
343
344* Deflation detection loop
345 IF ( kwtop .EQ. ilo .OR. s .EQ. czero ) THEN
346 kwbot = kwtop-1
347 ELSE
348 kwbot = ihi
349 k = 1
350 k2 = 1
351 DO WHILE ( k .LE. jw )
352* Try to deflate eigenvalue
353 tempr = abs( a( kwbot, kwbot ) )
354 IF( tempr .EQ. zero ) THEN
355 tempr = abs( s )
356 END IF
357 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
358 $ tempr, smlnum ) ) THEN
359* Deflatable
360 kwbot = kwbot-1
361 ELSE
362* Not deflatable, move out of the way
363 ifst = kwbot-kwtop+1
364 ilst = k2
365 CALL ztgexc( .true., .true., jw, a( kwtop, kwtop ),
366 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
367 $ zc, ldzc, ifst, ilst, ztgexc_info )
368 k2 = k2+1
369 END IF
370
371 k = k+1
372 END DO
373 END IF
374
375* Store eigenvalues
376 nd = ihi-kwbot
377 ns = jw-nd
378 k = kwtop
379 DO WHILE ( k .LE. ihi )
380 alpha( k ) = a( k, k )
381 beta( k ) = b( k, k )
382 k = k+1
383 END DO
384
385 IF ( kwtop .NE. ilo .AND. s .NE. czero ) THEN
386* Reflect spike back, this will create optimally packed bulges
387 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 ) *dconjg( qc( 1,
388 $ 1:jw-nd ) )
389 DO k = kwbot-1, kwtop, -1
390 CALL zlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
391 $ temp )
392 a( k, kwtop-1 ) = temp
393 a( k+1, kwtop-1 ) = czero
394 k2 = max( kwtop, k-1 )
395 CALL zrot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda,
396 $ c1,
397 $ s1 )
398 CALL zrot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1,
399 $ k-1 ),
400 $ ldb, c1, s1 )
401 CALL zrot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1,
402 $ k+1-kwtop+1 ),
403 $ 1, c1, dconjg( s1 ) )
404 END DO
405
406* Chase bulges down
407 istartm = kwtop
408 istopm = ihi
409 k = kwbot-1
410 DO WHILE ( k .GE. kwtop )
411
412* Move bulge down and remove it
413 DO k2 = k, kwbot-1
414 CALL zlaqz1( .true., .true., k2, kwtop, kwtop+jw-1,
415 $ kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
416 $ jw, kwtop, zc, ldzc )
417 END DO
418
419 k = k-1
420 END DO
421
422 END IF
423
424* Apply Qc and Zc to rest of the matrix
425 IF ( ilschur ) THEN
426 istartm = 1
427 istopm = n
428 ELSE
429 istartm = ilo
430 istopm = ihi
431 END IF
432
433 IF ( istopm-ihi > 0 ) THEN
434 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
435 $ a( kwtop, ihi+1 ), lda, czero, work, jw )
436 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
437 $ ihi+1 ), lda )
438 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
439 $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
440 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
441 $ ihi+1 ), ldb )
442 END IF
443 IF ( ilq ) THEN
444 CALL zgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq,
445 $ qc,
446 $ ldqc, czero, work, n )
447 CALL zlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
448 END IF
449
450 IF ( kwtop-1-istartm+1 > 0 ) THEN
451 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone,
452 $ a( istartm,
453 $ kwtop ), lda, zc, ldzc, czero, work,
454 $ kwtop-istartm )
455 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
456 $ a( istartm, kwtop ), lda )
457 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone,
458 $ b( istartm,
459 $ kwtop ), ldb, zc, ldzc, czero, work,
460 $ kwtop-istartm )
461 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
462 $ b( istartm, kwtop ), ldb )
463 END IF
464 IF ( ilz ) THEN
465 CALL zgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz,
466 $ zc,
467 $ ldzc, czero, work, n )
468 CALL zlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
469 END IF
470
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlaqz0
recursive subroutine zlaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, work, lwork, rwork, rec, info)
ZLAQZ0
Definition zlaqz0.f:283
zlaqz1
subroutine zlaqz1(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
ZLAQZ1
Definition zlaqz1.f:172
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:104
zrot
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition zrot.f:101
ztgexc
subroutine ztgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, info)
ZTGEXC
Definition ztgexc.f:198
Here is the call graph for this function:
Here is the caller graph for this function: