LAPACK  3.10.0
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

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 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 dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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: