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

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.```
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 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: