LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ dlaqz0()

recursive subroutine dlaqz0 ( character, intent(in)  WANTS,
character, intent(in)  WANTQ,
character, intent(in)  WANTZ,
integer, intent(in)  N,
integer, intent(in)  ILO,
integer, intent(in)  IHI,
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( * ), intent(inout)  ALPHAR,
double precision, dimension( * ), intent(inout)  ALPHAI,
double precision, dimension( * ), intent(inout)  BETA,
double precision, dimension( ldq, * ), intent(inout)  Q,
integer, intent(in)  LDQ,
double precision, dimension( ldz, * ), intent(inout)  Z,
integer, intent(in)  LDZ,
double precision, dimension( * ), intent(inout)  WORK,
integer, intent(in)  LWORK,
integer, intent(in)  REC,
integer, intent(out)  INFO 
)

DLAQZ0

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

Purpose:
 DLAQZ0 computes the eigenvalues of a real matrix pair (H,T),
 where H is an upper Hessenberg matrix and T is upper triangular,
 using the double-shift QZ method.
 Matrix pairs of this type are produced by the reduction to
 generalized upper Hessenberg form of a real matrix pair (A,B):

    A = Q1*H*Z1**T,  B = Q1*T*Z1**T,

 as computed by DGGHRD.

 If JOB='S', then the Hessenberg-triangular pair (H,T) is
 also reduced to generalized Schur form,

    H = Q*S*Z**T,  T = Q*P*Z**T,

 where Q and Z are orthogonal matrices, P is an upper triangular
 matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
 diagonal blocks.

 The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
 (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
 eigenvalues.

 Additionally, the 2-by-2 upper triangular diagonal blocks of P
 corresponding to 2-by-2 blocks of S are reduced to positive diagonal
 form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
 P(j,j) > 0, and P(j+1,j+1) > 0.

 Optionally, the orthogonal matrix Q from the generalized Schur
 factorization may be postmultiplied into an input matrix Q1, and the
 orthogonal matrix Z may be postmultiplied into an input matrix Z1.
 If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
 the matrix pair (A,B) to generalized upper Hessenberg form, then the
 output matrices Q1*Q and Z1*Z are the orthogonal factors from the
 generalized Schur factorization of (A,B):

    A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.

 To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
 of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
 complex and beta real.
 If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
 generalized nonsymmetric eigenvalue problem (GNEP)
    A*x = lambda*B*x
 and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
 alternate form of the GNEP
    mu*A*y = B*y.
 Real eigenvalues can be read directly from the generalized Schur
 form:
   alpha = S(i,i), beta = P(i,i).

 Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
      Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
      pp. 241--256.

 Ref: B. Kagstrom, D. Kressner, "Multishift Variants of the QZ
      Algorithm with Aggressive Early Deflation", SIAM J. Numer.
      Anal., 29(2006), pp. 199--227.

 Ref: T. Steel, D. Camps, K. Meerbergen, R. Vandebril "A multishift,
      multipole rational QZ method with agressive early deflation"
Parameters
[in]WANTS
          WANTS is CHARACTER*1
          = 'E': Compute eigenvalues only;
          = 'S': Compute eigenvalues and the Schur form.
[in]WANTQ
          WANTQ is CHARACTER*1
          = 'N': Left Schur vectors (Q) are not computed;
          = 'I': Q is initialized to the unit matrix and the matrix Q
                 of left Schur vectors of (A,B) is returned;
          = 'V': Q must contain an orthogonal matrix Q1 on entry and
                 the product Q1*Q is returned.
[in]WANTZ
          WANTZ is CHARACTER*1
          = 'N': Right Schur vectors (Z) are not computed;
          = 'I': Z is initialized to the unit matrix and the matrix Z
                 of right Schur vectors of (A,B) is returned;
          = 'V': Z must contain an orthogonal matrix Z1 on entry and
                 the product Z1*Z is returned.
[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 which are in
          Hessenberg form.  It is assumed that A is already upper
          triangular in rows and columns 1:ILO-1 and IHI+1:N.
          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA, N)
          On entry, the N-by-N upper Hessenberg matrix A.
          On exit, if JOB = 'S', A contains the upper quasi-triangular
          matrix S from the generalized Schur factorization.
          If JOB = 'E', the diagonal blocks of A match those of S, but
          the rest of A is unspecified.
[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)
          On entry, the N-by-N upper triangular matrix B.
          On exit, if JOB = 'S', B contains the upper triangular
          matrix P from the generalized Schur factorization;
          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
          are reduced to positive diagonal form, i.e., if A(j+1,j) is
          non-zero, then B(j+1,j) = B(j,j+1) = 0, B(j,j) > 0, and
          B(j+1,j+1) > 0.
          If JOB = 'E', the diagonal blocks of B match those of P, but
          the rest of B is unspecified.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B.  LDB >= max( 1, N ).
[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]Q
          Q is DOUBLE PRECISION array, dimension (LDQ, N)
          On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
          the reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
          vectors of (A,B), and if COMPQ = 'V', the orthogonal matrix
          of left Schur vectors of (A,B).
          Not referenced if COMPQ = 'N'.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  LDQ >= 1.
          If COMPQ='V' or 'I', then LDQ >= N.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ, N)
          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
          the reduction of (A,B) to generalized Hessenberg form.
          On exit, if COMPZ = 'I', the orthogonal matrix of
          right Schur vectors of (H,T), and if COMPZ = 'V', the
          orthogonal matrix of right Schur vectors of (A,B).
          Not referenced if COMPZ = 'N'.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1.
          If COMPZ='V' or 'I', then LDZ >= N.
[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
          = 1,...,N: the QZ iteration did not converge.  (A,B) is not
                     in Schur form, but ALPHAR(i), ALPHAI(i), and
                     BETA(i), i=INFO+1,...,N should be correct.
Author
Thijs Steel, KU Leuven
Date
May 2020

Definition at line 302 of file dlaqz0.f.

306  IMPLICIT NONE
307 
308 * Arguments
309  CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ
310  INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
311  $ REC
312 
313  INTEGER, INTENT( OUT ) :: INFO
314 
315  DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
316  $ Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ),
317  $ ALPHAI( * ), BETA( * ), WORK( * )
318 
319 * Parameters
320  DOUBLE PRECISION :: ZERO, ONE, HALF
321  parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
322 
323 * Local scalars
324  DOUBLE PRECISION :: SMLNUM, ULP, ESHIFT, SAFMIN, SAFMAX, C1, S1,
325  $ TEMP, SWAP
326  INTEGER :: ISTART, ISTOP, IITER, MAXIT, ISTART2, K, LD, NSHIFTS,
327  $ NBLOCK, NW, NMIN, NIBBLE, N_UNDEFLATED, N_DEFLATED,
328  $ NS, SWEEP_INFO, SHIFTPOS, LWORKREQ, K2, ISTARTM,
329  $ ISTOPM, IWANTS, IWANTQ, IWANTZ, NORM_INFO, AED_INFO,
330  $ NWR, NBR, NSR, ITEMP1, ITEMP2, RCOST, I
331  LOGICAL :: ILSCHUR, ILQ, ILZ
332  CHARACTER :: JBCMPZ*3
333 
334 * External Functions
335  EXTERNAL :: xerbla, dhgeqz, dlaset, dlaqz3, dlaqz4, dlabad,
336  $ dlartg, drot
337  DOUBLE PRECISION, EXTERNAL :: DLAMCH
338  LOGICAL, EXTERNAL :: LSAME
339  INTEGER, EXTERNAL :: ILAENV
340 
341 *
342 * Decode wantS,wantQ,wantZ
343 *
344  IF( lsame( wants, 'E' ) ) THEN
345  ilschur = .false.
346  iwants = 1
347  ELSE IF( lsame( wants, 'S' ) ) THEN
348  ilschur = .true.
349  iwants = 2
350  ELSE
351  iwants = 0
352  END IF
353 
354  IF( lsame( wantq, 'N' ) ) THEN
355  ilq = .false.
356  iwantq = 1
357  ELSE IF( lsame( wantq, 'V' ) ) THEN
358  ilq = .true.
359  iwantq = 2
360  ELSE IF( lsame( wantq, 'I' ) ) THEN
361  ilq = .true.
362  iwantq = 3
363  ELSE
364  iwantq = 0
365  END IF
366 
367  IF( lsame( wantz, 'N' ) ) THEN
368  ilz = .false.
369  iwantz = 1
370  ELSE IF( lsame( wantz, 'V' ) ) THEN
371  ilz = .true.
372  iwantz = 2
373  ELSE IF( lsame( wantz, 'I' ) ) THEN
374  ilz = .true.
375  iwantz = 3
376  ELSE
377  iwantz = 0
378  END IF
379 *
380 * Check Argument Values
381 *
382  info = 0
383  IF( iwants.EQ.0 ) THEN
384  info = -1
385  ELSE IF( iwantq.EQ.0 ) THEN
386  info = -2
387  ELSE IF( iwantz.EQ.0 ) THEN
388  info = -3
389  ELSE IF( n.LT.0 ) THEN
390  info = -4
391  ELSE IF( ilo.LT.1 ) THEN
392  info = -5
393  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
394  info = -6
395  ELSE IF( lda.LT.n ) THEN
396  info = -8
397  ELSE IF( ldb.LT.n ) THEN
398  info = -10
399  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
400  info = -15
401  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
402  info = -17
403  END IF
404  IF( info.NE.0 ) THEN
405  CALL xerbla( 'DLAQZ0', -info )
406  RETURN
407  END IF
408 
409 *
410 * Quick return if possible
411 *
412  IF( n.LE.0 ) THEN
413  work( 1 ) = dble( 1 )
414  RETURN
415  END IF
416 
417 *
418 * Get the parameters
419 *
420  jbcmpz( 1:1 ) = wants
421  jbcmpz( 2:2 ) = wantq
422  jbcmpz( 3:3 ) = wantz
423 
424  nmin = ilaenv( 12, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
425 
426  nwr = ilaenv( 13, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
427  nwr = max( 2, nwr )
428  nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
429 
430  nibble = ilaenv( 14, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
431 
432  nsr = ilaenv( 15, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
433  nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
434  nsr = max( 2, nsr-mod( nsr, 2 ) )
435 
436  rcost = ilaenv( 17, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
437  itemp1 = int( nsr/sqrt( 1+2*nsr/( dble( rcost )/100*n ) ) )
438  itemp1 = ( ( itemp1-1 )/4 )*4+4
439  nbr = nsr+itemp1
440 
441  IF( n .LT. nmin .OR. rec .GE. 2 ) THEN
442  CALL dhgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
443  $ alphar, alphai, beta, q, ldq, z, ldz, work,
444  $ lwork, info )
445  RETURN
446  END IF
447 
448 *
449 * Find out required workspace
450 *
451 
452 * Workspace query to dlaqz3
453  nw = max( nwr, nmin )
454  CALL dlaqz3( ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb,
455  $ q, ldq, z, ldz, n_undeflated, n_deflated, alphar,
456  $ alphai, beta, work, nw, work, nw, work, -1, rec,
457  $ aed_info )
458  itemp1 = int( work( 1 ) )
459 * Workspace query to dlaqz4
460  CALL dlaqz4( ilschur, ilq, ilz, n, ilo, ihi, nsr, nbr, alphar,
461  $ alphai, beta, a, lda, b, ldb, q, ldq, z, ldz, work,
462  $ nbr, work, nbr, work, -1, sweep_info )
463  itemp2 = int( work( 1 ) )
464 
465  lworkreq = max( itemp1+2*nw**2, itemp2+2*nbr**2 )
466  IF ( lwork .EQ.-1 ) THEN
467  work( 1 ) = dble( lworkreq )
468  RETURN
469  ELSE IF ( lwork .LT. lworkreq ) THEN
470  info = -19
471  END IF
472  IF( info.NE.0 ) THEN
473  CALL xerbla( 'DLAQZ0', info )
474  RETURN
475  END IF
476 *
477 * Initialize Q and Z
478 *
479  IF( iwantq.EQ.3 ) CALL dlaset( 'FULL', n, n, zero, one, q, ldq )
480  IF( iwantz.EQ.3 ) CALL dlaset( 'FULL', n, n, zero, one, z, ldz )
481 
482 * Get machine constants
483  safmin = dlamch( 'SAFE MINIMUM' )
484  safmax = one/safmin
485  CALL dlabad( safmin, safmax )
486  ulp = dlamch( 'PRECISION' )
487  smlnum = safmin*( dble( n )/ulp )
488 
489  istart = ilo
490  istop = ihi
491  maxit = 3*( ihi-ilo+1 )
492  ld = 0
493 
494  DO iiter = 1, maxit
495  IF( iiter .GE. maxit ) THEN
496  info = istop+1
497  GOTO 80
498  END IF
499  IF ( istart+1 .GE. istop ) THEN
500  istop = istart
501  EXIT
502  END IF
503 
504 * Check deflations at the end
505  IF ( abs( a( istop-1, istop-2 ) ) .LE. max( smlnum,
506  $ ulp*( abs( a( istop-1, istop-1 ) )+abs( a( istop-2,
507  $ istop-2 ) ) ) ) ) THEN
508  a( istop-1, istop-2 ) = zero
509  istop = istop-2
510  ld = 0
511  eshift = zero
512  ELSE IF ( abs( a( istop, istop-1 ) ) .LE. max( smlnum,
513  $ ulp*( abs( a( istop, istop ) )+abs( a( istop-1,
514  $ istop-1 ) ) ) ) ) THEN
515  a( istop, istop-1 ) = zero
516  istop = istop-1
517  ld = 0
518  eshift = zero
519  END IF
520 * Check deflations at the start
521  IF ( abs( a( istart+2, istart+1 ) ) .LE. max( smlnum,
522  $ ulp*( abs( a( istart+1, istart+1 ) )+abs( a( istart+2,
523  $ istart+2 ) ) ) ) ) THEN
524  a( istart+2, istart+1 ) = zero
525  istart = istart+2
526  ld = 0
527  eshift = zero
528  ELSE IF ( abs( a( istart+1, istart ) ) .LE. max( smlnum,
529  $ ulp*( abs( a( istart, istart ) )+abs( a( istart+1,
530  $ istart+1 ) ) ) ) ) THEN
531  a( istart+1, istart ) = zero
532  istart = istart+1
533  ld = 0
534  eshift = zero
535  END IF
536 
537  IF ( istart+1 .GE. istop ) THEN
538  EXIT
539  END IF
540 
541 * Check interior deflations
542  istart2 = istart
543  DO k = istop, istart+1, -1
544  IF ( abs( a( k, k-1 ) ) .LE. max( smlnum, ulp*( abs( a( k,
545  $ k ) )+abs( a( k-1, k-1 ) ) ) ) ) THEN
546  a( k, k-1 ) = zero
547  istart2 = k
548  EXIT
549  END IF
550  END DO
551 
552 * Get range to apply rotations to
553  IF ( ilschur ) THEN
554  istartm = 1
555  istopm = n
556  ELSE
557  istartm = istart2
558  istopm = istop
559  END IF
560 
561 * Check infinite eigenvalues, this is done without blocking so might
562 * slow down the method when many infinite eigenvalues are present
563  k = istop
564  DO WHILE ( k.GE.istart2 )
565  temp = zero
566  IF( k .LT. istop ) THEN
567  temp = temp+abs( b( k, k+1 ) )
568  END IF
569  IF( k .GT. istart2 ) THEN
570  temp = temp+abs( b( k-1, k ) )
571  END IF
572 
573  IF( abs( b( k, k ) ) .LT. max( smlnum, ulp*temp ) ) THEN
574 * A diagonal element of B is negligable, move it
575 * to the top and deflate it
576 
577  DO k2 = k, istart2+1, -1
578  CALL dlartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1, s1,
579  $ temp )
580  b( k2-1, k2 ) = temp
581  b( k2-1, k2-1 ) = zero
582 
583  CALL drot( k2-2-istartm+1, b( istartm, k2 ), 1,
584  $ b( istartm, k2-1 ), 1, c1, s1 )
585  CALL drot( min( k2+1, istop )-istartm+1, a( istartm,
586  $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
587  IF ( ilz ) THEN
588  CALL drot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1, c1,
589  $ s1 )
590  END IF
591 
592  IF( k2.LT.istop ) THEN
593  CALL dlartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
594  $ s1, temp )
595  a( k2, k2-1 ) = temp
596  a( k2+1, k2-1 ) = zero
597 
598  CALL drot( istopm-k2+1, a( k2, k2 ), lda, a( k2+1,
599  $ k2 ), lda, c1, s1 )
600  CALL drot( istopm-k2+1, b( k2, k2 ), ldb, b( k2+1,
601  $ k2 ), ldb, c1, s1 )
602  IF( ilq ) THEN
603  CALL drot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
604  $ c1, s1 )
605  END IF
606  END IF
607 
608  END DO
609 
610  IF( istart2.LT.istop )THEN
611  CALL dlartg( a( istart2, istart2 ), a( istart2+1,
612  $ istart2 ), c1, s1, temp )
613  a( istart2, istart2 ) = temp
614  a( istart2+1, istart2 ) = zero
615 
616  CALL drot( istopm-( istart2+1 )+1, a( istart2,
617  $ istart2+1 ), lda, a( istart2+1,
618  $ istart2+1 ), lda, c1, s1 )
619  CALL drot( istopm-( istart2+1 )+1, b( istart2,
620  $ istart2+1 ), ldb, b( istart2+1,
621  $ istart2+1 ), ldb, c1, s1 )
622  IF( ilq ) THEN
623  CALL drot( n, q( 1, istart2 ), 1, q( 1,
624  $ istart2+1 ), 1, c1, s1 )
625  END IF
626  END IF
627 
628  istart2 = istart2+1
629 
630  END IF
631  k = k-1
632  END DO
633 
634 * istart2 now points to the top of the bottom right
635 * unreduced Hessenberg block
636  IF ( istart2 .GE. istop ) THEN
637  istop = istart2-1
638  ld = 0
639  eshift = zero
640  cycle
641  END IF
642 
643  nw = nwr
644  nshifts = nsr
645  nblock = nbr
646 
647  IF ( istop-istart2+1 .LT. nmin ) THEN
648 * Setting nw to the size of the subblock will make AED deflate
649 * all the eigenvalues. This is slightly more efficient than just
650 * using DHGEQZ because the off diagonal part gets updated via BLAS.
651  IF ( istop-istart+1 .LT. nmin ) THEN
652  nw = istop-istart+1
653  istart2 = istart
654  ELSE
655  nw = istop-istart2+1
656  END IF
657  END IF
658 
659 *
660 * Time for AED
661 *
662  CALL dlaqz3( ilschur, ilq, ilz, n, istart2, istop, nw, a, lda,
663  $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
664  $ alphar, alphai, beta, work, nw, work( nw**2+1 ),
665  $ nw, work( 2*nw**2+1 ), lwork-2*nw**2, rec,
666  $ aed_info )
667 
668  IF ( n_deflated > 0 ) THEN
669  istop = istop-n_deflated
670  ld = 0
671  eshift = zero
672  END IF
673 
674  IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
675  $ istop-istart2+1 .LT. nmin ) THEN
676 * AED has uncovered many eigenvalues. Skip a QZ sweep and run
677 * AED again.
678  cycle
679  END IF
680 
681  ld = ld+1
682 
683  ns = min( nshifts, istop-istart2 )
684  ns = min( ns, n_undeflated )
685  shiftpos = istop-n_deflated-n_undeflated+1
686 *
687 * Shuffle shifts to put double shifts in front
688 * This ensures that we don't split up a double shift
689 *
690  DO i = shiftpos, shiftpos+n_undeflated-1, 2
691  IF( alphai( i ).NE.-alphai( i+1 ) ) THEN
692 *
693  swap = alphar( i )
694  alphar( i ) = alphar( i+1 )
695  alphar( i+1 ) = alphar( i+2 )
696  alphar( i+2 ) = swap
697 
698  swap = alphai( i )
699  alphai( i ) = alphai( i+1 )
700  alphai( i+1 ) = alphai( i+2 )
701  alphai( i+2 ) = swap
702 
703  swap = beta( i )
704  beta( i ) = beta( i+1 )
705  beta( i+1 ) = beta( i+2 )
706  beta( i+2 ) = swap
707  END IF
708  END DO
709 
710  IF ( mod( ld, 6 ) .EQ. 0 ) THEN
711 *
712 * Exceptional shift. Chosen for no particularly good reason.
713 *
714  IF( ( dble( maxit )*safmin )*abs( a( istop,
715  $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) ) THEN
716  eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
717  ELSE
718  eshift = eshift+one/( safmin*dble( maxit ) )
719  END IF
720  alphar( shiftpos ) = one
721  alphar( shiftpos+1 ) = zero
722  alphai( shiftpos ) = zero
723  alphai( shiftpos+1 ) = zero
724  beta( shiftpos ) = eshift
725  beta( shiftpos+1 ) = eshift
726  ns = 2
727  END IF
728 
729 *
730 * Time for a QZ sweep
731 *
732  CALL dlaqz4( ilschur, ilq, ilz, n, istart2, istop, ns, nblock,
733  $ alphar( shiftpos ), alphai( shiftpos ),
734  $ beta( shiftpos ), a, lda, b, ldb, q, ldq, z, ldz,
735  $ work, nblock, work( nblock**2+1 ), nblock,
736  $ work( 2*nblock**2+1 ), lwork-2*nblock**2,
737  $ sweep_info )
738 
739  END DO
740 
741 *
742 * Call DHGEQZ to normalize the eigenvalue blocks and set the eigenvalues
743 * If all the eigenvalues have been found, DHGEQZ will not do any iterations
744 * and only normalize the blocks. In case of a rare convergence failure,
745 * the single shift might perform better.
746 *
747  80 CALL dhgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
748  $ alphar, alphai, beta, q, ldq, z, ldz, work, lwork,
749  $ norm_info )
750 
751  info = norm_info
752 
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
Definition: dhgeqz.f:304
subroutine dlaqz4(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS, NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK, INFO)
DLAQZ4
Definition: dlaqz4.f:213
recursive subroutine dlaqz3(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHAR, ALPHAI, BETA, QC, LDQC, ZC, LDZC, WORK, LWORK, REC, INFO)
DLAQZ3
Definition: dlaqz3.f:239
Here is the call graph for this function:
Here is the caller graph for this function: