LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlaqz3()

recursive subroutine dlaqz3 ( 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,
double precision, dimension( lda, * ), intent(inout)  A,
integer, intent(in)  LDA,
double precision, dimension( ldb, * ), intent(inout)  B,
integer, intent(in)  LDB,
double precision, dimension( ldq, * ), intent(inout)  Q,
integer, intent(in)  LDQ,
double precision, dimension( ldz, * ), intent(inout)  Z,
integer, intent(in)  LDZ,
integer, intent(out)  NS,
integer, intent(out)  ND,
double precision, dimension( * ), intent(inout)  ALPHAR,
double precision, dimension( * ), intent(inout)  ALPHAI,
double precision, dimension( * ), intent(inout)  BETA,
double precision, dimension( ldqc, * )  QC,
integer, intent(in)  LDQC,
double precision, dimension( ldzc, * )  ZC,
integer, intent(in)  LDZC,
double precision, dimension( * )  WORK,
integer, intent(in)  LWORK,
integer, intent(in)  REC,
integer, intent(out)  INFO 
)

DLAQZ3

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

Purpose:
 DLAQZ3 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDQ, N)
[in]LDQ
          LDQ is INTEGER
[in,out]Z
          Z is DOUBLE PRECISION 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]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (N)
          The real parts of each scalar alpha defining an eigenvalue
          of GNEP.
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (N)
          The imaginary parts of each scalar alpha defining an
          eigenvalue of GNEP.
          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
          positive, then the j-th and (j+1)-st eigenvalues are a
          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (N)
          The scalars beta that define the eigenvalues of GNEP.
          Together, the quantities alpha = (ALPHAR(j),ALPHAI(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 DOUBLE PRECISION array, dimension (LDQC, NW)
[in]LDQC
          LDQC is INTEGER
[in,out]ZC
          ZC is DOUBLE PRECISION array, dimension (LDZC, NW)
[in]LDZC
          LDZ is INTEGER
[out]WORK
          WORK is DOUBLE PRECISION 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.
[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 235 of file dlaqz3.f.

239  IMPLICIT NONE
240 
241 * Arguments
242  LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
243  INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
244  $ LDQC, LDZC, LWORK, REC
245 
246  DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
247  $ Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ),
248  $ ALPHAI( * ), BETA( * )
249  INTEGER, INTENT( OUT ) :: NS, ND, INFO
250  DOUBLE PRECISION :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
251 
252 * Parameters
253  DOUBLE PRECISION :: ZERO, ONE, HALF
254  parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
255 
256 * Local Scalars
257  LOGICAL :: BULGE
258  INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, DTGEXC_INFO,
259  $ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
260  DOUBLE PRECISION :: S, SMLNUM, ULP, SAFMIN, SAFMAX, C1, S1, TEMP
261 
262 * External Functions
263  EXTERNAL :: xerbla, dtgexc, dlabad, dlaqz0, dlacpy, dlaset,
265  DOUBLE PRECISION, EXTERNAL :: DLAMCH
266 
267  info = 0
268 
269 * Set up deflation window
270  jw = min( nw, ihi-ilo+1 )
271  kwtop = ihi-jw+1
272  IF ( kwtop .EQ. ilo ) THEN
273  s = zero
274  ELSE
275  s = a( kwtop, kwtop-1 )
276  END IF
277 
278 * Determine required workspace
279  ifst = 1
280  ilst = jw
281  CALL dtgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
282  $ ldzc, ifst, ilst, work, -1, dtgexc_info )
283  lworkreq = int( work( 1 ) )
284  CALL dlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
285  $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
286  $ ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
287  lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
288  lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
289  IF ( lwork .EQ.-1 ) THEN
290 * workspace query, quick return
291  work( 1 ) = lworkreq
292  RETURN
293  ELSE IF ( lwork .LT. lworkreq ) THEN
294  info = -26
295  END IF
296 
297  IF( info.NE.0 ) THEN
298  CALL xerbla( 'DLAQZ3', -info )
299  RETURN
300  END IF
301 
302 * Get machine constants
303  safmin = dlamch( 'SAFE MINIMUM' )
304  safmax = one/safmin
305  CALL dlabad( safmin, safmax )
306  ulp = dlamch( 'PRECISION' )
307  smlnum = safmin*( dble( n )/ulp )
308 
309  IF ( ihi .EQ. kwtop ) THEN
310 * 1 by 1 deflation window, just try a regular deflation
311  alphar( kwtop ) = a( kwtop, kwtop )
312  alphai( kwtop ) = zero
313  beta( kwtop ) = b( kwtop, kwtop )
314  ns = 1
315  nd = 0
316  IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
317  $ kwtop ) ) ) ) THEN
318  ns = 0
319  nd = 1
320  IF ( kwtop .GT. ilo ) THEN
321  a( kwtop, kwtop-1 ) = zero
322  END IF
323  END IF
324  END IF
325 
326 
327 * Store window in case of convergence failure
328  CALL dlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
329  CALL dlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
330  $ 1 ), jw )
331 
332 * Transform window to real schur form
333  CALL dlaset( 'FULL', jw, jw, zero, one, qc, ldqc )
334  CALL dlaset( 'FULL', jw, jw, zero, one, zc, ldzc )
335  CALL dlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
336  $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
337  $ ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
338  $ rec+1, qz_small_info )
339 
340  IF( qz_small_info .NE. 0 ) THEN
341 * Convergence failure, restore the window and exit
342  nd = 0
343  ns = jw-qz_small_info
344  CALL dlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
345  CALL dlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
346  $ kwtop ), ldb )
347  RETURN
348  END IF
349 
350 * Deflation detection loop
351  IF ( kwtop .EQ. ilo .OR. s .EQ. zero ) THEN
352  kwbot = kwtop-1
353  ELSE
354  kwbot = ihi
355  k = 1
356  k2 = 1
357  DO WHILE ( k .LE. jw )
358  bulge = .false.
359  IF ( kwbot-kwtop+1 .GE. 2 ) THEN
360  bulge = a( kwbot, kwbot-1 ) .NE. zero
361  END IF
362  IF ( bulge ) THEN
363 
364 * Try to deflate complex conjugate eigenvalue pair
365  temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
366  $ kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
367  IF( temp .EQ. zero )THEN
368  temp = abs( s )
369  END IF
370  IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
371  $ kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
372  $ ulp*temp ) ) THEN
373 * Deflatable
374  kwbot = kwbot-2
375  ELSE
376 * Not deflatable, move out of the way
377  ifst = kwbot-kwtop+1
378  ilst = k2
379  CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
380  $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
381  $ zc, ldzc, ifst, ilst, work, lwork,
382  $ dtgexc_info )
383  k2 = k2+2
384  END IF
385  k = k+2
386  ELSE
387 
388 * Try to deflate real eigenvalue
389  temp = abs( a( kwbot, kwbot ) )
390  IF( temp .EQ. zero ) THEN
391  temp = abs( s )
392  END IF
393  IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
394  $ temp, smlnum ) ) THEN
395 * Deflatable
396  kwbot = kwbot-1
397  ELSE
398 * Not deflatable, move out of the way
399  ifst = kwbot-kwtop+1
400  ilst = k2
401  CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
402  $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
403  $ zc, ldzc, ifst, ilst, work, lwork,
404  $ dtgexc_info )
405  k2 = k2+1
406  END IF
407 
408  k = k+1
409 
410  END IF
411  END DO
412  END IF
413 
414 * Store eigenvalues
415  nd = ihi-kwbot
416  ns = jw-nd
417  k = kwtop
418  DO WHILE ( k .LE. ihi )
419  bulge = .false.
420  IF ( k .LT. ihi ) THEN
421  IF ( a( k+1, k ) .NE. zero ) THEN
422  bulge = .true.
423  END IF
424  END IF
425  IF ( bulge ) THEN
426 * 2x2 eigenvalue block
427  CALL dlag2( a( k, k ), lda, b( k, k ), ldb, safmin,
428  $ beta( k ), beta( k+1 ), alphar( k ),
429  $ alphar( k+1 ), alphai( k ) )
430  alphai( k+1 ) = -alphai( k )
431  k = k+2
432  ELSE
433 * 1x1 eigenvalue block
434  alphar( k ) = a( k, k )
435  alphai( k ) = zero
436  beta( k ) = b( k, k )
437  k = k+1
438  END IF
439  END DO
440 
441  IF ( kwtop .NE. ilo .AND. s .NE. zero ) THEN
442 * Reflect spike back, this will create optimally packed bulges
443  a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
444  $ 1:jw-nd )
445  DO k = kwbot-1, kwtop, -1
446  CALL dlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
447  $ temp )
448  a( k, kwtop-1 ) = temp
449  a( k+1, kwtop-1 ) = zero
450  k2 = max( kwtop, k-1 )
451  CALL drot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
452  $ s1 )
453  CALL drot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
454  $ ldb, c1, s1 )
455  CALL drot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
456  $ 1, c1, s1 )
457  END DO
458 
459 * Chase bulges down
460  istartm = kwtop
461  istopm = ihi
462  k = kwbot-1
463  DO WHILE ( k .GE. kwtop )
464  IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero ) THEN
465 
466 * Move double pole block down and remove it
467  DO k2 = k-1, kwbot-2
468  CALL dlaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
469  $ kwbot, a, lda, b, ldb, jw, kwtop, qc,
470  $ ldqc, jw, kwtop, zc, ldzc )
471  END DO
472 
473  k = k-2
474  ELSE
475 
476 * k points to single shift
477  DO k2 = k, kwbot-2
478 
479 * Move shift down
480  CALL dlartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
481  $ temp )
482  b( k2+1, k2+1 ) = temp
483  b( k2+1, k2 ) = zero
484  CALL drot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
485  $ a( istartm, k2 ), 1, c1, s1 )
486  CALL drot( k2-istartm+1, b( istartm, k2+1 ), 1,
487  $ b( istartm, k2 ), 1, c1, s1 )
488  CALL drot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
489  $ k2-kwtop+1 ), 1, c1, s1 )
490 
491  CALL dlartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
492  $ temp )
493  a( k2+1, k2 ) = temp
494  a( k2+2, k2 ) = zero
495  CALL drot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
496  $ k2+1 ), lda, c1, s1 )
497  CALL drot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
498  $ k2+1 ), ldb, c1, s1 )
499  CALL drot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
500  $ k2+2-kwtop+1 ), 1, c1, s1 )
501 
502  END DO
503 
504 * Remove the shift
505  CALL dlartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
506  $ s1, temp )
507  b( kwbot, kwbot ) = temp
508  b( kwbot, kwbot-1 ) = zero
509  CALL drot( kwbot-istartm, b( istartm, kwbot ), 1,
510  $ b( istartm, kwbot-1 ), 1, c1, s1 )
511  CALL drot( kwbot-istartm+1, a( istartm, kwbot ), 1,
512  $ a( istartm, kwbot-1 ), 1, c1, s1 )
513  CALL drot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
514  $ kwbot-1-kwtop+1 ), 1, c1, s1 )
515 
516  k = k-1
517  END IF
518  END DO
519 
520  END IF
521 
522 * Apply Qc and Zc to rest of the matrix
523  IF ( ilschur ) THEN
524  istartm = 1
525  istopm = n
526  ELSE
527  istartm = ilo
528  istopm = ihi
529  END IF
530 
531  IF ( istopm-ihi > 0 ) THEN
532  CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
533  $ a( kwtop, ihi+1 ), lda, zero, work, jw )
534  CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
535  $ ihi+1 ), lda )
536  CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
537  $ b( kwtop, ihi+1 ), ldb, zero, work, jw )
538  CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
539  $ ihi+1 ), ldb )
540  END IF
541  IF ( ilq ) THEN
542  CALL dgemm( 'N', 'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
543  $ ldqc, zero, work, n )
544  CALL dlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
545  END IF
546 
547  IF ( kwtop-1-istartm+1 > 0 ) THEN
548  CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, a( istartm,
549  $ kwtop ), lda, zc, ldzc, zero, work,
550  $ kwtop-istartm )
551  CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
552  $ a( istartm, kwtop ), lda )
553  CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, b( istartm,
554  $ kwtop ), ldb, zc, ldzc, zero, work,
555  $ kwtop-istartm )
556  CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
557  $ b( istartm, kwtop ), ldb )
558  END IF
559  IF ( ilz ) THEN
560  CALL dgemm( 'N', 'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
561  $ ldzc, zero, work, n )
562  CALL dlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
563  END IF
564 
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dtgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
DTGEXC
Definition: dtgexc.f:220
recursive subroutine dlaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, REC, INFO)
DLAQZ0
Definition: dlaqz0.f:306
subroutine dlaqz2(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
DLAQZ2
Definition: dlaqz2.f:174
subroutine dlag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: dlag2.f:156
Here is the call graph for this function:
Here is the caller graph for this function: