chgeqz.f

*> \brief \b CHGEQZ
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CHGEQZ + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chgeqz.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chgeqz.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chgeqz.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
*                          ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
*                          RWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          COMPQ, COMPZ, JOB
*       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
*       ..
*       .. Array Arguments ..
*       REAL               RWORK( * )
*       COMPLEX            ALPHA( * ), BETA( * ), H( LDH, * ),
*      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
*      $                   Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
*> where H is an upper Hessenberg matrix and T is upper triangular,
*> using the single-shift QZ method.
*> Matrix pairs of this type are produced by the reduction to
*> generalized upper Hessenberg form of a complex matrix pair (A,B):
*>
*>    A = Q1*H*Z1**H,  B = Q1*T*Z1**H,
*>
*> as computed by CGGHRD.
*>
*> If JOB='S', then the Hessenberg-triangular pair (H,T) is
*> also reduced to generalized Schur form,
*>
*>    H = Q*S*Z**H,  T = Q*P*Z**H,
*>
*> where Q and Z are unitary matrices and S and P are upper triangular.
*>
*> Optionally, the unitary matrix Q from the generalized Schur
*> factorization may be postmultiplied into an input matrix Q1, and the
*> unitary matrix Z may be postmultiplied into an input matrix Z1.
*> If Q1 and Z1 are the unitary matrices from CGGHRD that reduced
*> the matrix pair (A,B) to generalized Hessenberg form, then the output
*> matrices Q1*Q and Z1*Z are the unitary factors from the generalized
*> Schur factorization of (A,B):
*>
*>    A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.
*>
*> To avoid overflow, eigenvalues of the matrix pair (H,T)
*> (equivalently, of (A,B)) are computed as a pair of complex values
*> (alpha,beta).  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.
*> The values of alpha and beta for the i-th eigenvalue 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': Computer 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 a unitary 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': Q is initialized to the unit matrix and the matrix Z
*>                 of right Schur vectors of (H,T) is returned;
*>          = 'V': Z must contain a unitary 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 COMPLEX array, dimension (LDH, N)
*>          On entry, the N-by-N upper Hessenberg matrix H.
*>          On exit, if JOB = 'S', H contains the upper triangular
*>          matrix S from the generalized Schur factorization.
*>          If JOB = 'E', the diagonal of H matches that 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 COMPLEX 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.
*>          If JOB = 'E', the diagonal of T matches that 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] ALPHA
*> \verbatim
*>          ALPHA is COMPLEX array, dimension (N)
*>          The complex scalars alpha that define the eigenvalues of
*>          GNEP.  ALPHA(i) = S(i,i) in the generalized Schur
*>          factorization.
*> \endverbatim
*>
*> \param[out] BETA
*> \verbatim
*>          BETA is COMPLEX array, dimension (N)
*>          The real non-negative scalars beta that define the
*>          eigenvalues of GNEP.  BETA(i) = P(i,i) in the generalized
*>          Schur factorization.
*>
*>          Together, the quantities alpha = ALPHA(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 COMPLEX array, dimension (LDQ, N)
*>          On entry, if COMPQ = 'V', the unitary matrix Q1 used in the
*>          reduction of (A,B) to generalized Hessenberg form.
*>          On exit, if COMPQ = 'I', the unitary matrix of left Schur
*>          vectors of (H,T), and if COMPQ = 'V', the unitary 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 COMPLEX array, dimension (LDZ, N)
*>          On entry, if COMPZ = 'V', the unitary matrix Z1 used in the
*>          reduction of (A,B) to generalized Hessenberg form.
*>          On exit, if COMPZ = 'I', the unitary matrix of right Schur
*>          vectors of (H,T), and if COMPZ = 'V', the unitary 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 COMPLEX 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] RWORK
*> \verbatim
*>          RWORK is REAL array, dimension (N)
*> \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 ALPHA(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 ALPHA(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
*>
*>  We assume that complex ABS works as long as its value is less than
*>  overflow.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE chgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
     $                   ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
     $                   RWORK, 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 ..
      REAL               RWORK( * )
      COMPLEX            ALPHA( * ), BETA( * ), H( ldh, * ),
     $                   q( ldq, * ), t( ldt, * ), work( * ),
     $                   z( ldz, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX            CZERO, CONE
      parameter( czero = ( 0.0e+0, 0.0e+0 ),
     $                   cone = ( 1.0e+0, 0.0e+0 ) )
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
      REAL               HALF
      parameter( half = 0.5e+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
      INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
     $                   ilastm, in, ischur, istart, j, jc, jch, jiter,
     $                   jr, maxit
      REAL               ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
     $                   c, safmin, temp, temp2, tempr, ulp
      COMPLEX            ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
     $                   ctemp3, eshift, s, shift, signbc,
     $                   u12, x, abi12, y
*     ..
*     .. External Functions ..
      COMPLEX            CLADIV
      LOGICAL            LSAME
      REAL               CLANHS, SLAMCH
      EXTERNAL           cladiv, lsame, clanhs, slamch
*     ..
*     .. External Subroutines ..
      EXTERNAL           clartg, claset, crot, cscal, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, aimag, cmplx, conjg, max, min, REAL, SQRT
*     ..
*     .. Statement Functions ..
      REAL               ABS1
*     ..
*     .. Statement Function definitions ..
      abs1( x ) = abs( REAL( X ) ) + abs( AIMAG( x ) )
*     ..
*     .. 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
         ilschr = .true.
         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
         ilq = .true.
         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
         ilz = .true.
         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 = -14
      ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
         info = -16
      ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
         info = -18
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CHGEQZ', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
*     WORK( 1 ) = CMPLX( 1 )
      IF( n.LE.0 ) THEN
         work( 1 ) = cmplx( 1 )
         RETURN
      END IF
*
*     Initialize Q and Z
*
      IF( icompq.EQ.3 )
     $   CALL claset( 'Full', n, n, czero, cone, q, ldq )
      IF( icompz.EQ.3 )
     $   CALL claset( 'Full', n, n, czero, cone, z, ldz )
*
*     Machine Constants
*
      in = ihi + 1 - ilo
      safmin = slamch( 'S' )
      ulp = slamch( 'E' )*slamch( 'B' )
      anorm = clanhs( 'F', in, h( ilo, ilo ), ldh, rwork )
      bnorm = clanhs( 'F', in, t( ilo, ilo ), ldt, rwork )
      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 10 j = ihi + 1, n
         absb = abs( t( j, j ) )
         IF( absb.GT.safmin ) THEN
            signbc = conjg( t( j, j ) / absb )
            t( j, j ) = absb
            IF( ilschr ) THEN
               CALL cscal( j-1, signbc, t( 1, j ), 1 )
               CALL cscal( j, signbc, h( 1, j ), 1 )
            ELSE
               CALL cscal( 1, signbc, h( j, j ), 1 )
            END IF
            IF( ilz )
     $         CALL cscal( n, signbc, z( 1, j ), 1 )
         ELSE
            t( j, j ) = czero
         END IF
         alpha( j ) = h( j, j )
         beta( j ) = t( j, j )
   10 CONTINUE
*
*     If IHI < ILO, skip QZ steps
*
      IF( ihi.LT.ilo )
     $   GO TO 190
*
*     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 = czero
      maxit = 30*( ihi-ilo+1 )
*
      DO 170 jiter = 1, maxit
*
*        Check for too many iterations.
*
         IF( jiter.GT.maxit )
     $      GO TO 180
*
*        Split the matrix if possible.
*
*        Two tests:
*           1: H(j,j-1)=0  or  j=ILO
*           2: T(j,j)=0
*
*        Special case: j=ILAST
*
         IF( ilast.EQ.ilo ) THEN
            GO TO 60
         ELSE
            IF( abs1( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*( 
     $         abs1( h( ilast, ilast ) ) + abs1( h( ilast-1, ilast-1 ) 
     $         ) ) ) ) THEN
               h( ilast, ilast-1 ) = czero
               GO TO 60
            END IF
         END IF
*
         IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
            t( ilast, ilast ) = czero
            GO TO 50
         END IF
*
*        General case: j<ILAST
*
         DO 40 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( abs1( h( j, j-1 ) ).LE.max( safmin, ulp*( 
     $            abs1( h( j, j ) ) + abs1( h( j-1, j-1 ) ) 
     $            ) ) ) THEN
                  h( j, j-1 ) = czero
                  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 ) = czero
*
*              Test 1a: Check for 2 consecutive small subdiagonals in A
*
               ilazr2 = .false.
               IF( .NOT.ilazro ) THEN
                  IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
     $                j ) ) ).LE.abs1( h( j, j ) )*( 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 20 jch = j, ilast - 1
                     ctemp = h( jch, jch )
                     CALL clartg( ctemp, h( jch+1, jch ), c, s,
     $                            h( jch, jch ) )
                     h( jch+1, jch ) = czero
                     CALL crot( ilastm-jch, h( jch, jch+1 ), ldh,
     $                          h( jch+1, jch+1 ), ldh, c, s )
                     CALL crot( ilastm-jch, t( jch, jch+1 ), ldt,
     $                          t( jch+1, jch+1 ), ldt, c, s )
                     IF( ilq )
     $                  CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
     $                             c, conjg( s ) )
                     IF( ilazr2 )
     $                  h( jch, jch-1 ) = h( jch, jch-1 )*c
                     ilazr2 = .false.
                     IF( abs1( t( jch+1, jch+1 ) ).GE.btol ) THEN
                        IF( jch+1.GE.ilast ) THEN
                           GO TO 60
                        ELSE
                           ifirst = jch + 1
                           GO TO 70
                        END IF
                     END IF
                     t( jch+1, jch+1 ) = czero
   20             CONTINUE
                  GO TO 50
               ELSE
*
*                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
*                 Then process as in the case T(ILAST,ILAST)=0
*
                  DO 30 jch = j, ilast - 1
                     ctemp = t( jch, jch+1 )
                     CALL clartg( ctemp, t( jch+1, jch+1 ), c, s,
     $                            t( jch, jch+1 ) )
                     t( jch+1, jch+1 ) = czero
                     IF( jch.LT.ilastm-1 )
     $                  CALL crot( ilastm-jch-1, t( jch, jch+2 ), ldt,
     $                             t( jch+1, jch+2 ), ldt, c, s )
                     CALL crot( ilastm-jch+2, h( jch, jch-1 ), ldh,
     $                          h( jch+1, jch-1 ), ldh, c, s )
                     IF( ilq )
     $                  CALL crot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
     $                             c, conjg( s ) )
                     ctemp = h( jch+1, jch )
                     CALL clartg( ctemp, h( jch+1, jch-1 ), c, s,
     $                            h( jch+1, jch ) )
                     h( jch+1, jch-1 ) = czero
                     CALL crot( jch+1-ifrstm, h( ifrstm, jch ), 1,
     $                          h( ifrstm, jch-1 ), 1, c, s )
                     CALL crot( jch-ifrstm, t( ifrstm, jch ), 1,
     $                          t( ifrstm, jch-1 ), 1, c, s )
                     IF( ilz )
     $                  CALL crot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
     $                             c, s )
   30             CONTINUE
                  GO TO 50
               END IF
            ELSE IF( ilazro ) THEN
*
*              Only test 1 passed -- work on J:ILAST
*
               ifirst = j
               GO TO 70
            END IF
*
*           Neither test passed -- try next J
*
   40    CONTINUE
*
*        (Drop-through is "impossible")
*
         info = 2*n + 1
         GO TO 210
*
*        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
*        1x1 block.
*
   50    CONTINUE
         ctemp = h( ilast, ilast )
         CALL clartg( ctemp, h( ilast, ilast-1 ), c, s,
     $                h( ilast, ilast ) )
         h( ilast, ilast-1 ) = czero
         CALL crot( ilast-ifrstm, h( ifrstm, ilast ), 1,
     $              h( ifrstm, ilast-1 ), 1, c, s )
         CALL crot( ilast-ifrstm, t( ifrstm, ilast ), 1,
     $              t( ifrstm, ilast-1 ), 1, c, s )
         IF( ilz )
     $      CALL crot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
*
*        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
*
   60    CONTINUE
         absb = abs( t( ilast, ilast ) )
         IF( absb.GT.safmin ) THEN
            signbc = conjg( t( ilast, ilast ) / absb )
            t( ilast, ilast ) = absb
            IF( ilschr ) THEN
               CALL cscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
               CALL cscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
     $                     1 )
            ELSE
               CALL cscal( 1, signbc, h( ilast, ilast ), 1 )
            END IF
            IF( ilz )
     $         CALL cscal( n, signbc, z( 1, ilast ), 1 )
         ELSE
            t( ilast, ilast ) = czero
         END IF
         alpha( ilast ) = h( ilast, ilast )
         beta( ilast ) = t( ilast, ilast )
*
*        Go to next block -- exit if finished.
*
         ilast = ilast - 1
         IF( ilast.LT.ilo )
     $      GO TO 190
*
*        Reset counters
*
         iiter = 0
         eshift = czero
         IF( .NOT.ilschr ) THEN
            ilastm = ilast
            IF( ifrstm.GT.ilast )
     $         ifrstm = ilo
         END IF
         GO TO 160
*
*        QZ step
*
*        This iteration only involves rows/columns IFIRST:ILAST.  We
*        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
*
   70    CONTINUE
         iiter = iiter + 1
         IF( .NOT.ilschr ) THEN
            ifrstm = ifirst
         END IF
*
*        Compute the Shift.
*
*        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.NE.iiter ) THEN
*
*           The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
*           the bottom-right 2x2 block of A inv(B) which is nearest to
*           the bottom-right element.
*
*           We factor B as U*D, where U has unit diagonals, and
*           compute (A*inv(D))*inv(U).
*
            u12 = ( bscale*t( ilast-1, ilast ) ) /
     $            ( bscale*t( ilast, ilast ) )
            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 ) )
            abi22 = ad22 - u12*ad21
            abi12 = ad12 - u12*ad11
*
            shift = abi22
            ctemp = sqrt( abi12 )*sqrt( ad21 )
            temp = abs1( ctemp )
            IF( ctemp.NE.zero ) THEN
               x = half*( ad11-shift )
               temp2 = abs1( x )
               temp = max( temp, abs1( x ) )
               y = temp*sqrt( ( x / temp )**2+( ctemp / temp )**2 )
               IF( temp2.GT.zero ) THEN
                  IF( REAL( x / temp2 )*REAL( y )+
     $                aimag( x / temp2 )*aimag( y ).LT.zero )y = -y
               END IF
               shift = shift - ctemp*cladiv( ctemp, ( x+y ) )
            END IF
         ELSE
*
*           Exceptional shift.  Chosen for no particularly good reason.
*
            IF( ( iiter / 20 )*20.EQ.iiter .AND. 
     $         bscale*abs1(t( ilast, ilast )).GT.safmin ) THEN
               eshift = eshift + ( ascale*h( ilast,
     $            ilast ) )/( bscale*t( ilast, ilast ) )
            ELSE
               eshift = eshift + ( ascale*h( ilast,
     $            ilast-1 ) )/( bscale*t( ilast-1, ilast-1 ) )
            END IF
            shift = eshift
         END IF
*
*        Now check for two consecutive small subdiagonals.
*
         DO 80 j = ilast - 1, ifirst + 1, -1
            istart = j
            ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
            temp = abs1( ctemp )
            temp2 = ascale*abs1( h( j+1, j ) )
            tempr = max( temp, temp2 )
            IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
               temp = temp / tempr
               temp2 = temp2 / tempr
            END IF
            IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
     $         GO TO 90
   80    CONTINUE
*
         istart = ifirst
         ctemp = ascale*h( ifirst, ifirst ) -
     $           shift*( bscale*t( ifirst, ifirst ) )
   90    CONTINUE
*
*        Do an implicit-shift QZ sweep.
*
*        Initial Q
*
         ctemp2 = ascale*h( istart+1, istart )
         CALL clartg( ctemp, ctemp2, c, s, ctemp3 )
*
*        Sweep
*
         DO 150 j = istart, ilast - 1
            IF( j.GT.istart ) THEN
               ctemp = h( j, j-1 )
               CALL clartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
               h( j+1, j-1 ) = czero
            END IF
*
            DO 100 jc = j, ilastm
               ctemp = c*h( j, jc ) + s*h( j+1, jc )
               h( j+1, jc ) = -conjg( s )*h( j, jc ) + c*h( j+1, jc )
               h( j, jc ) = ctemp
               ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
               t( j+1, jc ) = -conjg( s )*t( j, jc ) + c*t( j+1, jc )
               t( j, jc ) = ctemp2
  100       CONTINUE
            IF( ilq ) THEN
               DO 110 jr = 1, n
                  ctemp = c*q( jr, j ) + conjg( s )*q( jr, j+1 )
                  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
                  q( jr, j ) = ctemp
  110          CONTINUE
            END IF
*
            ctemp = t( j+1, j+1 )
            CALL clartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
            t( j+1, j ) = czero
*
            DO 120 jr = ifrstm, min( j+2, ilast )
               ctemp = c*h( jr, j+1 ) + s*h( jr, j )
               h( jr, j ) = -conjg( s )*h( jr, j+1 ) + c*h( jr, j )
               h( jr, j+1 ) = ctemp
  120       CONTINUE
            DO 130 jr = ifrstm, j
               ctemp = c*t( jr, j+1 ) + s*t( jr, j )
               t( jr, j ) = -conjg( s )*t( jr, j+1 ) + c*t( jr, j )
               t( jr, j+1 ) = ctemp
  130       CONTINUE
            IF( ilz ) THEN
               DO 140 jr = 1, n
                  ctemp = c*z( jr, j+1 ) + s*z( jr, j )
                  z( jr, j ) = -conjg( s )*z( jr, j+1 ) + c*z( jr, j )
                  z( jr, j+1 ) = ctemp
  140          CONTINUE
            END IF
  150    CONTINUE
*
  160    CONTINUE
*
  170 CONTINUE
*
*     Drop-through = non-convergence
*
  180 CONTINUE
      info = ilast
      GO TO 210
*
*     Successful completion of all QZ steps
*
  190 CONTINUE
*
*     Set Eigenvalues 1:ILO-1
*
      DO 200 j = 1, ilo - 1
         absb = abs( t( j, j ) )
         IF( absb.GT.safmin ) THEN
            signbc = conjg( t( j, j ) / absb )
            t( j, j ) = absb
            IF( ilschr ) THEN
               CALL cscal( j-1, signbc, t( 1, j ), 1 )
               CALL cscal( j, signbc, h( 1, j ), 1 )
            ELSE
               CALL cscal( 1, signbc, h( j, j ), 1 )
            END IF
            IF( ilz )
     $         CALL cscal( n, signbc, z( 1, j ), 1 )
         ELSE
            t( j, j ) = czero
         END IF
         alpha( j ) = h( j, j )
         beta( j ) = t( j, j )
  200 CONTINUE
*
*     Normal Termination
*
      info = 0
*
*     Exit (other than argument error) -- return optimal workspace size
*
  210 CONTINUE
      work( 1 ) = cmplx( n )
      RETURN
*
*     End of CHGEQZ
*
      END