SRC/chgeqz.f

   1 *> \brief \b CHGEQZ
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CHGEQZ + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chgeqz.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chgeqz.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chgeqz.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
  22 *                          ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
  23 *                          RWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          COMPQ, COMPZ, JOB
  27 *       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       REAL               RWORK( * )
  31 *       COMPLEX            ALPHA( * ), BETA( * ), H( LDH, * ),
  32 *      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
  33 *      $                   Z( LDZ, * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> CHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
  43 *> where H is an upper Hessenberg matrix and T is upper triangular,
  44 *> using the single-shift QZ method.
  45 *> Matrix pairs of this type are produced by the reduction to
  46 *> generalized upper Hessenberg form of a complex matrix pair (A,B):
  47 *>
  48 *>    A = Q1*H*Z1**H,  B = Q1*T*Z1**H,
  49 *>
  50 *> as computed by CGGHRD.
  51 *>
  52 *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
  53 *> also reduced to generalized Schur form,
  54 *>
  55 *>    H = Q*S*Z**H,  T = Q*P*Z**H,
  56 *>
  57 *> where Q and Z are unitary matrices and S and P are upper triangular.
  58 *>
  59 *> Optionally, the unitary matrix Q from the generalized Schur
  60 *> factorization may be postmultiplied into an input matrix Q1, and the
  61 *> unitary matrix Z may be postmultiplied into an input matrix Z1.
  62 *> If Q1 and Z1 are the unitary matrices from CGGHRD that reduced
  63 *> the matrix pair (A,B) to generalized Hessenberg form, then the output
  64 *> matrices Q1*Q and Z1*Z are the unitary factors from the generalized
  65 *> Schur factorization of (A,B):
  66 *>
  67 *>    A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.
  68 *>
  69 *> To avoid overflow, eigenvalues of the matrix pair (H,T)
  70 *> (equivalently, of (A,B)) are computed as a pair of complex values
  71 *> (alpha,beta).  If beta is nonzero, lambda = alpha / beta is an
  72 *> eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP)
  73 *>    A*x = lambda*B*x
  74 *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
  75 *> alternate form of the GNEP
  76 *>    mu*A*y = B*y.
  77 *> The values of alpha and beta for the i-th eigenvalue can be read
  78 *> directly from the generalized Schur form:  alpha = S(i,i),
  79 *> beta = P(i,i).
  80 *>
  81 *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
  82 *>      Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
  83 *>      pp. 241--256.
  84 *> \endverbatim
  85 *
  86 *  Arguments:
  87 *  ==========
  88 *
  89 *> \param[in] JOB
  90 *> \verbatim
  91 *>          JOB is CHARACTER*1
  92 *>          = 'E': Compute eigenvalues only;
  93 *>          = 'S': Computer eigenvalues and the Schur form.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] COMPQ
  97 *> \verbatim
  98 *>          COMPQ is CHARACTER*1
  99 *>          = 'N': Left Schur vectors (Q) are not computed;
 100 *>          = 'I': Q is initialized to the unit matrix and the matrix Q
 101 *>                 of left Schur vectors of (H,T) is returned;
 102 *>          = 'V': Q must contain a unitary matrix Q1 on entry and
 103 *>                 the product Q1*Q is returned.
 104 *> \endverbatim
 105 *>
 106 *> \param[in] COMPZ
 107 *> \verbatim
 108 *>          COMPZ is CHARACTER*1
 109 *>          = 'N': Right Schur vectors (Z) are not computed;
 110 *>          = 'I': Q is initialized to the unit matrix and the matrix Z
 111 *>                 of right Schur vectors of (H,T) is returned;
 112 *>          = 'V': Z must contain a unitary matrix Z1 on entry and
 113 *>                 the product Z1*Z is returned.
 114 *> \endverbatim
 115 *>
 116 *> \param[in] N
 117 *> \verbatim
 118 *>          N is INTEGER
 119 *>          The order of the matrices H, T, Q, and Z.  N >= 0.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] ILO
 123 *> \verbatim
 124 *>          ILO is INTEGER
 125 *> \endverbatim
 126 *>
 127 *> \param[in] IHI
 128 *> \verbatim
 129 *>          IHI is INTEGER
 130 *>          ILO and IHI mark the rows and columns of H which are in
 131 *>          Hessenberg form.  It is assumed that A is already upper
 132 *>          triangular in rows and columns 1:ILO-1 and IHI+1:N.
 133 *>          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
 134 *> \endverbatim
 135 *>
 136 *> \param[in,out] H
 137 *> \verbatim
 138 *>          H is COMPLEX array, dimension (LDH, N)
 139 *>          On entry, the N-by-N upper Hessenberg matrix H.
 140 *>          On exit, if JOB = 'S', H contains the upper triangular
 141 *>          matrix S from the generalized Schur factorization.
 142 *>          If JOB = 'E', the diagonal of H matches that of S, but
 143 *>          the rest of H is unspecified.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] LDH
 147 *> \verbatim
 148 *>          LDH is INTEGER
 149 *>          The leading dimension of the array H.  LDH >= max( 1, N ).
 150 *> \endverbatim
 151 *>
 152 *> \param[in,out] T
 153 *> \verbatim
 154 *>          T is COMPLEX array, dimension (LDT, N)
 155 *>          On entry, the N-by-N upper triangular matrix T.
 156 *>          On exit, if JOB = 'S', T contains the upper triangular
 157 *>          matrix P from the generalized Schur factorization.
 158 *>          If JOB = 'E', the diagonal of T matches that of P, but
 159 *>          the rest of T is unspecified.
 160 *> \endverbatim
 161 *>
 162 *> \param[in] LDT
 163 *> \verbatim
 164 *>          LDT is INTEGER
 165 *>          The leading dimension of the array T.  LDT >= max( 1, N ).
 166 *> \endverbatim
 167 *>
 168 *> \param[out] ALPHA
 169 *> \verbatim
 170 *>          ALPHA is COMPLEX array, dimension (N)
 171 *>          The complex scalars alpha that define the eigenvalues of
 172 *>          GNEP.  ALPHA(i) = S(i,i) in the generalized Schur
 173 *>          factorization.
 174 *> \endverbatim
 175 *>
 176 *> \param[out] BETA
 177 *> \verbatim
 178 *>          BETA is COMPLEX array, dimension (N)
 179 *>          The real non-negative scalars beta that define the
 180 *>          eigenvalues of GNEP.  BETA(i) = P(i,i) in the generalized
 181 *>          Schur factorization.
 182 *>
 183 *>          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
 184 *>          represent the j-th eigenvalue of the matrix pair (A,B), in
 185 *>          one of the forms lambda = alpha/beta or mu = beta/alpha.
 186 *>          Since either lambda or mu may overflow, they should not,
 187 *>          in general, be computed.
 188 *> \endverbatim
 189 *>
 190 *> \param[in,out] Q
 191 *> \verbatim
 192 *>          Q is COMPLEX array, dimension (LDQ, N)
 193 *>          On entry, if COMPQ = 'V', the unitary matrix Q1 used in the
 194 *>          reduction of (A,B) to generalized Hessenberg form.
 195 *>          On exit, if COMPQ = 'I', the unitary matrix of left Schur
 196 *>          vectors of (H,T), and if COMPQ = 'V', the unitary matrix of
 197 *>          left Schur vectors of (A,B).
 198 *>          Not referenced if COMPQ = 'N'.
 199 *> \endverbatim
 200 *>
 201 *> \param[in] LDQ
 202 *> \verbatim
 203 *>          LDQ is INTEGER
 204 *>          The leading dimension of the array Q.  LDQ >= 1.
 205 *>          If COMPQ='V' or 'I', then LDQ >= N.
 206 *> \endverbatim
 207 *>
 208 *> \param[in,out] Z
 209 *> \verbatim
 210 *>          Z is COMPLEX array, dimension (LDZ, N)
 211 *>          On entry, if COMPZ = 'V', the unitary matrix Z1 used in the
 212 *>          reduction of (A,B) to generalized Hessenberg form.
 213 *>          On exit, if COMPZ = 'I', the unitary matrix of right Schur
 214 *>          vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
 215 *>          right Schur vectors of (A,B).
 216 *>          Not referenced if COMPZ = 'N'.
 217 *> \endverbatim
 218 *>
 219 *> \param[in] LDZ
 220 *> \verbatim
 221 *>          LDZ is INTEGER
 222 *>          The leading dimension of the array Z.  LDZ >= 1.
 223 *>          If COMPZ='V' or 'I', then LDZ >= N.
 224 *> \endverbatim
 225 *>
 226 *> \param[out] WORK
 227 *> \verbatim
 228 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 229 *>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
 230 *> \endverbatim
 231 *>
 232 *> \param[in] LWORK
 233 *> \verbatim
 234 *>          LWORK is INTEGER
 235 *>          The dimension of the array WORK.  LWORK >= max(1,N).
 236 *>
 237 *>          If LWORK = -1, then a workspace query is assumed; the routine
 238 *>          only calculates the optimal size of the WORK array, returns
 239 *>          this value as the first entry of the WORK array, and no error
 240 *>          message related to LWORK is issued by XERBLA.
 241 *> \endverbatim
 242 *>
 243 *> \param[out] RWORK
 244 *> \verbatim
 245 *>          RWORK is REAL array, dimension (N)
 246 *> \endverbatim
 247 *>
 248 *> \param[out] INFO
 249 *> \verbatim
 250 *>          INFO is INTEGER
 251 *>          = 0: successful exit
 252 *>          < 0: if INFO = -i, the i-th argument had an illegal value
 253 *>          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
 254 *>                     in Schur form, but ALPHA(i) and BETA(i),
 255 *>                     i=INFO+1,...,N should be correct.
 256 *>          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
 257 *>                     in Schur form, but ALPHA(i) and BETA(i),
 258 *>                     i=INFO-N+1,...,N should be correct.
 259 *> \endverbatim
 260 *
 261 *  Authors:
 262 *  ========
 263 *
 264 *> \author Univ. of Tennessee
 265 *> \author Univ. of California Berkeley
 266 *> \author Univ. of Colorado Denver
 267 *> \author NAG Ltd.
 268 *
 269 *> \date April 2012
 270 *
 271 *> \ingroup complexGEcomputational
 272 *
 273 *> \par Further Details:
 274 *  =====================
 275 *>
 276 *> \verbatim
 277 *>
 278 *>  We assume that complex ABS works as long as its value is less than
 279 *>  overflow.
 280 *> \endverbatim
 281 *>
 282 *  =====================================================================
 283       SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
 284      $                   ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
 285      $                   RWORK, INFO )
 286 *
 287 *  -- LAPACK computational routine (version 3.6.1) --
 288 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 289 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 290 *     April 2012
 291 *
 292 *     .. Scalar Arguments ..
 293       CHARACTER          COMPQ, COMPZ, JOB
 294       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
 295 *     ..
 296 *     .. Array Arguments ..
 297       REAL               RWORK( * )
 298       COMPLEX            ALPHA( * ), BETA( * ), H( LDH, * ),
 299      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
 300      $                   Z( LDZ, * )
 301 *     ..
 302 *
 303 *  =====================================================================
 304 *
 305 *     .. Parameters ..
 306       COMPLEX            CZERO, CONE
 307       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
 308      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
 309       REAL               ZERO, ONE
 310       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 311       REAL               HALF
 312       PARAMETER          ( HALF = 0.5E+0 )
 313 *     ..
 314 *     .. Local Scalars ..
 315       LOGICAL            ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
 316       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
 317      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
 318      $                   JR, MAXIT
 319       REAL               ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
 320      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
 321       COMPLEX            ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
 322      $                   CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
 323      $                   U12, X
 324 *     ..
 325 *     .. External Functions ..
 326       LOGICAL            LSAME
 327       REAL               CLANHS, SLAMCH
 328       EXTERNAL           LSAME, CLANHS, SLAMCH
 329 *     ..
 330 *     .. External Subroutines ..
 331       EXTERNAL           CLARTG, CLASET, CROT, CSCAL, XERBLA
 332 *     ..
 333 *     .. Intrinsic Functions ..
 334       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT
 335 *     ..
 336 *     .. Statement Functions ..
 337       REAL               ABS1
 338 *     ..
 339 *     .. Statement Function definitions ..
 340       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
 341 *     ..
 342 *     .. Executable Statements ..
 343 *
 344 *     Decode JOB, COMPQ, COMPZ
 345 *
 346       IF( LSAME( JOB, 'E' ) ) THEN
 347          ILSCHR = .FALSE.
 348          ISCHUR = 1
 349       ELSE IF( LSAME( JOB, 'S' ) ) THEN
 350          ILSCHR = .TRUE.
 351          ISCHUR = 2
 352       ELSE
 353          ISCHUR = 0
 354       END IF
 355 *
 356       IF( LSAME( COMPQ, 'N' ) ) THEN
 357          ILQ = .FALSE.
 358          ICOMPQ = 1
 359       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
 360          ILQ = .TRUE.
 361          ICOMPQ = 2
 362       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
 363          ILQ = .TRUE.
 364          ICOMPQ = 3
 365       ELSE
 366          ICOMPQ = 0
 367       END IF
 368 *
 369       IF( LSAME( COMPZ, 'N' ) ) THEN
 370          ILZ = .FALSE.
 371          ICOMPZ = 1
 372       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 373          ILZ = .TRUE.
 374          ICOMPZ = 2
 375       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 376          ILZ = .TRUE.
 377          ICOMPZ = 3
 378       ELSE
 379          ICOMPZ = 0
 380       END IF
 381 *
 382 *     Check Argument Values
 383 *
 384       INFO = 0
 385       WORK( 1 ) = MAX( 1, N )
 386       LQUERY = ( LWORK.EQ.-1 )
 387       IF( ISCHUR.EQ.0 ) THEN
 388          INFO = -1
 389       ELSE IF( ICOMPQ.EQ.0 ) THEN
 390          INFO = -2
 391       ELSE IF( ICOMPZ.EQ.0 ) THEN
 392          INFO = -3
 393       ELSE IF( N.LT.0 ) THEN
 394          INFO = -4
 395       ELSE IF( ILO.LT.1 ) THEN
 396          INFO = -5
 397       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
 398          INFO = -6
 399       ELSE IF( LDH.LT.N ) THEN
 400          INFO = -8
 401       ELSE IF( LDT.LT.N ) THEN
 402          INFO = -10
 403       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
 404          INFO = -14
 405       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
 406          INFO = -16
 407       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 408          INFO = -18
 409       END IF
 410       IF( INFO.NE.0 ) THEN
 411          CALL XERBLA( 'CHGEQZ', -INFO )
 412          RETURN
 413       ELSE IF( LQUERY ) THEN
 414          RETURN
 415       END IF
 416 *
 417 *     Quick return if possible
 418 *
 419 *     WORK( 1 ) = CMPLX( 1 )
 420       IF( N.LE.0 ) THEN
 421          WORK( 1 ) = CMPLX( 1 )
 422          RETURN
 423       END IF
 424 *
 425 *     Initialize Q and Z
 426 *
 427       IF( ICOMPQ.EQ.3 )
 428      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
 429       IF( ICOMPZ.EQ.3 )
 430      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
 431 *
 432 *     Machine Constants
 433 *
 434       IN = IHI + 1 - ILO
 435       SAFMIN = SLAMCH( 'S' )
 436       ULP = SLAMCH( 'E' )*SLAMCH( 'B' )
 437       ANORM = CLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK )
 438       BNORM = CLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK )
 439       ATOL = MAX( SAFMIN, ULP*ANORM )
 440       BTOL = MAX( SAFMIN, ULP*BNORM )
 441       ASCALE = ONE / MAX( SAFMIN, ANORM )
 442       BSCALE = ONE / MAX( SAFMIN, BNORM )
 443 *
 444 *
 445 *     Set Eigenvalues IHI+1:N
 446 *
 447       DO 10 J = IHI + 1, N
 448          ABSB = ABS( T( J, J ) )
 449          IF( ABSB.GT.SAFMIN ) THEN
 450             SIGNBC = CONJG( T( J, J ) / ABSB )
 451             T( J, J ) = ABSB
 452             IF( ILSCHR ) THEN
 453                CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
 454                CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
 455             ELSE
 456                CALL CSCAL( 1, SIGNBC, H( J, J ), 1 )
 457             END IF
 458             IF( ILZ )
 459      $         CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
 460          ELSE
 461             T( J, J ) = CZERO
 462          END IF
 463          ALPHA( J ) = H( J, J )
 464          BETA( J ) = T( J, J )
 465    10 CONTINUE
 466 *
 467 *     If IHI < ILO, skip QZ steps
 468 *
 469       IF( IHI.LT.ILO )
 470      $   GO TO 190
 471 *
 472 *     MAIN QZ ITERATION LOOP
 473 *
 474 *     Initialize dynamic indices
 475 *
 476 *     Eigenvalues ILAST+1:N have been found.
 477 *        Column operations modify rows IFRSTM:whatever
 478 *        Row operations modify columns whatever:ILASTM
 479 *
 480 *     If only eigenvalues are being computed, then
 481 *        IFRSTM is the row of the last splitting row above row ILAST;
 482 *        this is always at least ILO.
 483 *     IITER counts iterations since the last eigenvalue was found,
 484 *        to tell when to use an extraordinary shift.
 485 *     MAXIT is the maximum number of QZ sweeps allowed.
 486 *
 487       ILAST = IHI
 488       IF( ILSCHR ) THEN
 489          IFRSTM = 1
 490          ILASTM = N
 491       ELSE
 492          IFRSTM = ILO
 493          ILASTM = IHI
 494       END IF
 495       IITER = 0
 496       ESHIFT = CZERO
 497       MAXIT = 30*( IHI-ILO+1 )
 498 *
 499       DO 170 JITER = 1, MAXIT
 500 *
 501 *        Check for too many iterations.
 502 *
 503          IF( JITER.GT.MAXIT )
 504      $      GO TO 180
 505 *
 506 *        Split the matrix if possible.
 507 *
 508 *        Two tests:
 509 *           1: H(j,j-1)=0  or  j=ILO
 510 *           2: T(j,j)=0
 511 *
 512 *        Special case: j=ILAST
 513 *
 514          IF( ILAST.EQ.ILO ) THEN
 515             GO TO 60
 516          ELSE
 517             IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
 518                H( ILAST, ILAST-1 ) = CZERO
 519                GO TO 60
 520             END IF
 521          END IF
 522 *
 523          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
 524             T( ILAST, ILAST ) = CZERO
 525             GO TO 50
 526          END IF
 527 *
 528 *        General case: j<ILAST
 529 *
 530          DO 40 J = ILAST - 1, ILO, -1
 531 *
 532 *           Test 1: for H(j,j-1)=0 or j=ILO
 533 *
 534             IF( J.EQ.ILO ) THEN
 535                ILAZRO = .TRUE.
 536             ELSE
 537                IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN
 538                   H( J, J-1 ) = CZERO
 539                   ILAZRO = .TRUE.
 540                ELSE
 541                   ILAZRO = .FALSE.
 542                END IF
 543             END IF
 544 *
 545 *           Test 2: for T(j,j)=0
 546 *
 547             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
 548                T( J, J ) = CZERO
 549 *
 550 *              Test 1a: Check for 2 consecutive small subdiagonals in A
 551 *
 552                ILAZR2 = .FALSE.
 553                IF( .NOT.ILAZRO ) THEN
 554                   IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1,
 555      $                J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) )
 556      $                ILAZR2 = .TRUE.
 557                END IF
 558 *
 559 *              If both tests pass (1 & 2), i.e., the leading diagonal
 560 *              element of B in the block is zero, split a 1x1 block off
 561 *              at the top. (I.e., at the J-th row/column) The leading
 562 *              diagonal element of the remainder can also be zero, so
 563 *              this may have to be done repeatedly.
 564 *
 565                IF( ILAZRO .OR. ILAZR2 ) THEN
 566                   DO 20 JCH = J, ILAST - 1
 567                      CTEMP = H( JCH, JCH )
 568                      CALL CLARTG( CTEMP, H( JCH+1, JCH ), C, S,
 569      $                            H( JCH, JCH ) )
 570                      H( JCH+1, JCH ) = CZERO
 571                      CALL CROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
 572      $                          H( JCH+1, JCH+1 ), LDH, C, S )
 573                      CALL CROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
 574      $                          T( JCH+1, JCH+1 ), LDT, C, S )
 575                      IF( ILQ )
 576      $                  CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
 577      $                             C, CONJG( S ) )
 578                      IF( ILAZR2 )
 579      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
 580                      ILAZR2 = .FALSE.
 581                      IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
 582                         IF( JCH+1.GE.ILAST ) THEN
 583                            GO TO 60
 584                         ELSE
 585                            IFIRST = JCH + 1
 586                            GO TO 70
 587                         END IF
 588                      END IF
 589                      T( JCH+1, JCH+1 ) = CZERO
 590    20             CONTINUE
 591                   GO TO 50
 592                ELSE
 593 *
 594 *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
 595 *                 Then process as in the case T(ILAST,ILAST)=0
 596 *
 597                   DO 30 JCH = J, ILAST - 1
 598                      CTEMP = T( JCH, JCH+1 )
 599                      CALL CLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S,
 600      $                            T( JCH, JCH+1 ) )
 601                      T( JCH+1, JCH+1 ) = CZERO
 602                      IF( JCH.LT.ILASTM-1 )
 603      $                  CALL CROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
 604      $                             T( JCH+1, JCH+2 ), LDT, C, S )
 605                      CALL CROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
 606      $                          H( JCH+1, JCH-1 ), LDH, C, S )
 607                      IF( ILQ )
 608      $                  CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
 609      $                             C, CONJG( S ) )
 610                      CTEMP = H( JCH+1, JCH )
 611                      CALL CLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S,
 612      $                            H( JCH+1, JCH ) )
 613                      H( JCH+1, JCH-1 ) = CZERO
 614                      CALL CROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
 615      $                          H( IFRSTM, JCH-1 ), 1, C, S )
 616                      CALL CROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
 617      $                          T( IFRSTM, JCH-1 ), 1, C, S )
 618                      IF( ILZ )
 619      $                  CALL CROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
 620      $                             C, S )
 621    30             CONTINUE
 622                   GO TO 50
 623                END IF
 624             ELSE IF( ILAZRO ) THEN
 625 *
 626 *              Only test 1 passed -- work on J:ILAST
 627 *
 628                IFIRST = J
 629                GO TO 70
 630             END IF
 631 *
 632 *           Neither test passed -- try next J
 633 *
 634    40    CONTINUE
 635 *
 636 *        (Drop-through is "impossible")
 637 *
 638          INFO = 2*N + 1
 639          GO TO 210
 640 *
 641 *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
 642 *        1x1 block.
 643 *
 644    50    CONTINUE
 645          CTEMP = H( ILAST, ILAST )
 646          CALL CLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S,
 647      $                H( ILAST, ILAST ) )
 648          H( ILAST, ILAST-1 ) = CZERO
 649          CALL CROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
 650      $              H( IFRSTM, ILAST-1 ), 1, C, S )
 651          CALL CROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
 652      $              T( IFRSTM, ILAST-1 ), 1, C, S )
 653          IF( ILZ )
 654      $      CALL CROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
 655 *
 656 *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
 657 *
 658    60    CONTINUE
 659          ABSB = ABS( T( ILAST, ILAST ) )
 660          IF( ABSB.GT.SAFMIN ) THEN
 661             SIGNBC = CONJG( T( ILAST, ILAST ) / ABSB )
 662             T( ILAST, ILAST ) = ABSB
 663             IF( ILSCHR ) THEN
 664                CALL CSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 )
 665                CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
 666      $                     1 )
 667             ELSE
 668                CALL CSCAL( 1, SIGNBC, H( ILAST, ILAST ), 1 )
 669             END IF
 670             IF( ILZ )
 671      $         CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
 672          ELSE
 673             T( ILAST, ILAST ) = CZERO
 674          END IF
 675          ALPHA( ILAST ) = H( ILAST, ILAST )
 676          BETA( ILAST ) = T( ILAST, ILAST )
 677 *
 678 *        Go to next block -- exit if finished.
 679 *
 680          ILAST = ILAST - 1
 681          IF( ILAST.LT.ILO )
 682      $      GO TO 190
 683 *
 684 *        Reset counters
 685 *
 686          IITER = 0
 687          ESHIFT = CZERO
 688          IF( .NOT.ILSCHR ) THEN
 689             ILASTM = ILAST
 690             IF( IFRSTM.GT.ILAST )
 691      $         IFRSTM = ILO
 692          END IF
 693          GO TO 160
 694 *
 695 *        QZ step
 696 *
 697 *        This iteration only involves rows/columns IFIRST:ILAST.  We
 698 *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
 699 *
 700    70    CONTINUE
 701          IITER = IITER + 1
 702          IF( .NOT.ILSCHR ) THEN
 703             IFRSTM = IFIRST
 704          END IF
 705 *
 706 *        Compute the Shift.
 707 *
 708 *        At this point, IFIRST < ILAST, and the diagonal elements of
 709 *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
 710 *        magnitude)
 711 *
 712          IF( ( IITER / 10 )*10.NE.IITER ) THEN
 713 *
 714 *           The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
 715 *           the bottom-right 2x2 block of A inv(B) which is nearest to
 716 *           the bottom-right element.
 717 *
 718 *           We factor B as U*D, where U has unit diagonals, and
 719 *           compute (A*inv(D))*inv(U).
 720 *
 721             U12 = ( BSCALE*T( ILAST-1, ILAST ) ) /
 722      $            ( BSCALE*T( ILAST, ILAST ) )
 723             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
 724      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
 725             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
 726      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
 727             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
 728      $             ( BSCALE*T( ILAST, ILAST ) )
 729             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
 730      $             ( BSCALE*T( ILAST, ILAST ) )
 731             ABI22 = AD22 - U12*AD21
 732 *
 733             T1 = HALF*( AD11+ABI22 )
 734             RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
 735             TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) +
 736      $             AIMAG( T1-ABI22 )*AIMAG( RTDISC )
 737             IF( TEMP.LE.ZERO ) THEN
 738                SHIFT = T1 + RTDISC
 739             ELSE
 740                SHIFT = T1 - RTDISC
 741             END IF
 742          ELSE
 743 *
 744 *           Exceptional shift.  Chosen for no particularly good reason.
 745 *
 746             ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
 747      $                        (BSCALE*T(ILAST-1,ILAST-1))
 748             SHIFT = ESHIFT
 749          END IF
 750 *
 751 *        Now check for two consecutive small subdiagonals.
 752 *
 753          DO 80 J = ILAST - 1, IFIRST + 1, -1
 754             ISTART = J
 755             CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) )
 756             TEMP = ABS1( CTEMP )
 757             TEMP2 = ASCALE*ABS1( H( J+1, J ) )
 758             TEMPR = MAX( TEMP, TEMP2 )
 759             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
 760                TEMP = TEMP / TEMPR
 761                TEMP2 = TEMP2 / TEMPR
 762             END IF
 763             IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL )
 764      $         GO TO 90
 765    80    CONTINUE
 766 *
 767          ISTART = IFIRST
 768          CTEMP = ASCALE*H( IFIRST, IFIRST ) -
 769      $           SHIFT*( BSCALE*T( IFIRST, IFIRST ) )
 770    90    CONTINUE
 771 *
 772 *        Do an implicit-shift QZ sweep.
 773 *
 774 *        Initial Q
 775 *
 776          CTEMP2 = ASCALE*H( ISTART+1, ISTART )
 777          CALL CLARTG( CTEMP, CTEMP2, C, S, CTEMP3 )
 778 *
 779 *        Sweep
 780 *
 781          DO 150 J = ISTART, ILAST - 1
 782             IF( J.GT.ISTART ) THEN
 783                CTEMP = H( J, J-1 )
 784                CALL CLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
 785                H( J+1, J-1 ) = CZERO
 786             END IF
 787 *
 788             DO 100 JC = J, ILASTM
 789                CTEMP = C*H( J, JC ) + S*H( J+1, JC )
 790                H( J+1, JC ) = -CONJG( S )*H( J, JC ) + C*H( J+1, JC )
 791                H( J, JC ) = CTEMP
 792                CTEMP2 = C*T( J, JC ) + S*T( J+1, JC )
 793                T( J+1, JC ) = -CONJG( S )*T( J, JC ) + C*T( J+1, JC )
 794                T( J, JC ) = CTEMP2
 795   100       CONTINUE
 796             IF( ILQ ) THEN
 797                DO 110 JR = 1, N
 798                   CTEMP = C*Q( JR, J ) + CONJG( S )*Q( JR, J+1 )
 799                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
 800                   Q( JR, J ) = CTEMP
 801   110          CONTINUE
 802             END IF
 803 *
 804             CTEMP = T( J+1, J+1 )
 805             CALL CLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
 806             T( J+1, J ) = CZERO
 807 *
 808             DO 120 JR = IFRSTM, MIN( J+2, ILAST )
 809                CTEMP = C*H( JR, J+1 ) + S*H( JR, J )
 810                H( JR, J ) = -CONJG( S )*H( JR, J+1 ) + C*H( JR, J )
 811                H( JR, J+1 ) = CTEMP
 812   120       CONTINUE
 813             DO 130 JR = IFRSTM, J
 814                CTEMP = C*T( JR, J+1 ) + S*T( JR, J )
 815                T( JR, J ) = -CONJG( S )*T( JR, J+1 ) + C*T( JR, J )
 816                T( JR, J+1 ) = CTEMP
 817   130       CONTINUE
 818             IF( ILZ ) THEN
 819                DO 140 JR = 1, N
 820                   CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
 821                   Z( JR, J ) = -CONJG( S )*Z( JR, J+1 ) + C*Z( JR, J )
 822                   Z( JR, J+1 ) = CTEMP
 823   140          CONTINUE
 824             END IF
 825   150    CONTINUE
 826 *
 827   160    CONTINUE
 828 *
 829   170 CONTINUE
 830 *
 831 *     Drop-through = non-convergence
 832 *
 833   180 CONTINUE
 834       INFO = ILAST
 835       GO TO 210
 836 *
 837 *     Successful completion of all QZ steps
 838 *
 839   190 CONTINUE
 840 *
 841 *     Set Eigenvalues 1:ILO-1
 842 *
 843       DO 200 J = 1, ILO - 1
 844          ABSB = ABS( T( J, J ) )
 845          IF( ABSB.GT.SAFMIN ) THEN
 846             SIGNBC = CONJG( T( J, J ) / ABSB )
 847             T( J, J ) = ABSB
 848             IF( ILSCHR ) THEN
 849                CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
 850                CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
 851             ELSE
 852                CALL CSCAL( 1, SIGNBC, H( J, J ), 1 )
 853             END IF
 854             IF( ILZ )
 855      $         CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
 856          ELSE
 857             T( J, J ) = CZERO
 858          END IF
 859          ALPHA( J ) = H( J, J )
 860          BETA( J ) = T( J, J )
 861   200 CONTINUE
 862 *
 863 *     Normal Termination
 864 *
 865       INFO = 0
 866 *
 867 *     Exit (other than argument error) -- return optimal workspace size
 868 *
 869   210 CONTINUE
 870       WORK( 1 ) = CMPLX( N )
 871       RETURN
 872 *
 873 *     End of CHGEQZ
 874 *
 875       END