LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dhgeqz ( character  JOB,
character  COMPQ,
character  COMPZ,
integer  N,
integer  ILO,
integer  IHI,
double precision, dimension( ldh, * )  H,
integer  LDH,
double precision, dimension( ldt, * )  T,
integer  LDT,
double precision, dimension( * )  ALPHAR,
double precision, dimension( * )  ALPHAI,
double precision, dimension( * )  BETA,
double precision, dimension( ldq, * )  Q,
integer  LDQ,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DHGEQZ

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

Purpose:
 DHGEQZ 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.
Parameters
[in]JOB
          JOB is CHARACTER*1
          = 'E': Compute eigenvalues only;
          = 'S': Compute eigenvalues and the Schur form. 
[in]COMPQ
          COMPQ 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 (H,T) is returned;
          = 'V': Q must contain an orthogonal matrix Q1 on entry and
                 the product Q1*Q is returned.
[in]COMPZ
          COMPZ 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 (H,T) 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 H, T, Q, and Z.  N >= 0.
[in]ILO
          ILO is INTEGER
[in]IHI
          IHI is INTEGER
          ILO and IHI mark the rows and columns of H 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]H
          H is DOUBLE PRECISION array, dimension (LDH, N)
          On entry, the N-by-N upper Hessenberg matrix H.
          On exit, if JOB = 'S', H contains the upper quasi-triangular
          matrix S from the generalized Schur factorization.
          If JOB = 'E', the diagonal blocks of H match those of S, but
          the rest of H is unspecified.
[in]LDH
          LDH is INTEGER
          The leading dimension of the array H.  LDH >= max( 1, N ).
[in,out]T
          T is DOUBLE PRECISION array, dimension (LDT, N)
          On entry, the N-by-N upper triangular matrix T.
          On exit, if JOB = 'S', T 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 H(j+1,j) is
          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
          T(j+1,j+1) > 0.
          If JOB = 'E', the diagonal blocks of T match those of P, but
          the rest of T is unspecified.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T.  LDT >= 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 (H,T), 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.
[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.  (H,T) is not
                     in Schur form, but ALPHAR(i), ALPHAI(i), and
                     BETA(i), i=INFO+1,...,N should be correct.
          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
                     in Schur form, but ALPHAR(i), ALPHAI(i), and
                     BETA(i), i=INFO-N+1,...,N should be correct.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016
Further Details:
  Iteration counters:

  JITER  -- counts iterations.
  IITER  -- counts iterations run since ILAST was last
            changed.  This is therefore reset only when a 1-by-1 or
            2-by-2 block deflates off the bottom.

Definition at line 306 of file dhgeqz.f.

306 *
307 * -- LAPACK computational routine (version 3.6.1) --
308 * -- LAPACK is a software package provided by Univ. of Tennessee, --
309 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310 * June 2016
311 *
312 * .. Scalar Arguments ..
313  CHARACTER compq, compz, job
314  INTEGER ihi, ilo, info, ldh, ldq, ldt, ldz, lwork, n
315 * ..
316 * .. Array Arguments ..
317  DOUBLE PRECISION alphai( * ), alphar( * ), beta( * ),
318  $ h( ldh, * ), q( ldq, * ), t( ldt, * ),
319  $ work( * ), z( ldz, * )
320 * ..
321 *
322 * =====================================================================
323 *
324 * .. Parameters ..
325 * $ SAFETY = 1.0E+0 )
326  DOUBLE PRECISION half, zero, one, safety
327  parameter ( half = 0.5d+0, zero = 0.0d+0, one = 1.0d+0,
328  $ safety = 1.0d+2 )
329 * ..
330 * .. Local Scalars ..
331  LOGICAL ilazr2, ilazro, ilpivt, ilq, ilschr, ilz,
332  $ lquery
333  INTEGER icompq, icompz, ifirst, ifrstm, iiter, ilast,
334  $ ilastm, in, ischur, istart, j, jc, jch, jiter,
335  $ jr, maxit
336  DOUBLE PRECISION a11, a12, a1i, a1r, a21, a22, a2i, a2r, ad11,
337  $ ad11l, ad12, ad12l, ad21, ad21l, ad22, ad22l,
338  $ ad32l, an, anorm, ascale, atol, b11, b1a, b1i,
339  $ b1r, b22, b2a, b2i, b2r, bn, bnorm, bscale,
340  $ btol, c, c11i, c11r, c12, c21, c22i, c22r, cl,
341  $ cq, cr, cz, eshift, s, s1, s1inv, s2, safmax,
342  $ safmin, scale, sl, sqi, sqr, sr, szi, szr, t1,
343  $ tau, temp, temp2, tempi, tempr, u1, u12, u12l,
344  $ u2, ulp, vs, w11, w12, w21, w22, wabs, wi, wr,
345  $ wr2
346 * ..
347 * .. Local Arrays ..
348  DOUBLE PRECISION v( 3 )
349 * ..
350 * .. External Functions ..
351  LOGICAL lsame
352  DOUBLE PRECISION dlamch, dlanhs, dlapy2, dlapy3
353  EXTERNAL lsame, dlamch, dlanhs, dlapy2, dlapy3
354 * ..
355 * .. External Subroutines ..
356  EXTERNAL dlag2, dlarfg, dlartg, dlaset, dlasv2, drot,
357  $ xerbla
358 * ..
359 * .. Intrinsic Functions ..
360  INTRINSIC abs, dble, max, min, sqrt
361 * ..
362 * .. Executable Statements ..
363 *
364 * Decode JOB, COMPQ, COMPZ
365 *
366  IF( lsame( job, 'E' ) ) THEN
367  ilschr = .false.
368  ischur = 1
369  ELSE IF( lsame( job, 'S' ) ) THEN
370  ilschr = .true.
371  ischur = 2
372  ELSE
373  ischur = 0
374  END IF
375 *
376  IF( lsame( compq, 'N' ) ) THEN
377  ilq = .false.
378  icompq = 1
379  ELSE IF( lsame( compq, 'V' ) ) THEN
380  ilq = .true.
381  icompq = 2
382  ELSE IF( lsame( compq, 'I' ) ) THEN
383  ilq = .true.
384  icompq = 3
385  ELSE
386  icompq = 0
387  END IF
388 *
389  IF( lsame( compz, 'N' ) ) THEN
390  ilz = .false.
391  icompz = 1
392  ELSE IF( lsame( compz, 'V' ) ) THEN
393  ilz = .true.
394  icompz = 2
395  ELSE IF( lsame( compz, 'I' ) ) THEN
396  ilz = .true.
397  icompz = 3
398  ELSE
399  icompz = 0
400  END IF
401 *
402 * Check Argument Values
403 *
404  info = 0
405  work( 1 ) = max( 1, n )
406  lquery = ( lwork.EQ.-1 )
407  IF( ischur.EQ.0 ) THEN
408  info = -1
409  ELSE IF( icompq.EQ.0 ) THEN
410  info = -2
411  ELSE IF( icompz.EQ.0 ) THEN
412  info = -3
413  ELSE IF( n.LT.0 ) THEN
414  info = -4
415  ELSE IF( ilo.LT.1 ) THEN
416  info = -5
417  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
418  info = -6
419  ELSE IF( ldh.LT.n ) THEN
420  info = -8
421  ELSE IF( ldt.LT.n ) THEN
422  info = -10
423  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
424  info = -15
425  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
426  info = -17
427  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
428  info = -19
429  END IF
430  IF( info.NE.0 ) THEN
431  CALL xerbla( 'DHGEQZ', -info )
432  RETURN
433  ELSE IF( lquery ) THEN
434  RETURN
435  END IF
436 *
437 * Quick return if possible
438 *
439  IF( n.LE.0 ) THEN
440  work( 1 ) = dble( 1 )
441  RETURN
442  END IF
443 *
444 * Initialize Q and Z
445 *
446  IF( icompq.EQ.3 )
447  $ CALL dlaset( 'Full', n, n, zero, one, q, ldq )
448  IF( icompz.EQ.3 )
449  $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
450 *
451 * Machine Constants
452 *
453  in = ihi + 1 - ilo
454  safmin = dlamch( 'S' )
455  safmax = one / safmin
456  ulp = dlamch( 'E' )*dlamch( 'B' )
457  anorm = dlanhs( 'F', in, h( ilo, ilo ), ldh, work )
458  bnorm = dlanhs( 'F', in, t( ilo, ilo ), ldt, work )
459  atol = max( safmin, ulp*anorm )
460  btol = max( safmin, ulp*bnorm )
461  ascale = one / max( safmin, anorm )
462  bscale = one / max( safmin, bnorm )
463 *
464 * Set Eigenvalues IHI+1:N
465 *
466  DO 30 j = ihi + 1, n
467  IF( t( j, j ).LT.zero ) THEN
468  IF( ilschr ) THEN
469  DO 10 jr = 1, j
470  h( jr, j ) = -h( jr, j )
471  t( jr, j ) = -t( jr, j )
472  10 CONTINUE
473  ELSE
474  h( j, j ) = -h( j, j )
475  t( j, j ) = -t( j, j )
476  END IF
477  IF( ilz ) THEN
478  DO 20 jr = 1, n
479  z( jr, j ) = -z( jr, j )
480  20 CONTINUE
481  END IF
482  END IF
483  alphar( j ) = h( j, j )
484  alphai( j ) = zero
485  beta( j ) = t( j, j )
486  30 CONTINUE
487 *
488 * If IHI < ILO, skip QZ steps
489 *
490  IF( ihi.LT.ilo )
491  $ GO TO 380
492 *
493 * MAIN QZ ITERATION LOOP
494 *
495 * Initialize dynamic indices
496 *
497 * Eigenvalues ILAST+1:N have been found.
498 * Column operations modify rows IFRSTM:whatever.
499 * Row operations modify columns whatever:ILASTM.
500 *
501 * If only eigenvalues are being computed, then
502 * IFRSTM is the row of the last splitting row above row ILAST;
503 * this is always at least ILO.
504 * IITER counts iterations since the last eigenvalue was found,
505 * to tell when to use an extraordinary shift.
506 * MAXIT is the maximum number of QZ sweeps allowed.
507 *
508  ilast = ihi
509  IF( ilschr ) THEN
510  ifrstm = 1
511  ilastm = n
512  ELSE
513  ifrstm = ilo
514  ilastm = ihi
515  END IF
516  iiter = 0
517  eshift = zero
518  maxit = 30*( ihi-ilo+1 )
519 *
520  DO 360 jiter = 1, maxit
521 *
522 * Split the matrix if possible.
523 *
524 * Two tests:
525 * 1: H(j,j-1)=0 or j=ILO
526 * 2: T(j,j)=0
527 *
528  IF( ilast.EQ.ilo ) THEN
529 *
530 * Special case: j=ILAST
531 *
532  GO TO 80
533  ELSE
534  IF( abs( h( ilast, ilast-1 ) ).LE.atol ) THEN
535  h( ilast, ilast-1 ) = zero
536  GO TO 80
537  END IF
538  END IF
539 *
540  IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
541  t( ilast, ilast ) = zero
542  GO TO 70
543  END IF
544 *
545 * General case: j<ILAST
546 *
547  DO 60 j = ilast - 1, ilo, -1
548 *
549 * Test 1: for H(j,j-1)=0 or j=ILO
550 *
551  IF( j.EQ.ilo ) THEN
552  ilazro = .true.
553  ELSE
554  IF( abs( h( j, j-1 ) ).LE.atol ) THEN
555  h( j, j-1 ) = zero
556  ilazro = .true.
557  ELSE
558  ilazro = .false.
559  END IF
560  END IF
561 *
562 * Test 2: for T(j,j)=0
563 *
564  IF( abs( t( j, j ) ).LT.btol ) THEN
565  t( j, j ) = zero
566 *
567 * Test 1a: Check for 2 consecutive small subdiagonals in A
568 *
569  ilazr2 = .false.
570  IF( .NOT.ilazro ) THEN
571  temp = abs( h( j, j-1 ) )
572  temp2 = abs( h( j, j ) )
573  tempr = max( temp, temp2 )
574  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
575  temp = temp / tempr
576  temp2 = temp2 / tempr
577  END IF
578  IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
579  $ ( ascale*atol ) )ilazr2 = .true.
580  END IF
581 *
582 * If both tests pass (1 & 2), i.e., the leading diagonal
583 * element of B in the block is zero, split a 1x1 block off
584 * at the top. (I.e., at the J-th row/column) The leading
585 * diagonal element of the remainder can also be zero, so
586 * this may have to be done repeatedly.
587 *
588  IF( ilazro .OR. ilazr2 ) THEN
589  DO 40 jch = j, ilast - 1
590  temp = h( jch, jch )
591  CALL dlartg( temp, h( jch+1, jch ), c, s,
592  $ h( jch, jch ) )
593  h( jch+1, jch ) = zero
594  CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
595  $ h( jch+1, jch+1 ), ldh, c, s )
596  CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
597  $ t( jch+1, jch+1 ), ldt, c, s )
598  IF( ilq )
599  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
600  $ c, s )
601  IF( ilazr2 )
602  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
603  ilazr2 = .false.
604  IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
605  IF( jch+1.GE.ilast ) THEN
606  GO TO 80
607  ELSE
608  ifirst = jch + 1
609  GO TO 110
610  END IF
611  END IF
612  t( jch+1, jch+1 ) = zero
613  40 CONTINUE
614  GO TO 70
615  ELSE
616 *
617 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
618 * Then process as in the case T(ILAST,ILAST)=0
619 *
620  DO 50 jch = j, ilast - 1
621  temp = t( jch, jch+1 )
622  CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
623  $ t( jch, jch+1 ) )
624  t( jch+1, jch+1 ) = zero
625  IF( jch.LT.ilastm-1 )
626  $ CALL drot( ilastm-jch-1, t( jch, jch+2 ), ldt,
627  $ t( jch+1, jch+2 ), ldt, c, s )
628  CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
629  $ h( jch+1, jch-1 ), ldh, c, s )
630  IF( ilq )
631  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
632  $ c, s )
633  temp = h( jch+1, jch )
634  CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
635  $ h( jch+1, jch ) )
636  h( jch+1, jch-1 ) = zero
637  CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
638  $ h( ifrstm, jch-1 ), 1, c, s )
639  CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
640  $ t( ifrstm, jch-1 ), 1, c, s )
641  IF( ilz )
642  $ CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
643  $ c, s )
644  50 CONTINUE
645  GO TO 70
646  END IF
647  ELSE IF( ilazro ) THEN
648 *
649 * Only test 1 passed -- work on J:ILAST
650 *
651  ifirst = j
652  GO TO 110
653  END IF
654 *
655 * Neither test passed -- try next J
656 *
657  60 CONTINUE
658 *
659 * (Drop-through is "impossible")
660 *
661  info = n + 1
662  GO TO 420
663 *
664 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
665 * 1x1 block.
666 *
667  70 CONTINUE
668  temp = h( ilast, ilast )
669  CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
670  $ h( ilast, ilast ) )
671  h( ilast, ilast-1 ) = zero
672  CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
673  $ h( ifrstm, ilast-1 ), 1, c, s )
674  CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
675  $ t( ifrstm, ilast-1 ), 1, c, s )
676  IF( ilz )
677  $ CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
678 *
679 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
680 * and BETA
681 *
682  80 CONTINUE
683  IF( t( ilast, ilast ).LT.zero ) THEN
684  IF( ilschr ) THEN
685  DO 90 j = ifrstm, ilast
686  h( j, ilast ) = -h( j, ilast )
687  t( j, ilast ) = -t( j, ilast )
688  90 CONTINUE
689  ELSE
690  h( ilast, ilast ) = -h( ilast, ilast )
691  t( ilast, ilast ) = -t( ilast, ilast )
692  END IF
693  IF( ilz ) THEN
694  DO 100 j = 1, n
695  z( j, ilast ) = -z( j, ilast )
696  100 CONTINUE
697  END IF
698  END IF
699  alphar( ilast ) = h( ilast, ilast )
700  alphai( ilast ) = zero
701  beta( ilast ) = t( ilast, ilast )
702 *
703 * Go to next block -- exit if finished.
704 *
705  ilast = ilast - 1
706  IF( ilast.LT.ilo )
707  $ GO TO 380
708 *
709 * Reset counters
710 *
711  iiter = 0
712  eshift = zero
713  IF( .NOT.ilschr ) THEN
714  ilastm = ilast
715  IF( ifrstm.GT.ilast )
716  $ ifrstm = ilo
717  END IF
718  GO TO 350
719 *
720 * QZ step
721 *
722 * This iteration only involves rows/columns IFIRST:ILAST. We
723 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
724 *
725  110 CONTINUE
726  iiter = iiter + 1
727  IF( .NOT.ilschr ) THEN
728  ifrstm = ifirst
729  END IF
730 *
731 * Compute single shifts.
732 *
733 * At this point, IFIRST < ILAST, and the diagonal elements of
734 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
735 * magnitude)
736 *
737  IF( ( iiter / 10 )*10.EQ.iiter ) THEN
738 *
739 * Exceptional shift. Chosen for no particularly good reason.
740 * (Single shift only.)
741 *
742  IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
743  $ abs( t( ilast-1, ilast-1 ) ) ) THEN
744  eshift = h( ilast, ilast-1 ) /
745  $ t( ilast-1, ilast-1 )
746  ELSE
747  eshift = eshift + one / ( safmin*dble( maxit ) )
748  END IF
749  s1 = one
750  wr = eshift
751 *
752  ELSE
753 *
754 * Shifts based on the generalized eigenvalues of the
755 * bottom-right 2x2 block of A and B. The first eigenvalue
756 * returned by DLAG2 is the Wilkinson shift (AEP p.512),
757 *
758  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
759  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
760  $ s2, wr, wr2, wi )
761 *
762  IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
763  $ .GT. abs( (wr2/s2)*t( ilast, ilast )
764  $ - h( ilast, ilast ) ) ) THEN
765  temp = wr
766  wr = wr2
767  wr2 = temp
768  temp = s1
769  s1 = s2
770  s2 = temp
771  END IF
772  temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
773  IF( wi.NE.zero )
774  $ GO TO 200
775  END IF
776 *
777 * Fiddle with shift to avoid overflow
778 *
779  temp = min( ascale, one )*( half*safmax )
780  IF( s1.GT.temp ) THEN
781  scale = temp / s1
782  ELSE
783  scale = one
784  END IF
785 *
786  temp = min( bscale, one )*( half*safmax )
787  IF( abs( wr ).GT.temp )
788  $ scale = min( scale, temp / abs( wr ) )
789  s1 = scale*s1
790  wr = scale*wr
791 *
792 * Now check for two consecutive small subdiagonals.
793 *
794  DO 120 j = ilast - 1, ifirst + 1, -1
795  istart = j
796  temp = abs( s1*h( j, j-1 ) )
797  temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
798  tempr = max( temp, temp2 )
799  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
800  temp = temp / tempr
801  temp2 = temp2 / tempr
802  END IF
803  IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
804  $ temp2 )GO TO 130
805  120 CONTINUE
806 *
807  istart = ifirst
808  130 CONTINUE
809 *
810 * Do an implicit single-shift QZ sweep.
811 *
812 * Initial Q
813 *
814  temp = s1*h( istart, istart ) - wr*t( istart, istart )
815  temp2 = s1*h( istart+1, istart )
816  CALL dlartg( temp, temp2, c, s, tempr )
817 *
818 * Sweep
819 *
820  DO 190 j = istart, ilast - 1
821  IF( j.GT.istart ) THEN
822  temp = h( j, j-1 )
823  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
824  h( j+1, j-1 ) = zero
825  END IF
826 *
827  DO 140 jc = j, ilastm
828  temp = c*h( j, jc ) + s*h( j+1, jc )
829  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
830  h( j, jc ) = temp
831  temp2 = c*t( j, jc ) + s*t( j+1, jc )
832  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
833  t( j, jc ) = temp2
834  140 CONTINUE
835  IF( ilq ) THEN
836  DO 150 jr = 1, n
837  temp = c*q( jr, j ) + s*q( jr, j+1 )
838  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
839  q( jr, j ) = temp
840  150 CONTINUE
841  END IF
842 *
843  temp = t( j+1, j+1 )
844  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
845  t( j+1, j ) = zero
846 *
847  DO 160 jr = ifrstm, min( j+2, ilast )
848  temp = c*h( jr, j+1 ) + s*h( jr, j )
849  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
850  h( jr, j+1 ) = temp
851  160 CONTINUE
852  DO 170 jr = ifrstm, j
853  temp = c*t( jr, j+1 ) + s*t( jr, j )
854  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
855  t( jr, j+1 ) = temp
856  170 CONTINUE
857  IF( ilz ) THEN
858  DO 180 jr = 1, n
859  temp = c*z( jr, j+1 ) + s*z( jr, j )
860  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
861  z( jr, j+1 ) = temp
862  180 CONTINUE
863  END IF
864  190 CONTINUE
865 *
866  GO TO 350
867 *
868 * Use Francis double-shift
869 *
870 * Note: the Francis double-shift should work with real shifts,
871 * but only if the block is at least 3x3.
872 * This code may break if this point is reached with
873 * a 2x2 block with real eigenvalues.
874 *
875  200 CONTINUE
876  IF( ifirst+1.EQ.ilast ) THEN
877 *
878 * Special case -- 2x2 block with complex eigenvectors
879 *
880 * Step 1: Standardize, that is, rotate so that
881 *
882 * ( B11 0 )
883 * B = ( ) with B11 non-negative.
884 * ( 0 B22 )
885 *
886  CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
887  $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
888 *
889  IF( b11.LT.zero ) THEN
890  cr = -cr
891  sr = -sr
892  b11 = -b11
893  b22 = -b22
894  END IF
895 *
896  CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
897  $ h( ilast, ilast-1 ), ldh, cl, sl )
898  CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
899  $ h( ifrstm, ilast ), 1, cr, sr )
900 *
901  IF( ilast.LT.ilastm )
902  $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
903  $ t( ilast, ilast+1 ), ldt, cl, sl )
904  IF( ifrstm.LT.ilast-1 )
905  $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
906  $ t( ifrstm, ilast ), 1, cr, sr )
907 *
908  IF( ilq )
909  $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
910  $ sl )
911  IF( ilz )
912  $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
913  $ sr )
914 *
915  t( ilast-1, ilast-1 ) = b11
916  t( ilast-1, ilast ) = zero
917  t( ilast, ilast-1 ) = zero
918  t( ilast, ilast ) = b22
919 *
920 * If B22 is negative, negate column ILAST
921 *
922  IF( b22.LT.zero ) THEN
923  DO 210 j = ifrstm, ilast
924  h( j, ilast ) = -h( j, ilast )
925  t( j, ilast ) = -t( j, ilast )
926  210 CONTINUE
927 *
928  IF( ilz ) THEN
929  DO 220 j = 1, n
930  z( j, ilast ) = -z( j, ilast )
931  220 CONTINUE
932  END IF
933  b22 = -b22
934  END IF
935 *
936 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
937 *
938 * Recompute shift
939 *
940  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
941  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
942  $ temp, wr, temp2, wi )
943 *
944 * If standardization has perturbed the shift onto real line,
945 * do another (real single-shift) QR step.
946 *
947  IF( wi.EQ.zero )
948  $ GO TO 350
949  s1inv = one / s1
950 *
951 * Do EISPACK (QZVAL) computation of alpha and beta
952 *
953  a11 = h( ilast-1, ilast-1 )
954  a21 = h( ilast, ilast-1 )
955  a12 = h( ilast-1, ilast )
956  a22 = h( ilast, ilast )
957 *
958 * Compute complex Givens rotation on right
959 * (Assume some element of C = (sA - wB) > unfl )
960 * __
961 * (sA - wB) ( CZ -SZ )
962 * ( SZ CZ )
963 *
964  c11r = s1*a11 - wr*b11
965  c11i = -wi*b11
966  c12 = s1*a12
967  c21 = s1*a21
968  c22r = s1*a22 - wr*b22
969  c22i = -wi*b22
970 *
971  IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
972  $ abs( c22r )+abs( c22i ) ) THEN
973  t1 = dlapy3( c12, c11r, c11i )
974  cz = c12 / t1
975  szr = -c11r / t1
976  szi = -c11i / t1
977  ELSE
978  cz = dlapy2( c22r, c22i )
979  IF( cz.LE.safmin ) THEN
980  cz = zero
981  szr = one
982  szi = zero
983  ELSE
984  tempr = c22r / cz
985  tempi = c22i / cz
986  t1 = dlapy2( cz, c21 )
987  cz = cz / t1
988  szr = -c21*tempr / t1
989  szi = c21*tempi / t1
990  END IF
991  END IF
992 *
993 * Compute Givens rotation on left
994 *
995 * ( CQ SQ )
996 * ( __ ) A or B
997 * ( -SQ CQ )
998 *
999  an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1000  bn = abs( b11 ) + abs( b22 )
1001  wabs = abs( wr ) + abs( wi )
1002  IF( s1*an.GT.wabs*bn ) THEN
1003  cq = cz*b11
1004  sqr = szr*b22
1005  sqi = -szi*b22
1006  ELSE
1007  a1r = cz*a11 + szr*a12
1008  a1i = szi*a12
1009  a2r = cz*a21 + szr*a22
1010  a2i = szi*a22
1011  cq = dlapy2( a1r, a1i )
1012  IF( cq.LE.safmin ) THEN
1013  cq = zero
1014  sqr = one
1015  sqi = zero
1016  ELSE
1017  tempr = a1r / cq
1018  tempi = a1i / cq
1019  sqr = tempr*a2r + tempi*a2i
1020  sqi = tempi*a2r - tempr*a2i
1021  END IF
1022  END IF
1023  t1 = dlapy3( cq, sqr, sqi )
1024  cq = cq / t1
1025  sqr = sqr / t1
1026  sqi = sqi / t1
1027 *
1028 * Compute diagonal elements of QBZ
1029 *
1030  tempr = sqr*szr - sqi*szi
1031  tempi = sqr*szi + sqi*szr
1032  b1r = cq*cz*b11 + tempr*b22
1033  b1i = tempi*b22
1034  b1a = dlapy2( b1r, b1i )
1035  b2r = cq*cz*b22 + tempr*b11
1036  b2i = -tempi*b11
1037  b2a = dlapy2( b2r, b2i )
1038 *
1039 * Normalize so beta > 0, and Im( alpha1 ) > 0
1040 *
1041  beta( ilast-1 ) = b1a
1042  beta( ilast ) = b2a
1043  alphar( ilast-1 ) = ( wr*b1a )*s1inv
1044  alphai( ilast-1 ) = ( wi*b1a )*s1inv
1045  alphar( ilast ) = ( wr*b2a )*s1inv
1046  alphai( ilast ) = -( wi*b2a )*s1inv
1047 *
1048 * Step 3: Go to next block -- exit if finished.
1049 *
1050  ilast = ifirst - 1
1051  IF( ilast.LT.ilo )
1052  $ GO TO 380
1053 *
1054 * Reset counters
1055 *
1056  iiter = 0
1057  eshift = zero
1058  IF( .NOT.ilschr ) THEN
1059  ilastm = ilast
1060  IF( ifrstm.GT.ilast )
1061  $ ifrstm = ilo
1062  END IF
1063  GO TO 350
1064  ELSE
1065 *
1066 * Usual case: 3x3 or larger block, using Francis implicit
1067 * double-shift
1068 *
1069 * 2
1070 * Eigenvalue equation is w - c w + d = 0,
1071 *
1072 * -1 2 -1
1073 * so compute 1st column of (A B ) - c A B + d
1074 * using the formula in QZIT (from EISPACK)
1075 *
1076 * We assume that the block is at least 3x3
1077 *
1078  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1079  $ ( bscale*t( ilast-1, ilast-1 ) )
1080  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1081  $ ( bscale*t( ilast-1, ilast-1 ) )
1082  ad12 = ( ascale*h( ilast-1, ilast ) ) /
1083  $ ( bscale*t( ilast, ilast ) )
1084  ad22 = ( ascale*h( ilast, ilast ) ) /
1085  $ ( bscale*t( ilast, ilast ) )
1086  u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1087  ad11l = ( ascale*h( ifirst, ifirst ) ) /
1088  $ ( bscale*t( ifirst, ifirst ) )
1089  ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1090  $ ( bscale*t( ifirst, ifirst ) )
1091  ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1092  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1093  ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1094  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1095  ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1096  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1097  u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1098 *
1099  v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1100  $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1101  v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1102  $ ( ad22-ad11l )+ad21*u12 )*ad21l
1103  v( 3 ) = ad32l*ad21l
1104 *
1105  istart = ifirst
1106 *
1107  CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1108  v( 1 ) = one
1109 *
1110 * Sweep
1111 *
1112  DO 290 j = istart, ilast - 2
1113 *
1114 * All but last elements: use 3x3 Householder transforms.
1115 *
1116 * Zero (j-1)st column of A
1117 *
1118  IF( j.GT.istart ) THEN
1119  v( 1 ) = h( j, j-1 )
1120  v( 2 ) = h( j+1, j-1 )
1121  v( 3 ) = h( j+2, j-1 )
1122 *
1123  CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1124  v( 1 ) = one
1125  h( j+1, j-1 ) = zero
1126  h( j+2, j-1 ) = zero
1127  END IF
1128 *
1129  DO 230 jc = j, ilastm
1130  temp = tau*( h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1131  $ h( j+2, jc ) )
1132  h( j, jc ) = h( j, jc ) - temp
1133  h( j+1, jc ) = h( j+1, jc ) - temp*v( 2 )
1134  h( j+2, jc ) = h( j+2, jc ) - temp*v( 3 )
1135  temp2 = tau*( t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1136  $ t( j+2, jc ) )
1137  t( j, jc ) = t( j, jc ) - temp2
1138  t( j+1, jc ) = t( j+1, jc ) - temp2*v( 2 )
1139  t( j+2, jc ) = t( j+2, jc ) - temp2*v( 3 )
1140  230 CONTINUE
1141  IF( ilq ) THEN
1142  DO 240 jr = 1, n
1143  temp = tau*( q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1144  $ q( jr, j+2 ) )
1145  q( jr, j ) = q( jr, j ) - temp
1146  q( jr, j+1 ) = q( jr, j+1 ) - temp*v( 2 )
1147  q( jr, j+2 ) = q( jr, j+2 ) - temp*v( 3 )
1148  240 CONTINUE
1149  END IF
1150 *
1151 * Zero j-th column of B (see DLAGBC for details)
1152 *
1153 * Swap rows to pivot
1154 *
1155  ilpivt = .false.
1156  temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1157  temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1158  IF( max( temp, temp2 ).LT.safmin ) THEN
1159  scale = zero
1160  u1 = one
1161  u2 = zero
1162  GO TO 250
1163  ELSE IF( temp.GE.temp2 ) THEN
1164  w11 = t( j+1, j+1 )
1165  w21 = t( j+2, j+1 )
1166  w12 = t( j+1, j+2 )
1167  w22 = t( j+2, j+2 )
1168  u1 = t( j+1, j )
1169  u2 = t( j+2, j )
1170  ELSE
1171  w21 = t( j+1, j+1 )
1172  w11 = t( j+2, j+1 )
1173  w22 = t( j+1, j+2 )
1174  w12 = t( j+2, j+2 )
1175  u2 = t( j+1, j )
1176  u1 = t( j+2, j )
1177  END IF
1178 *
1179 * Swap columns if nec.
1180 *
1181  IF( abs( w12 ).GT.abs( w11 ) ) THEN
1182  ilpivt = .true.
1183  temp = w12
1184  temp2 = w22
1185  w12 = w11
1186  w22 = w21
1187  w11 = temp
1188  w21 = temp2
1189  END IF
1190 *
1191 * LU-factor
1192 *
1193  temp = w21 / w11
1194  u2 = u2 - temp*u1
1195  w22 = w22 - temp*w12
1196  w21 = zero
1197 *
1198 * Compute SCALE
1199 *
1200  scale = one
1201  IF( abs( w22 ).LT.safmin ) THEN
1202  scale = zero
1203  u2 = one
1204  u1 = -w12 / w11
1205  GO TO 250
1206  END IF
1207  IF( abs( w22 ).LT.abs( u2 ) )
1208  $ scale = abs( w22 / u2 )
1209  IF( abs( w11 ).LT.abs( u1 ) )
1210  $ scale = min( scale, abs( w11 / u1 ) )
1211 *
1212 * Solve
1213 *
1214  u2 = ( scale*u2 ) / w22
1215  u1 = ( scale*u1-w12*u2 ) / w11
1216 *
1217  250 CONTINUE
1218  IF( ilpivt ) THEN
1219  temp = u2
1220  u2 = u1
1221  u1 = temp
1222  END IF
1223 *
1224 * Compute Householder Vector
1225 *
1226  t1 = sqrt( scale**2+u1**2+u2**2 )
1227  tau = one + scale / t1
1228  vs = -one / ( scale+t1 )
1229  v( 1 ) = one
1230  v( 2 ) = vs*u1
1231  v( 3 ) = vs*u2
1232 *
1233 * Apply transformations from the right.
1234 *
1235  DO 260 jr = ifrstm, min( j+3, ilast )
1236  temp = tau*( h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1237  $ h( jr, j+2 ) )
1238  h( jr, j ) = h( jr, j ) - temp
1239  h( jr, j+1 ) = h( jr, j+1 ) - temp*v( 2 )
1240  h( jr, j+2 ) = h( jr, j+2 ) - temp*v( 3 )
1241  260 CONTINUE
1242  DO 270 jr = ifrstm, j + 2
1243  temp = tau*( t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1244  $ t( jr, j+2 ) )
1245  t( jr, j ) = t( jr, j ) - temp
1246  t( jr, j+1 ) = t( jr, j+1 ) - temp*v( 2 )
1247  t( jr, j+2 ) = t( jr, j+2 ) - temp*v( 3 )
1248  270 CONTINUE
1249  IF( ilz ) THEN
1250  DO 280 jr = 1, n
1251  temp = tau*( z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1252  $ z( jr, j+2 ) )
1253  z( jr, j ) = z( jr, j ) - temp
1254  z( jr, j+1 ) = z( jr, j+1 ) - temp*v( 2 )
1255  z( jr, j+2 ) = z( jr, j+2 ) - temp*v( 3 )
1256  280 CONTINUE
1257  END IF
1258  t( j+1, j ) = zero
1259  t( j+2, j ) = zero
1260  290 CONTINUE
1261 *
1262 * Last elements: Use Givens rotations
1263 *
1264 * Rotations from the left
1265 *
1266  j = ilast - 1
1267  temp = h( j, j-1 )
1268  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1269  h( j+1, j-1 ) = zero
1270 *
1271  DO 300 jc = j, ilastm
1272  temp = c*h( j, jc ) + s*h( j+1, jc )
1273  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1274  h( j, jc ) = temp
1275  temp2 = c*t( j, jc ) + s*t( j+1, jc )
1276  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1277  t( j, jc ) = temp2
1278  300 CONTINUE
1279  IF( ilq ) THEN
1280  DO 310 jr = 1, n
1281  temp = c*q( jr, j ) + s*q( jr, j+1 )
1282  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1283  q( jr, j ) = temp
1284  310 CONTINUE
1285  END IF
1286 *
1287 * Rotations from the right.
1288 *
1289  temp = t( j+1, j+1 )
1290  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1291  t( j+1, j ) = zero
1292 *
1293  DO 320 jr = ifrstm, ilast
1294  temp = c*h( jr, j+1 ) + s*h( jr, j )
1295  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1296  h( jr, j+1 ) = temp
1297  320 CONTINUE
1298  DO 330 jr = ifrstm, ilast - 1
1299  temp = c*t( jr, j+1 ) + s*t( jr, j )
1300  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1301  t( jr, j+1 ) = temp
1302  330 CONTINUE
1303  IF( ilz ) THEN
1304  DO 340 jr = 1, n
1305  temp = c*z( jr, j+1 ) + s*z( jr, j )
1306  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1307  z( jr, j+1 ) = temp
1308  340 CONTINUE
1309  END IF
1310 *
1311 * End of Double-Shift code
1312 *
1313  END IF
1314 *
1315  GO TO 350
1316 *
1317 * End of iteration loop
1318 *
1319  350 CONTINUE
1320  360 CONTINUE
1321 *
1322 * Drop-through = non-convergence
1323 *
1324  info = ilast
1325  GO TO 420
1326 *
1327 * Successful completion of all QZ steps
1328 *
1329  380 CONTINUE
1330 *
1331 * Set Eigenvalues 1:ILO-1
1332 *
1333  DO 410 j = 1, ilo - 1
1334  IF( t( j, j ).LT.zero ) THEN
1335  IF( ilschr ) THEN
1336  DO 390 jr = 1, j
1337  h( jr, j ) = -h( jr, j )
1338  t( jr, j ) = -t( jr, j )
1339  390 CONTINUE
1340  ELSE
1341  h( j, j ) = -h( j, j )
1342  t( j, j ) = -t( j, j )
1343  END IF
1344  IF( ilz ) THEN
1345  DO 400 jr = 1, n
1346  z( jr, j ) = -z( jr, j )
1347  400 CONTINUE
1348  END IF
1349  END IF
1350  alphar( j ) = h( j, j )
1351  alphai( j ) = zero
1352  beta( j ) = t( j, j )
1353  410 CONTINUE
1354 *
1355 * Normal Termination
1356 *
1357  info = 0
1358 *
1359 * Exit (other than argument error) -- return optimal workspace size
1360 *
1361  420 CONTINUE
1362  work( 1 ) = dble( n )
1363  RETURN
1364 *
1365 * End of DHGEQZ
1366 *
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:112
double precision function dlanhs(NORM, N, A, LDA, WORK)
DLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlanhs.f:110
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:158
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:108
subroutine dlasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: dlasv2.f:140
double precision function dlapy3(X, Y, Z)
DLAPY3 returns sqrt(x2+y2+z2).
Definition: dlapy3.f:70
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: