dhgeqz.f

Go to the documentation of this file.

*> \brief \b DHGEQZ
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DHGEQZ + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhgeqz.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhgeqz.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhgeqz.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
*                          ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
*                          LWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          COMPQ, COMPZ, JOB
*       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ),
*      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
*      $                   WORK( * ), Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOB
*> \verbatim
*>          JOB is CHARACTER*1
*>          = 'E': Compute eigenvalues only;
*>          = 'S': Compute eigenvalues and the Schur form.
*> \endverbatim
*>
*> \param[in] COMPQ
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] COMPZ
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices H, T, Q, and Z.  N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in,out] H
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDH
*> \verbatim
*>          LDH is INTEGER
*>          The leading dimension of the array H.  LDH >= max( 1, N ).
*> \endverbatim
*>
*> \param[in,out] T
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDT
*> \verbatim
*>          LDT is INTEGER
*>          The leading dimension of the array T.  LDT >= max( 1, N ).
*> \endverbatim
*>
*> \param[out] ALPHAR
*> \verbatim
*>          ALPHAR is DOUBLE PRECISION array, dimension (N)
*>          The real parts of each scalar alpha defining an eigenvalue
*>          of GNEP.
*> \endverbatim
*>
*> \param[out] ALPHAI
*> \verbatim
*>          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).
*> \endverbatim
*>
*> \param[out] BETA
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          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'.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q.  LDQ >= 1.
*>          If COMPQ='V' or 'I', then LDQ >= N.
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          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'.
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>          The leading dimension of the array Z.  LDZ >= 1.
*>          If COMPZ='V' or 'I', then LDZ >= N.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
*>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          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.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hgeqz
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  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.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE dhgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
     $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
     $                   LWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          COMPQ, COMPZ, JOB
      INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ),
     $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
     $                   work( * ), z( ldz, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
*    $                     SAFETY = 1.0E+0 )
      DOUBLE PRECISION   HALF, ZERO, ONE, SAFETY
      PARAMETER          ( HALF = 0.5d+0, zero = 0.0d+0, one = 1.0d+0,
     $                   safety = 1.0d+2 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
     $                   LQUERY
      INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
     $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
     $                   jr, maxit
      DOUBLE PRECISION   A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
     $                   AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
     $                   ad32l, an, anorm, ascale, atol, b11, b1a, b1i,
     $                   b1r, b22, b2a, b2i, b2r, bn, bnorm, bscale,
     $                   btol, c, c11i, c11r, c12, c21, c22i, c22r, cl,
     $                   cq, cr, cz, eshift, s, s1, s1inv, s2, safmax,
     $                   safmin, scale, sl, sqi, sqr, sr, szi, szr, t1,
     $                   t2, t3, tau, temp, temp2, tempi, tempr, u1,
     $                   u12, u12l, u2, ulp, vs, w11, w12, w21, w22,
     $                   wabs, wi, wr, wr2
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   V( 3 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, DLANHS, DLAPY2, DLAPY3
      EXTERNAL           lsame, dlamch, dlanhs, dlapy2, dlapy3
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlag2, dlarfg, dlartg, dlaset, dlasv2, drot,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, max, min, sqrt
*     ..
*     .. Executable Statements ..
*
*     Decode JOB, COMPQ, COMPZ
*
      IF( lsame( job, 'E' ) ) THEN
         ilschr = .false.
         ischur = 1
      ELSE IF( lsame( job, 'S' ) ) THEN
         ilschr = .true.
         ischur = 2
      ELSE
         ischur = 0
      END IF
*
      IF( lsame( compq, 'N' ) ) THEN
         ilq = .false.
         icompq = 1
      ELSE IF( lsame( compq, 'V' ) ) THEN
         ilq = .true.
         icompq = 2
      ELSE IF( lsame( compq, 'I' ) ) THEN
         ilq = .true.
         icompq = 3
      ELSE
         icompq = 0
      END IF
*
      IF( lsame( compz, 'N' ) ) THEN
         ilz = .false.
         icompz = 1
      ELSE IF( lsame( compz, 'V' ) ) THEN
         ilz = .true.
         icompz = 2
      ELSE IF( lsame( compz, 'I' ) ) THEN
         ilz = .true.
         icompz = 3
      ELSE
         icompz = 0
      END IF
*
*     Check Argument Values
*
      info = 0
      work( 1 ) = max( 1, n )
      lquery = ( lwork.EQ.-1 )
      IF( ischur.EQ.0 ) THEN
         info = -1
      ELSE IF( icompq.EQ.0 ) THEN
         info = -2
      ELSE IF( icompz.EQ.0 ) THEN
         info = -3
      ELSE IF( n.LT.0 ) THEN
         info = -4
      ELSE IF( ilo.LT.1 ) THEN
         info = -5
      ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
         info = -6
      ELSE IF( ldh.LT.n ) THEN
         info = -8
      ELSE IF( ldt.LT.n ) THEN
         info = -10
      ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
         info = -15
      ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
         info = -17
      ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
         info = -19
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DHGEQZ', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( n.LE.0 ) THEN
         work( 1 ) = dble( 1 )
         RETURN
      END IF
*
*     Initialize Q and Z
*
      IF( icompq.EQ.3 )
     $   CALL dlaset( 'Full', n, n, zero, one, q, ldq )
      IF( icompz.EQ.3 )
     $   CALL dlaset( 'Full', n, n, zero, one, z, ldz )
*
*     Machine Constants
*
      in = ihi + 1 - ilo
      safmin = dlamch( 'S' )
      safmax = one / safmin
      ulp = dlamch( 'E' )*dlamch( 'B' )
      anorm = dlanhs( 'F', in, h( ilo, ilo ), ldh, work )
      bnorm = dlanhs( 'F', in, t( ilo, ilo ), ldt, work )
      atol = max( safmin, ulp*anorm )
      btol = max( safmin, ulp*bnorm )
      ascale = one / max( safmin, anorm )
      bscale = one / max( safmin, bnorm )
*
*     Set Eigenvalues IHI+1:N
*
      DO 30 j = ihi + 1, n
         IF( t( j, j ).LT.zero ) THEN
            IF( ilschr ) THEN
               DO 10 jr = 1, j
                  h( jr, j ) = -h( jr, j )
                  t( jr, j ) = -t( jr, j )
   10          CONTINUE
            ELSE
               h( j, j ) = -h( j, j )
               t( j, j ) = -t( j, j )
            END IF
            IF( ilz ) THEN
               DO 20 jr = 1, n
                  z( jr, j ) = -z( jr, j )
   20          CONTINUE
            END IF
         END IF
         alphar( j ) = h( j, j )
         alphai( j ) = zero
         beta( j ) = t( j, j )
   30 CONTINUE
*
*     If IHI < ILO, skip QZ steps
*
      IF( ihi.LT.ilo )
     $   GO TO 380
*
*     MAIN QZ ITERATION LOOP
*
*     Initialize dynamic indices
*
*     Eigenvalues ILAST+1:N have been found.
*        Column operations modify rows IFRSTM:whatever.
*        Row operations modify columns whatever:ILASTM.
*
*     If only eigenvalues are being computed, then
*        IFRSTM is the row of the last splitting row above row ILAST;
*        this is always at least ILO.
*     IITER counts iterations since the last eigenvalue was found,
*        to tell when to use an extraordinary shift.
*     MAXIT is the maximum number of QZ sweeps allowed.
*
      ilast = ihi
      IF( ilschr ) THEN
         ifrstm = 1
         ilastm = n
      ELSE
         ifrstm = ilo
         ilastm = ihi
      END IF
      iiter = 0
      eshift = zero
      maxit = 30*( ihi-ilo+1 )
*
      DO 360 jiter = 1, maxit
*
*        Split the matrix if possible.
*
*        Two tests:
*           1: H(j,j-1)=0  or  j=ILO
*           2: T(j,j)=0
*
         IF( ilast.EQ.ilo ) THEN
*
*           Special case: j=ILAST
*
            GO TO 80
         ELSE
            IF( abs( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*( 
     $         abs( h( ilast, ilast ) ) + abs( h( ilast-1, ilast-1 ) ) 
     $         ) ) ) THEN
               h( ilast, ilast-1 ) = zero
               GO TO 80
            END IF
         END IF
*
         IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
            t( ilast, ilast ) = zero
            GO TO 70
         END IF
*
*        General case: j<ILAST
*
         DO 60 j = ilast - 1, ilo, -1
*
*           Test 1: for H(j,j-1)=0 or j=ILO
*
            IF( j.EQ.ilo ) THEN
               ilazro = .true.
            ELSE
               IF( abs( h( j, j-1 ) ).LE.max( safmin, ulp*( 
     $         abs( h( j, j ) ) + abs( h( j-1, j-1 ) ) 
     $         ) ) ) THEN
                  h( j, j-1 ) = zero
                  ilazro = .true.
               ELSE
                  ilazro = .false.
               END IF
            END IF
*
*           Test 2: for T(j,j)=0
*
            IF( abs( t( j, j ) ).LT.btol ) THEN
               t( j, j ) = zero
*
*              Test 1a: Check for 2 consecutive small subdiagonals in A
*
               ilazr2 = .false.
               IF( .NOT.ilazro ) THEN
                  temp = abs( h( j, j-1 ) )
                  temp2 = abs( h( j, j ) )
                  tempr = max( temp, temp2 )
                  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
                     temp = temp / tempr
                     temp2 = temp2 / tempr
                  END IF
                  IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
     $                ( ascale*atol ) )ilazr2 = .true.
               END IF
*
*              If both tests pass (1 & 2), i.e., the leading diagonal
*              element of B in the block is zero, split a 1x1 block off
*              at the top. (I.e., at the J-th row/column) The leading
*              diagonal element of the remainder can also be zero, so
*              this may have to be done repeatedly.
*
               IF( ilazro .OR. ilazr2 ) THEN
                  DO 40 jch = j, ilast - 1
                     temp = h( jch, jch )
                     CALL dlartg( temp, h( jch+1, jch ), c, s,
     $                            h( jch, jch ) )
                     h( jch+1, jch ) = zero
                     CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
     $                          h( jch+1, jch+1 ), ldh, c, s )
                     CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
     $                          t( jch+1, jch+1 ), ldt, c, s )
                     IF( ilq )
     $                  CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
     $                             c, s )
                     IF( ilazr2 )
     $                  h( jch, jch-1 ) = h( jch, jch-1 )*c
                     ilazr2 = .false.
                     IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
                        IF( jch+1.GE.ilast ) THEN
                           GO TO 80
                        ELSE
                           ifirst = jch + 1
                           GO TO 110
                        END IF
                     END IF
                     t( jch+1, jch+1 ) = zero
   40             CONTINUE
                  GO TO 70
               ELSE
*
*                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
*                 Then process as in the case T(ILAST,ILAST)=0
*
                  DO 50 jch = j, ilast - 1
                     temp = t( jch, jch+1 )
                     CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
     $                            t( jch, jch+1 ) )
                     t( jch+1, jch+1 ) = zero
                     IF( jch.LT.ilastm-1 )
     $                  CALL drot( ilastm-jch-1, t( jch, jch+2 ), ldt,
     $                             t( jch+1, jch+2 ), ldt, c, s )
                     CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
     $                          h( jch+1, jch-1 ), ldh, c, s )
                     IF( ilq )
     $                  CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
     $                             c, s )
                     temp = h( jch+1, jch )
                     CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
     $                            h( jch+1, jch ) )
                     h( jch+1, jch-1 ) = zero
                     CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
     $                          h( ifrstm, jch-1 ), 1, c, s )
                     CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
     $                          t( ifrstm, jch-1 ), 1, c, s )
                     IF( ilz )
     $                  CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
     $                             c, s )
   50             CONTINUE
                  GO TO 70
               END IF
            ELSE IF( ilazro ) THEN
*
*              Only test 1 passed -- work on J:ILAST
*
               ifirst = j
               GO TO 110
            END IF
*
*           Neither test passed -- try next J
*
   60    CONTINUE
*
*        (Drop-through is "impossible")
*
         info = n + 1
         GO TO 420
*
*        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
*        1x1 block.
*
   70    CONTINUE
         temp = h( ilast, ilast )
         CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
     $                h( ilast, ilast ) )
         h( ilast, ilast-1 ) = zero
         CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
     $              h( ifrstm, ilast-1 ), 1, c, s )
         CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
     $              t( ifrstm, ilast-1 ), 1, c, s )
         IF( ilz )
     $      CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
*
*        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
*                              and BETA
*
   80    CONTINUE
         IF( t( ilast, ilast ).LT.zero ) THEN
            IF( ilschr ) THEN
               DO 90 j = ifrstm, ilast
                  h( j, ilast ) = -h( j, ilast )
                  t( j, ilast ) = -t( j, ilast )
   90          CONTINUE
            ELSE
               h( ilast, ilast ) = -h( ilast, ilast )
               t( ilast, ilast ) = -t( ilast, ilast )
            END IF
            IF( ilz ) THEN
               DO 100 j = 1, n
                  z( j, ilast ) = -z( j, ilast )
  100          CONTINUE
            END IF
         END IF
         alphar( ilast ) = h( ilast, ilast )
         alphai( ilast ) = zero
         beta( ilast ) = t( ilast, ilast )
*
*        Go to next block -- exit if finished.
*
         ilast = ilast - 1
         IF( ilast.LT.ilo )
     $      GO TO 380
*
*        Reset counters
*
         iiter = 0
         eshift = zero
         IF( .NOT.ilschr ) THEN
            ilastm = ilast
            IF( ifrstm.GT.ilast )
     $         ifrstm = ilo
         END IF
         GO TO 350
*
*        QZ step
*
*        This iteration only involves rows/columns IFIRST:ILAST. We
*        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
*
  110    CONTINUE
         iiter = iiter + 1
         IF( .NOT.ilschr ) THEN
            ifrstm = ifirst
         END IF
*
*        Compute single shifts.
*
*        At this point, IFIRST < ILAST, and the diagonal elements of
*        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
*        magnitude)
*
         IF( ( iiter / 10 )*10.EQ.iiter ) THEN
*
*           Exceptional shift.  Chosen for no particularly good reason.
*           (Single shift only.)
*
            IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
     $          abs( t( ilast-1, ilast-1 ) ) ) THEN
               eshift = h( ilast, ilast-1 ) /
     $                  t( ilast-1, ilast-1 )
            ELSE
               eshift = eshift + one / ( safmin*dble( maxit ) )
            END IF
            s1 = one
            wr = eshift
*
         ELSE
*
*           Shifts based on the generalized eigenvalues of the
*           bottom-right 2x2 block of A and B. The first eigenvalue
*           returned by DLAG2 is the Wilkinson shift (AEP p.512),
*
            CALL dlag2( h( ilast-1, ilast-1 ), ldh,
     $                  t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
     $                  s2, wr, wr2, wi )
*
            IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
     $         .GT. abs( (wr2/s2)*t( ilast, ilast )
     $         - h( ilast, ilast ) ) ) THEN
               temp = wr
               wr = wr2
               wr2 = temp
               temp = s1
               s1 = s2
               s2 = temp
            END IF
            temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
            IF( wi.NE.zero )
     $         GO TO 200
         END IF
*
*        Fiddle with shift to avoid overflow
*
         temp = min( ascale, one )*( half*safmax )
         IF( s1.GT.temp ) THEN
            scale = temp / s1
         ELSE
            scale = one
         END IF
*
         temp = min( bscale, one )*( half*safmax )
         IF( abs( wr ).GT.temp )
     $      scale = min( scale, temp / abs( wr ) )
         s1 = scale*s1
         wr = scale*wr
*
*        Now check for two consecutive small subdiagonals.
*
         DO 120 j = ilast - 1, ifirst + 1, -1
            istart = j
            temp = abs( s1*h( j, j-1 ) )
            temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
            tempr = max( temp, temp2 )
            IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
               temp = temp / tempr
               temp2 = temp2 / tempr
            END IF
            IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
     $          temp2 )GO TO 130
  120    CONTINUE
*
         istart = ifirst
  130    CONTINUE
*
*        Do an implicit single-shift QZ sweep.
*
*        Initial Q
*
         temp = s1*h( istart, istart ) - wr*t( istart, istart )
         temp2 = s1*h( istart+1, istart )
         CALL dlartg( temp, temp2, c, s, tempr )
*
*        Sweep
*
         DO 190 j = istart, ilast - 1
            IF( j.GT.istart ) THEN
               temp = h( j, j-1 )
               CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
               h( j+1, j-1 ) = zero
            END IF
*
            DO 140 jc = j, ilastm
               temp = c*h( j, jc ) + s*h( j+1, jc )
               h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
               h( j, jc ) = temp
               temp2 = c*t( j, jc ) + s*t( j+1, jc )
               t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
               t( j, jc ) = temp2
  140       CONTINUE
            IF( ilq ) THEN
               DO 150 jr = 1, n
                  temp = c*q( jr, j ) + s*q( jr, j+1 )
                  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
                  q( jr, j ) = temp
  150          CONTINUE
            END IF
*
            temp = t( j+1, j+1 )
            CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
            t( j+1, j ) = zero
*
            DO 160 jr = ifrstm, min( j+2, ilast )
               temp = c*h( jr, j+1 ) + s*h( jr, j )
               h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
               h( jr, j+1 ) = temp
  160       CONTINUE
            DO 170 jr = ifrstm, j
               temp = c*t( jr, j+1 ) + s*t( jr, j )
               t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
               t( jr, j+1 ) = temp
  170       CONTINUE
            IF( ilz ) THEN
               DO 180 jr = 1, n
                  temp = c*z( jr, j+1 ) + s*z( jr, j )
                  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
                  z( jr, j+1 ) = temp
  180          CONTINUE
            END IF
  190    CONTINUE
*
         GO TO 350
*
*        Use Francis double-shift
*
*        Note: the Francis double-shift should work with real shifts,
*              but only if the block is at least 3x3.
*              This code may break if this point is reached with
*              a 2x2 block with real eigenvalues.
*
  200    CONTINUE
         IF( ifirst+1.EQ.ilast ) THEN
*
*           Special case -- 2x2 block with complex eigenvectors
*
*           Step 1: Standardize, that is, rotate so that
*
*                       ( B11  0  )
*                   B = (         )  with B11 non-negative.
*                       (  0  B22 )
*
            CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
     $                   t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
*
            IF( b11.LT.zero ) THEN
               cr = -cr
               sr = -sr
               b11 = -b11
               b22 = -b22
            END IF
*
            CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
     $                 h( ilast, ilast-1 ), ldh, cl, sl )
            CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
     $                 h( ifrstm, ilast ), 1, cr, sr )
*
            IF( ilast.LT.ilastm )
     $         CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
     $                    t( ilast, ilast+1 ), ldt, cl, sl )
            IF( ifrstm.LT.ilast-1 )
     $         CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
     $                    t( ifrstm, ilast ), 1, cr, sr )
*
            IF( ilq )
     $         CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
     $                    sl )
            IF( ilz )
     $         CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
     $                    sr )
*
            t( ilast-1, ilast-1 ) = b11
            t( ilast-1, ilast ) = zero
            t( ilast, ilast-1 ) = zero
            t( ilast, ilast ) = b22
*
*           If B22 is negative, negate column ILAST
*
            IF( b22.LT.zero ) THEN
               DO 210 j = ifrstm, ilast
                  h( j, ilast ) = -h( j, ilast )
                  t( j, ilast ) = -t( j, ilast )
  210          CONTINUE
*
               IF( ilz ) THEN
                  DO 220 j = 1, n
                     z( j, ilast ) = -z( j, ilast )
  220             CONTINUE
               END IF
               b22 = -b22
            END IF
*
*           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
*
*           Recompute shift
*
            CALL dlag2( h( ilast-1, ilast-1 ), ldh,
     $                  t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
     $                  temp, wr, temp2, wi )
*
*           If standardization has perturbed the shift onto real line,
*           do another (real single-shift) QR step.
*
            IF( wi.EQ.zero )
     $         GO TO 350
            s1inv = one / s1
*
*           Do EISPACK (QZVAL) computation of alpha and beta
*
            a11 = h( ilast-1, ilast-1 )
            a21 = h( ilast, ilast-1 )
            a12 = h( ilast-1, ilast )
            a22 = h( ilast, ilast )
*
*           Compute complex Givens rotation on right
*           (Assume some element of C = (sA - wB) > unfl )
*                            __
*           (sA - wB) ( CZ   -SZ )
*                     ( SZ    CZ )
*
            c11r = s1*a11 - wr*b11
            c11i = -wi*b11
            c12 = s1*a12
            c21 = s1*a21
            c22r = s1*a22 - wr*b22
            c22i = -wi*b22
*
            IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
     $          abs( c22r )+abs( c22i ) ) THEN
               t1 = dlapy3( c12, c11r, c11i )
               cz = c12 / t1
               szr = -c11r / t1
               szi = -c11i / t1
            ELSE
               cz = dlapy2( c22r, c22i )
               IF( cz.LE.safmin ) THEN
                  cz = zero
                  szr = one
                  szi = zero
               ELSE
                  tempr = c22r / cz
                  tempi = c22i / cz
                  t1 = dlapy2( cz, c21 )
                  cz = cz / t1
                  szr = -c21*tempr / t1
                  szi = c21*tempi / t1
               END IF
            END IF
*
*           Compute Givens rotation on left
*
*           (  CQ   SQ )
*           (  __      )  A or B
*           ( -SQ   CQ )
*
            an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
            bn = abs( b11 ) + abs( b22 )
            wabs = abs( wr ) + abs( wi )
            IF( s1*an.GT.wabs*bn ) THEN
               cq = cz*b11
               sqr = szr*b22
               sqi = -szi*b22
            ELSE
               a1r = cz*a11 + szr*a12
               a1i = szi*a12
               a2r = cz*a21 + szr*a22
               a2i = szi*a22
               cq = dlapy2( a1r, a1i )
               IF( cq.LE.safmin ) THEN
                  cq = zero
                  sqr = one
                  sqi = zero
               ELSE
                  tempr = a1r / cq
                  tempi = a1i / cq
                  sqr = tempr*a2r + tempi*a2i
                  sqi = tempi*a2r - tempr*a2i
               END IF
            END IF
            t1 = dlapy3( cq, sqr, sqi )
            cq = cq / t1
            sqr = sqr / t1
            sqi = sqi / t1
*
*           Compute diagonal elements of QBZ
*
            tempr = sqr*szr - sqi*szi
            tempi = sqr*szi + sqi*szr
            b1r = cq*cz*b11 + tempr*b22
            b1i = tempi*b22
            b1a = dlapy2( b1r, b1i )
            b2r = cq*cz*b22 + tempr*b11
            b2i = -tempi*b11
            b2a = dlapy2( b2r, b2i )
*
*           Normalize so beta > 0, and Im( alpha1 ) > 0
*
            beta( ilast-1 ) = b1a
            beta( ilast ) = b2a
            alphar( ilast-1 ) = ( wr*b1a )*s1inv
            alphai( ilast-1 ) = ( wi*b1a )*s1inv
            alphar( ilast ) = ( wr*b2a )*s1inv
            alphai( ilast ) = -( wi*b2a )*s1inv
*
*           Step 3: Go to next block -- exit if finished.
*
            ilast = ifirst - 1
            IF( ilast.LT.ilo )
     $         GO TO 380
*
*           Reset counters
*
            iiter = 0
            eshift = zero
            IF( .NOT.ilschr ) THEN
               ilastm = ilast
               IF( ifrstm.GT.ilast )
     $            ifrstm = ilo
            END IF
            GO TO 350
         ELSE
*
*           Usual case: 3x3 or larger block, using Francis implicit
*                       double-shift
*
*                                    2
*           Eigenvalue equation is  w  - c w + d = 0,
*
*                                         -1 2        -1
*           so compute 1st column of  (A B  )  - c A B   + d
*           using the formula in QZIT (from EISPACK)
*
*           We assume that the block is at least 3x3
*
            ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
     $             ( bscale*t( ilast-1, ilast-1 ) )
            ad21 = ( ascale*h( ilast, ilast-1 ) ) /
     $             ( bscale*t( ilast-1, ilast-1 ) )
            ad12 = ( ascale*h( ilast-1, ilast ) ) /
     $             ( bscale*t( ilast, ilast ) )
            ad22 = ( ascale*h( ilast, ilast ) ) /
     $             ( bscale*t( ilast, ilast ) )
            u12 = t( ilast-1, ilast ) / t( ilast, ilast )
            ad11l = ( ascale*h( ifirst, ifirst ) ) /
     $              ( bscale*t( ifirst, ifirst ) )
            ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
     $              ( bscale*t( ifirst, ifirst ) )
            ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
     $              ( bscale*t( ifirst+1, ifirst+1 ) )
            ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
     $              ( bscale*t( ifirst+1, ifirst+1 ) )
            ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
     $              ( bscale*t( ifirst+1, ifirst+1 ) )
            u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
*
            v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
     $               ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
            v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
     $               ( ad22-ad11l )+ad21*u12 )*ad21l
            v( 3 ) = ad32l*ad21l
*
            istart = ifirst
*
            CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
            v( 1 ) = one
*
*           Sweep
*
            DO 290 j = istart, ilast - 2
*
*              All but last elements: use 3x3 Householder transforms.
*
*              Zero (j-1)st column of A
*
               IF( j.GT.istart ) THEN
                  v( 1 ) = h( j, j-1 )
                  v( 2 ) = h( j+1, j-1 )
                  v( 3 ) = h( j+2, j-1 )
*
                  CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
                  v( 1 ) = one
                  h( j+1, j-1 ) = zero
                  h( j+2, j-1 ) = zero
               END IF
*
               t2 = tau*v( 2 )
               t3 = tau*v( 3 )
               DO 230 jc = j, ilastm
                  temp = h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
     $                   h( j+2, jc )
                  h( j, jc ) = h( j, jc ) - temp*tau
                  h( j+1, jc ) = h( j+1, jc ) - temp*t2
                  h( j+2, jc ) = h( j+2, jc ) - temp*t3
                  temp2 = t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
     $                    t( j+2, jc )
                  t( j, jc ) = t( j, jc ) - temp2*tau
                  t( j+1, jc ) = t( j+1, jc ) - temp2*t2
                  t( j+2, jc ) = t( j+2, jc ) - temp2*t3
  230          CONTINUE
               IF( ilq ) THEN
                  DO 240 jr = 1, n
                     temp = q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
     $                      q( jr, j+2 )
                     q( jr, j ) = q( jr, j ) - temp*tau
                     q( jr, j+1 ) = q( jr, j+1 ) - temp*t2
                     q( jr, j+2 ) = q( jr, j+2 ) - temp*t3
  240             CONTINUE
               END IF
*
*              Zero j-th column of B (see DLAGBC for details)
*
*              Swap rows to pivot
*
               ilpivt = .false.
               temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
               temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
               IF( max( temp, temp2 ).LT.safmin ) THEN
                  scale = zero
                  u1 = one
                  u2 = zero
                  GO TO 250
               ELSE IF( temp.GE.temp2 ) THEN
                  w11 = t( j+1, j+1 )
                  w21 = t( j+2, j+1 )
                  w12 = t( j+1, j+2 )
                  w22 = t( j+2, j+2 )
                  u1 = t( j+1, j )
                  u2 = t( j+2, j )
               ELSE
                  w21 = t( j+1, j+1 )
                  w11 = t( j+2, j+1 )
                  w22 = t( j+1, j+2 )
                  w12 = t( j+2, j+2 )
                  u2 = t( j+1, j )
                  u1 = t( j+2, j )
               END IF
*
*              Swap columns if nec.
*
               IF( abs( w12 ).GT.abs( w11 ) ) THEN
                  ilpivt = .true.
                  temp = w12
                  temp2 = w22
                  w12 = w11
                  w22 = w21
                  w11 = temp
                  w21 = temp2
               END IF
*
*              LU-factor
*
               temp = w21 / w11
               u2 = u2 - temp*u1
               w22 = w22 - temp*w12
               w21 = zero
*
*              Compute SCALE
*
               scale = one
               IF( abs( w22 ).LT.safmin ) THEN
                  scale = zero
                  u2 = one
                  u1 = -w12 / w11
                  GO TO 250
               END IF
               IF( abs( w22 ).LT.abs( u2 ) )
     $            scale = abs( w22 / u2 )
               IF( abs( w11 ).LT.abs( u1 ) )
     $            scale = min( scale, abs( w11 / u1 ) )
*
*              Solve
*
               u2 = ( scale*u2 ) / w22
               u1 = ( scale*u1-w12*u2 ) / w11
*
  250          CONTINUE
               IF( ilpivt ) THEN
                  temp = u2
                  u2 = u1
                  u1 = temp
               END IF
*
*              Compute Householder Vector
*
               t1 = sqrt( scale**2+u1**2+u2**2 )
               tau = one + scale / t1
               vs = -one / ( scale+t1 )
               v( 1 ) = one
               v( 2 ) = vs*u1
               v( 3 ) = vs*u2
*
*              Apply transformations from the right.
*
               t2 = tau*v(2)
               t3 = tau*v(3)
               DO 260 jr = ifrstm, min( j+3, ilast )
                  temp = h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
     $                   h( jr, j+2 )
                  h( jr, j ) = h( jr, j ) - temp*tau
                  h( jr, j+1 ) = h( jr, j+1 ) - temp*t2
                  h( jr, j+2 ) = h( jr, j+2 ) - temp*t3
  260          CONTINUE
               DO 270 jr = ifrstm, j + 2
                  temp = t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
     $                   t( jr, j+2 )
                  t( jr, j ) = t( jr, j ) - temp*tau
                  t( jr, j+1 ) = t( jr, j+1 ) - temp*t2
                  t( jr, j+2 ) = t( jr, j+2 ) - temp*t3
  270          CONTINUE
               IF( ilz ) THEN
                  DO 280 jr = 1, n
                     temp = z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
     $                      z( jr, j+2 )
                     z( jr, j ) = z( jr, j ) - temp*tau
                     z( jr, j+1 ) = z( jr, j+1 ) - temp*t2
                     z( jr, j+2 ) = z( jr, j+2 ) - temp*t3
  280             CONTINUE
               END IF
               t( j+1, j ) = zero
               t( j+2, j ) = zero
  290       CONTINUE
*
*           Last elements: Use Givens rotations
*
*           Rotations from the left
*
            j = ilast - 1
            temp = h( j, j-1 )
            CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
            h( j+1, j-1 ) = zero
*
            DO 300 jc = j, ilastm
               temp = c*h( j, jc ) + s*h( j+1, jc )
               h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
               h( j, jc ) = temp
               temp2 = c*t( j, jc ) + s*t( j+1, jc )
               t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
               t( j, jc ) = temp2
  300       CONTINUE
            IF( ilq ) THEN
               DO 310 jr = 1, n
                  temp = c*q( jr, j ) + s*q( jr, j+1 )
                  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
                  q( jr, j ) = temp
  310          CONTINUE
            END IF
*
*           Rotations from the right.
*
            temp = t( j+1, j+1 )
            CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
            t( j+1, j ) = zero
*
            DO 320 jr = ifrstm, ilast
               temp = c*h( jr, j+1 ) + s*h( jr, j )
               h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
               h( jr, j+1 ) = temp
  320       CONTINUE
            DO 330 jr = ifrstm, ilast - 1
               temp = c*t( jr, j+1 ) + s*t( jr, j )
               t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
               t( jr, j+1 ) = temp
  330       CONTINUE
            IF( ilz ) THEN
               DO 340 jr = 1, n
                  temp = c*z( jr, j+1 ) + s*z( jr, j )
                  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
                  z( jr, j+1 ) = temp
  340          CONTINUE
            END IF
*
*           End of Double-Shift code
*
         END IF
*
         GO TO 350
*
*        End of iteration loop
*
  350    CONTINUE
  360 CONTINUE
*
*     Drop-through = non-convergence
*
      info = ilast
      GO TO 420
*
*     Successful completion of all QZ steps
*
  380 CONTINUE
*
*     Set Eigenvalues 1:ILO-1
*
      DO 410 j = 1, ilo - 1
         IF( t( j, j ).LT.zero ) THEN
            IF( ilschr ) THEN
               DO 390 jr = 1, j
                  h( jr, j ) = -h( jr, j )
                  t( jr, j ) = -t( jr, j )
  390          CONTINUE
            ELSE
               h( j, j ) = -h( j, j )
               t( j, j ) = -t( j, j )
            END IF
            IF( ilz ) THEN
               DO 400 jr = 1, n
                  z( jr, j ) = -z( jr, j )
  400          CONTINUE
            END IF
         END IF
         alphar( j ) = h( j, j )
         alphai( j ) = zero
         beta( j ) = t( j, j )
  410 CONTINUE
*
*     Normal Termination
*
      info = 0
*
*     Exit (other than argument error) -- return optimal workspace size
*
  420 CONTINUE
      work( 1 ) = dble( n )
      RETURN
*
*     End of DHGEQZ
*
      END