SRC/zhgeqz.f

   1 *> \brief \b ZHGEQZ
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHGEQZ + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhgeqz.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhgeqz.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhgeqz.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHGEQZ( 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 *       DOUBLE PRECISION   RWORK( * )
  31 *       COMPLEX*16         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 *> ZHGEQZ 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 ZGGHRD.
  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 ZGGHRD 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*16 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*16 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*16 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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 complex16GEcomputational
 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 ZHGEQZ( 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       DOUBLE PRECISION   RWORK( * )
 298       COMPLEX*16         ALPHA( * ), BETA( * ), H( LDH, * ),
 299      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
 300      $                   Z( LDZ, * )
 301 *     ..
 302 *
 303 *  =====================================================================
 304 *
 305 *     .. Parameters ..
 306       COMPLEX*16         CZERO, CONE
 307       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 308      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 309       DOUBLE PRECISION   ZERO, ONE
 310       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 311       DOUBLE PRECISION   HALF
 312       PARAMETER          ( HALF = 0.5D+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       DOUBLE PRECISION   ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
 320      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
 321       COMPLEX*16         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       DOUBLE PRECISION   DLAMCH, ZLANHS
 328       EXTERNAL           LSAME, DLAMCH, ZLANHS
 329 *     ..
 330 *     .. External Subroutines ..
 331       EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL
 332 *     ..
 333 *     .. Intrinsic Functions ..
 334       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN,
 335      $                   SQRT
 336 *     ..
 337 *     .. Statement Functions ..
 338       DOUBLE PRECISION   ABS1
 339 *     ..
 340 *     .. Statement Function definitions ..
 341       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
 342 *     ..
 343 *     .. Executable Statements ..
 344 *
 345 *     Decode JOB, COMPQ, COMPZ
 346 *
 347       IF( LSAME( JOB, 'E' ) ) THEN
 348          ILSCHR = .FALSE.
 349          ISCHUR = 1
 350       ELSE IF( LSAME( JOB, 'S' ) ) THEN
 351          ILSCHR = .TRUE.
 352          ISCHUR = 2
 353       ELSE
 354          ISCHUR = 0
 355       END IF
 356 *
 357       IF( LSAME( COMPQ, 'N' ) ) THEN
 358          ILQ = .FALSE.
 359          ICOMPQ = 1
 360       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
 361          ILQ = .TRUE.
 362          ICOMPQ = 2
 363       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
 364          ILQ = .TRUE.
 365          ICOMPQ = 3
 366       ELSE
 367          ICOMPQ = 0
 368       END IF
 369 *
 370       IF( LSAME( COMPZ, 'N' ) ) THEN
 371          ILZ = .FALSE.
 372          ICOMPZ = 1
 373       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 374          ILZ = .TRUE.
 375          ICOMPZ = 2
 376       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 377          ILZ = .TRUE.
 378          ICOMPZ = 3
 379       ELSE
 380          ICOMPZ = 0
 381       END IF
 382 *
 383 *     Check Argument Values
 384 *
 385       INFO = 0
 386       WORK( 1 ) = MAX( 1, N )
 387       LQUERY = ( LWORK.EQ.-1 )
 388       IF( ISCHUR.EQ.0 ) THEN
 389          INFO = -1
 390       ELSE IF( ICOMPQ.EQ.0 ) THEN
 391          INFO = -2
 392       ELSE IF( ICOMPZ.EQ.0 ) THEN
 393          INFO = -3
 394       ELSE IF( N.LT.0 ) THEN
 395          INFO = -4
 396       ELSE IF( ILO.LT.1 ) THEN
 397          INFO = -5
 398       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
 399          INFO = -6
 400       ELSE IF( LDH.LT.N ) THEN
 401          INFO = -8
 402       ELSE IF( LDT.LT.N ) THEN
 403          INFO = -10
 404       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
 405          INFO = -14
 406       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
 407          INFO = -16
 408       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 409          INFO = -18
 410       END IF
 411       IF( INFO.NE.0 ) THEN
 412          CALL XERBLA( 'ZHGEQZ', -INFO )
 413          RETURN
 414       ELSE IF( LQUERY ) THEN
 415          RETURN
 416       END IF
 417 *
 418 *     Quick return if possible
 419 *
 420 *     WORK( 1 ) = CMPLX( 1 )
 421       IF( N.LE.0 ) THEN
 422          WORK( 1 ) = DCMPLX( 1 )
 423          RETURN
 424       END IF
 425 *
 426 *     Initialize Q and Z
 427 *
 428       IF( ICOMPQ.EQ.3 )
 429      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
 430       IF( ICOMPZ.EQ.3 )
 431      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
 432 *
 433 *     Machine Constants
 434 *
 435       IN = IHI + 1 - ILO
 436       SAFMIN = DLAMCH( 'S' )
 437       ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
 438       ANORM = ZLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK )
 439       BNORM = ZLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK )
 440       ATOL = MAX( SAFMIN, ULP*ANORM )
 441       BTOL = MAX( SAFMIN, ULP*BNORM )
 442       ASCALE = ONE / MAX( SAFMIN, ANORM )
 443       BSCALE = ONE / MAX( SAFMIN, BNORM )
 444 *
 445 *
 446 *     Set Eigenvalues IHI+1:N
 447 *
 448       DO 10 J = IHI + 1, N
 449          ABSB = ABS( T( J, J ) )
 450          IF( ABSB.GT.SAFMIN ) THEN
 451             SIGNBC = DCONJG( T( J, J ) / ABSB )
 452             T( J, J ) = ABSB
 453             IF( ILSCHR ) THEN
 454                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
 455                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
 456             ELSE
 457                CALL ZSCAL( 1, SIGNBC, H( J, J ), 1 )
 458             END IF
 459             IF( ILZ )
 460      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
 461          ELSE
 462             T( J, J ) = CZERO
 463          END IF
 464          ALPHA( J ) = H( J, J )
 465          BETA( J ) = T( J, J )
 466    10 CONTINUE
 467 *
 468 *     If IHI < ILO, skip QZ steps
 469 *
 470       IF( IHI.LT.ILO )
 471      $   GO TO 190
 472 *
 473 *     MAIN QZ ITERATION LOOP
 474 *
 475 *     Initialize dynamic indices
 476 *
 477 *     Eigenvalues ILAST+1:N have been found.
 478 *        Column operations modify rows IFRSTM:whatever
 479 *        Row operations modify columns whatever:ILASTM
 480 *
 481 *     If only eigenvalues are being computed, then
 482 *        IFRSTM is the row of the last splitting row above row ILAST;
 483 *        this is always at least ILO.
 484 *     IITER counts iterations since the last eigenvalue was found,
 485 *        to tell when to use an extraordinary shift.
 486 *     MAXIT is the maximum number of QZ sweeps allowed.
 487 *
 488       ILAST = IHI
 489       IF( ILSCHR ) THEN
 490          IFRSTM = 1
 491          ILASTM = N
 492       ELSE
 493          IFRSTM = ILO
 494          ILASTM = IHI
 495       END IF
 496       IITER = 0
 497       ESHIFT = CZERO
 498       MAXIT = 30*( IHI-ILO+1 )
 499 *
 500       DO 170 JITER = 1, MAXIT
 501 *
 502 *        Check for too many iterations.
 503 *
 504          IF( JITER.GT.MAXIT )
 505      $      GO TO 180
 506 *
 507 *        Split the matrix if possible.
 508 *
 509 *        Two tests:
 510 *           1: H(j,j-1)=0  or  j=ILO
 511 *           2: T(j,j)=0
 512 *
 513 *        Special case: j=ILAST
 514 *
 515          IF( ILAST.EQ.ILO ) THEN
 516             GO TO 60
 517          ELSE
 518             IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
 519                H( ILAST, ILAST-1 ) = CZERO
 520                GO TO 60
 521             END IF
 522          END IF
 523 *
 524          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
 525             T( ILAST, ILAST ) = CZERO
 526             GO TO 50
 527          END IF
 528 *
 529 *        General case: j<ILAST
 530 *
 531          DO 40 J = ILAST - 1, ILO, -1
 532 *
 533 *           Test 1: for H(j,j-1)=0 or j=ILO
 534 *
 535             IF( J.EQ.ILO ) THEN
 536                ILAZRO = .TRUE.
 537             ELSE
 538                IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN
 539                   H( J, J-1 ) = CZERO
 540                   ILAZRO = .TRUE.
 541                ELSE
 542                   ILAZRO = .FALSE.
 543                END IF
 544             END IF
 545 *
 546 *           Test 2: for T(j,j)=0
 547 *
 548             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
 549                T( J, J ) = CZERO
 550 *
 551 *              Test 1a: Check for 2 consecutive small subdiagonals in A
 552 *
 553                ILAZR2 = .FALSE.
 554                IF( .NOT.ILAZRO ) THEN
 555                   IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1,
 556      $                J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) )
 557      $                ILAZR2 = .TRUE.
 558                END IF
 559 *
 560 *              If both tests pass (1 & 2), i.e., the leading diagonal
 561 *              element of B in the block is zero, split a 1x1 block off
 562 *              at the top. (I.e., at the J-th row/column) The leading
 563 *              diagonal element of the remainder can also be zero, so
 564 *              this may have to be done repeatedly.
 565 *
 566                IF( ILAZRO .OR. ILAZR2 ) THEN
 567                   DO 20 JCH = J, ILAST - 1
 568                      CTEMP = H( JCH, JCH )
 569                      CALL ZLARTG( CTEMP, H( JCH+1, JCH ), C, S,
 570      $                            H( JCH, JCH ) )
 571                      H( JCH+1, JCH ) = CZERO
 572                      CALL ZROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
 573      $                          H( JCH+1, JCH+1 ), LDH, C, S )
 574                      CALL ZROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
 575      $                          T( JCH+1, JCH+1 ), LDT, C, S )
 576                      IF( ILQ )
 577      $                  CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
 578      $                             C, DCONJG( S ) )
 579                      IF( ILAZR2 )
 580      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
 581                      ILAZR2 = .FALSE.
 582                      IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
 583                         IF( JCH+1.GE.ILAST ) THEN
 584                            GO TO 60
 585                         ELSE
 586                            IFIRST = JCH + 1
 587                            GO TO 70
 588                         END IF
 589                      END IF
 590                      T( JCH+1, JCH+1 ) = CZERO
 591    20             CONTINUE
 592                   GO TO 50
 593                ELSE
 594 *
 595 *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
 596 *                 Then process as in the case T(ILAST,ILAST)=0
 597 *
 598                   DO 30 JCH = J, ILAST - 1
 599                      CTEMP = T( JCH, JCH+1 )
 600                      CALL ZLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S,
 601      $                            T( JCH, JCH+1 ) )
 602                      T( JCH+1, JCH+1 ) = CZERO
 603                      IF( JCH.LT.ILASTM-1 )
 604      $                  CALL ZROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
 605      $                             T( JCH+1, JCH+2 ), LDT, C, S )
 606                      CALL ZROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
 607      $                          H( JCH+1, JCH-1 ), LDH, C, S )
 608                      IF( ILQ )
 609      $                  CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
 610      $                             C, DCONJG( S ) )
 611                      CTEMP = H( JCH+1, JCH )
 612                      CALL ZLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S,
 613      $                            H( JCH+1, JCH ) )
 614                      H( JCH+1, JCH-1 ) = CZERO
 615                      CALL ZROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
 616      $                          H( IFRSTM, JCH-1 ), 1, C, S )
 617                      CALL ZROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
 618      $                          T( IFRSTM, JCH-1 ), 1, C, S )
 619                      IF( ILZ )
 620      $                  CALL ZROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
 621      $                             C, S )
 622    30             CONTINUE
 623                   GO TO 50
 624                END IF
 625             ELSE IF( ILAZRO ) THEN
 626 *
 627 *              Only test 1 passed -- work on J:ILAST
 628 *
 629                IFIRST = J
 630                GO TO 70
 631             END IF
 632 *
 633 *           Neither test passed -- try next J
 634 *
 635    40    CONTINUE
 636 *
 637 *        (Drop-through is "impossible")
 638 *
 639          INFO = 2*N + 1
 640          GO TO 210
 641 *
 642 *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
 643 *        1x1 block.
 644 *
 645    50    CONTINUE
 646          CTEMP = H( ILAST, ILAST )
 647          CALL ZLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S,
 648      $                H( ILAST, ILAST ) )
 649          H( ILAST, ILAST-1 ) = CZERO
 650          CALL ZROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
 651      $              H( IFRSTM, ILAST-1 ), 1, C, S )
 652          CALL ZROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
 653      $              T( IFRSTM, ILAST-1 ), 1, C, S )
 654          IF( ILZ )
 655      $      CALL ZROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
 656 *
 657 *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
 658 *
 659    60    CONTINUE
 660          ABSB = ABS( T( ILAST, ILAST ) )
 661          IF( ABSB.GT.SAFMIN ) THEN
 662             SIGNBC = DCONJG( T( ILAST, ILAST ) / ABSB )
 663             T( ILAST, ILAST ) = ABSB
 664             IF( ILSCHR ) THEN
 665                CALL ZSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 )
 666                CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
 667      $                     1 )
 668             ELSE
 669                CALL ZSCAL( 1, SIGNBC, H( ILAST, ILAST ), 1 )
 670             END IF
 671             IF( ILZ )
 672      $         CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
 673          ELSE
 674             T( ILAST, ILAST ) = CZERO
 675          END IF
 676          ALPHA( ILAST ) = H( ILAST, ILAST )
 677          BETA( ILAST ) = T( ILAST, ILAST )
 678 *
 679 *        Go to next block -- exit if finished.
 680 *
 681          ILAST = ILAST - 1
 682          IF( ILAST.LT.ILO )
 683      $      GO TO 190
 684 *
 685 *        Reset counters
 686 *
 687          IITER = 0
 688          ESHIFT = CZERO
 689          IF( .NOT.ILSCHR ) THEN
 690             ILASTM = ILAST
 691             IF( IFRSTM.GT.ILAST )
 692      $         IFRSTM = ILO
 693          END IF
 694          GO TO 160
 695 *
 696 *        QZ step
 697 *
 698 *        This iteration only involves rows/columns IFIRST:ILAST.  We
 699 *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
 700 *
 701    70    CONTINUE
 702          IITER = IITER + 1
 703          IF( .NOT.ILSCHR ) THEN
 704             IFRSTM = IFIRST
 705          END IF
 706 *
 707 *        Compute the Shift.
 708 *
 709 *        At this point, IFIRST < ILAST, and the diagonal elements of
 710 *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
 711 *        magnitude)
 712 *
 713          IF( ( IITER / 10 )*10.NE.IITER ) THEN
 714 *
 715 *           The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
 716 *           the bottom-right 2x2 block of A inv(B) which is nearest to
 717 *           the bottom-right element.
 718 *
 719 *           We factor B as U*D, where U has unit diagonals, and
 720 *           compute (A*inv(D))*inv(U).
 721 *
 722             U12 = ( BSCALE*T( ILAST-1, ILAST ) ) /
 723      $            ( BSCALE*T( ILAST, ILAST ) )
 724             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
 725      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
 726             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
 727      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
 728             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
 729      $             ( BSCALE*T( ILAST, ILAST ) )
 730             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
 731      $             ( BSCALE*T( ILAST, ILAST ) )
 732             ABI22 = AD22 - U12*AD21
 733 *
 734             T1 = HALF*( AD11+ABI22 )
 735             RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
 736             TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) +
 737      $             DIMAG( T1-ABI22 )*DIMAG( RTDISC )
 738             IF( TEMP.LE.ZERO ) THEN
 739                SHIFT = T1 + RTDISC
 740             ELSE
 741                SHIFT = T1 - RTDISC
 742             END IF
 743          ELSE
 744 *
 745 *           Exceptional shift.  Chosen for no particularly good reason.
 746 *
 747             ESHIFT = ESHIFT + (ASCALE*H(ILAST,ILAST-1))/
 748      $                        (BSCALE*T(ILAST-1,ILAST-1))
 749             SHIFT = ESHIFT
 750          END IF
 751 *
 752 *        Now check for two consecutive small subdiagonals.
 753 *
 754          DO 80 J = ILAST - 1, IFIRST + 1, -1
 755             ISTART = J
 756             CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) )
 757             TEMP = ABS1( CTEMP )
 758             TEMP2 = ASCALE*ABS1( H( J+1, J ) )
 759             TEMPR = MAX( TEMP, TEMP2 )
 760             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
 761                TEMP = TEMP / TEMPR
 762                TEMP2 = TEMP2 / TEMPR
 763             END IF
 764             IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL )
 765      $         GO TO 90
 766    80    CONTINUE
 767 *
 768          ISTART = IFIRST
 769          CTEMP = ASCALE*H( IFIRST, IFIRST ) -
 770      $           SHIFT*( BSCALE*T( IFIRST, IFIRST ) )
 771    90    CONTINUE
 772 *
 773 *        Do an implicit-shift QZ sweep.
 774 *
 775 *        Initial Q
 776 *
 777          CTEMP2 = ASCALE*H( ISTART+1, ISTART )
 778          CALL ZLARTG( CTEMP, CTEMP2, C, S, CTEMP3 )
 779 *
 780 *        Sweep
 781 *
 782          DO 150 J = ISTART, ILAST - 1
 783             IF( J.GT.ISTART ) THEN
 784                CTEMP = H( J, J-1 )
 785                CALL ZLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
 786                H( J+1, J-1 ) = CZERO
 787             END IF
 788 *
 789             DO 100 JC = J, ILASTM
 790                CTEMP = C*H( J, JC ) + S*H( J+1, JC )
 791                H( J+1, JC ) = -DCONJG( S )*H( J, JC ) + C*H( J+1, JC )
 792                H( J, JC ) = CTEMP
 793                CTEMP2 = C*T( J, JC ) + S*T( J+1, JC )
 794                T( J+1, JC ) = -DCONJG( S )*T( J, JC ) + C*T( J+1, JC )
 795                T( J, JC ) = CTEMP2
 796   100       CONTINUE
 797             IF( ILQ ) THEN
 798                DO 110 JR = 1, N
 799                   CTEMP = C*Q( JR, J ) + DCONJG( S )*Q( JR, J+1 )
 800                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
 801                   Q( JR, J ) = CTEMP
 802   110          CONTINUE
 803             END IF
 804 *
 805             CTEMP = T( J+1, J+1 )
 806             CALL ZLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
 807             T( J+1, J ) = CZERO
 808 *
 809             DO 120 JR = IFRSTM, MIN( J+2, ILAST )
 810                CTEMP = C*H( JR, J+1 ) + S*H( JR, J )
 811                H( JR, J ) = -DCONJG( S )*H( JR, J+1 ) + C*H( JR, J )
 812                H( JR, J+1 ) = CTEMP
 813   120       CONTINUE
 814             DO 130 JR = IFRSTM, J
 815                CTEMP = C*T( JR, J+1 ) + S*T( JR, J )
 816                T( JR, J ) = -DCONJG( S )*T( JR, J+1 ) + C*T( JR, J )
 817                T( JR, J+1 ) = CTEMP
 818   130       CONTINUE
 819             IF( ILZ ) THEN
 820                DO 140 JR = 1, N
 821                   CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
 822                   Z( JR, J ) = -DCONJG( S )*Z( JR, J+1 ) + C*Z( JR, J )
 823                   Z( JR, J+1 ) = CTEMP
 824   140          CONTINUE
 825             END IF
 826   150    CONTINUE
 827 *
 828   160    CONTINUE
 829 *
 830   170 CONTINUE
 831 *
 832 *     Drop-through = non-convergence
 833 *
 834   180 CONTINUE
 835       INFO = ILAST
 836       GO TO 210
 837 *
 838 *     Successful completion of all QZ steps
 839 *
 840   190 CONTINUE
 841 *
 842 *     Set Eigenvalues 1:ILO-1
 843 *
 844       DO 200 J = 1, ILO - 1
 845          ABSB = ABS( T( J, J ) )
 846          IF( ABSB.GT.SAFMIN ) THEN
 847             SIGNBC = DCONJG( T( J, J ) / ABSB )
 848             T( J, J ) = ABSB
 849             IF( ILSCHR ) THEN
 850                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
 851                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
 852             ELSE
 853                CALL ZSCAL( 1, SIGNBC, H( J, J ), 1 )
 854             END IF
 855             IF( ILZ )
 856      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
 857          ELSE
 858             T( J, J ) = CZERO
 859          END IF
 860          ALPHA( J ) = H( J, J )
 861          BETA( J ) = T( J, J )
 862   200 CONTINUE
 863 *
 864 *     Normal Termination
 865 *
 866       INFO = 0
 867 *
 868 *     Exit (other than argument error) -- return optimal workspace size
 869 *
 870   210 CONTINUE
 871       WORK( 1 ) = DCMPLX( N )
 872       RETURN
 873 *
 874 *     End of ZHGEQZ
 875 *
 876       END