 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ 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

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```
Date
May 2020

Definition at line 230 of file zlaqz2.f.

234  IMPLICIT NONE
235
236 * Arguments
237  LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
238  INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
239  \$ LDQC, LDZC, LWORK, REC
240
241  COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
242  \$ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
243  INTEGER, INTENT( OUT ) :: NS, ND, INFO
244  COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
245  DOUBLE PRECISION :: RWORK( * )
246
247 * Parameters
248  COMPLEX*16 CZERO, CONE
249  parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
250  \$ 0.0d+0 ) )
251  DOUBLE PRECISION :: ZERO, ONE, HALF
252  parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
253
254 * Local Scalars
255  INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, ZTGEXC_INFO,
256  \$ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
257  DOUBLE PRECISION ::SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
258  COMPLEX*16 :: S, S1, TEMP
259
260 * External Functions
261  EXTERNAL :: xerbla, zlaqz0, zlaqz1, dlabad, zlacpy, zlaset, zgemm,
262  \$ ztgexc, zlartg, zrot
263  DOUBLE PRECISION, EXTERNAL :: DLAMCH
264
265  info = 0
266
267 * Set up deflation window
268  jw = min( nw, ihi-ilo+1 )
269  kwtop = ihi-jw+1
270  IF ( kwtop .EQ. ilo ) THEN
271  s = czero
272  ELSE
273  s = a( kwtop, kwtop-1 )
274  END IF
275
276 * Determine required workspace
277  ifst = 1
278  ilst = jw
279  CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
280  \$ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
281  \$ ldzc, work, -1, rwork, rec+1, qz_small_info )
282  lworkreq = int( work( 1 ) )+2*jw**2
283  lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
284  IF ( lwork .EQ.-1 ) THEN
285 * workspace query, quick return
286  work( 1 ) = lworkreq
287  RETURN
288  ELSE IF ( lwork .LT. lworkreq ) THEN
289  info = -26
290  END IF
291
292  IF( info.NE.0 ) THEN
293  CALL xerbla( 'ZLAQZ2', -info )
294  RETURN
295  END IF
296
297 * Get machine constants
298  safmin = dlamch( 'SAFE MINIMUM' )
299  safmax = one/safmin
300  CALL dlabad( safmin, safmax )
301  ulp = dlamch( 'PRECISION' )
302  smlnum = safmin*( dble( n )/ulp )
303
304  IF ( ihi .EQ. kwtop ) THEN
305 * 1 by 1 deflation window, just try a regular deflation
306  alpha( kwtop ) = a( kwtop, kwtop )
307  beta( kwtop ) = b( kwtop, kwtop )
308  ns = 1
309  nd = 0
310  IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
311  \$ kwtop ) ) ) ) THEN
312  ns = 0
313  nd = 1
314  IF ( kwtop .GT. ilo ) THEN
315  a( kwtop, kwtop-1 ) = czero
316  END IF
317  END IF
318  END IF
319
320
321 * Store window in case of convergence failure
322  CALL zlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
323  CALL zlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
324  \$ 1 ), jw )
325
326 * Transform window to real schur form
327  CALL zlaset( 'FULL', jw, jw, czero, cone, qc, ldqc )
328  CALL zlaset( 'FULL', jw, jw, czero, cone, zc, ldzc )
329  CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
330  \$ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
331  \$ ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
332  \$ rec+1, qz_small_info )
333
334  IF( qz_small_info .NE. 0 ) THEN
335 * Convergence failure, restore the window and exit
336  nd = 0
337  ns = jw-qz_small_info
338  CALL zlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), 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, c1,
396  \$ s1 )
397  CALL zrot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
398  \$ ldb, c1, s1 )
399  CALL zrot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
400  \$ 1, c1, dconjg( s1 ) )
401  END DO
402
403 * Chase bulges down
404  istartm = kwtop
405  istopm = ihi
406  k = kwbot-1
407  DO WHILE ( k .GE. kwtop )
408
409 * Move bulge down and remove it
410  DO k2 = k, kwbot-1
411  CALL zlaqz1( .true., .true., k2, kwtop, kwtop+jw-1,
412  \$ kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
413  \$ jw, kwtop, zc, ldzc )
414  END DO
415
416  k = k-1
417  END DO
418
419  END IF
420
421 * Apply Qc and Zc to rest of the matrix
422  IF ( ilschur ) THEN
423  istartm = 1
424  istopm = n
425  ELSE
426  istartm = ilo
427  istopm = ihi
428  END IF
429
430  IF ( istopm-ihi > 0 ) THEN
431  CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
432  \$ a( kwtop, ihi+1 ), lda, czero, work, jw )
433  CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
434  \$ ihi+1 ), lda )
435  CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
436  \$ b( kwtop, ihi+1 ), ldb, czero, work, jw )
437  CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
438  \$ ihi+1 ), ldb )
439  END IF
440  IF ( ilq ) THEN
441  CALL zgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
442  \$ ldqc, czero, work, n )
443  CALL zlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
444  END IF
445
446  IF ( kwtop-1-istartm+1 > 0 ) THEN
447  CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
448  \$ kwtop ), lda, zc, ldzc, czero, work,
449  \$ kwtop-istartm )
450  CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
451  \$ a( istartm, kwtop ), lda )
452  CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
453  \$ kwtop ), ldb, zc, ldzc, czero, work,
454  \$ kwtop-istartm )
455  CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
456  \$ b( istartm, kwtop ), ldb )
457  END IF
458  IF ( ilz ) THEN
459  CALL zgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
460  \$ ldzc, czero, work, n )
461  CALL zlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
462  END IF
463
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zlaqz1(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
ZLAQZ1
Definition: zlaqz1.f:173
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:284
subroutine ztgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
ZTGEXC
Definition: ztgexc.f:200
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:103
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
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:106
Here is the call graph for this function:
Here is the caller graph for this function: