LAPACK: SRC/zhgeqz.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
00002      $                   ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
00003      $                   RWORK, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          COMPQ, COMPZ, JOB
00012       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
00013 *     ..
00014 *     .. Array Arguments ..
00015       DOUBLE PRECISION   RWORK( * )
00016       COMPLEX*16         ALPHA( * ), BETA( * ), H( LDH, * ),
00017      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
00018      $                   Z( LDZ, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  ZHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
00025 *  where H is an upper Hessenberg matrix and T is upper triangular,
00026 *  using the single-shift QZ method.
00027 *  Matrix pairs of this type are produced by the reduction to
00028 *  generalized upper Hessenberg form of a complex matrix pair (A,B):
00029 *  
00030 *     A = Q1*H*Z1**H,  B = Q1*T*Z1**H,
00031 *  
00032 *  as computed by ZGGHRD.
00033 *  
00034 *  If JOB='S', then the Hessenberg-triangular pair (H,T) is
00035 *  also reduced to generalized Schur form,
00036 *  
00037 *     H = Q*S*Z**H,  T = Q*P*Z**H,
00038 *  
00039 *  where Q and Z are unitary matrices and S and P are upper triangular.
00040 *  
00041 *  Optionally, the unitary matrix Q from the generalized Schur
00042 *  factorization may be postmultiplied into an input matrix Q1, and the
00043 *  unitary matrix Z may be postmultiplied into an input matrix Z1.
00044 *  If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced
00045 *  the matrix pair (A,B) to generalized Hessenberg form, then the output
00046 *  matrices Q1*Q and Z1*Z are the unitary factors from the generalized
00047 *  Schur factorization of (A,B):
00048 *  
00049 *     A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.
00050 *  
00051 *  To avoid overflow, eigenvalues of the matrix pair (H,T)
00052 *  (equivalently, of (A,B)) are computed as a pair of complex values
00053 *  (alpha,beta).  If beta is nonzero, lambda = alpha / beta is an
00054 *  eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP)
00055 *     A*x = lambda*B*x
00056 *  and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
00057 *  alternate form of the GNEP
00058 *     mu*A*y = B*y.
00059 *  The values of alpha and beta for the i-th eigenvalue can be read
00060 *  directly from the generalized Schur form:  alpha = S(i,i),
00061 *  beta = P(i,i).
00062 *
00063 *  Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
00064 *       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
00065 *       pp. 241--256.
00066 *
00067 *  Arguments
00068 *  =========
00069 *
00070 *  JOB     (input) CHARACTER*1
00071 *          = 'E': Compute eigenvalues only;
00072 *          = 'S': Computer eigenvalues and the Schur form.
00073 *
00074 *  COMPQ   (input) CHARACTER*1
00075 *          = 'N': Left Schur vectors (Q) are not computed;
00076 *          = 'I': Q is initialized to the unit matrix and the matrix Q
00077 *                 of left Schur vectors of (H,T) is returned;
00078 *          = 'V': Q must contain a unitary matrix Q1 on entry and
00079 *                 the product Q1*Q is returned.
00080 *
00081 *  COMPZ   (input) CHARACTER*1
00082 *          = 'N': Right Schur vectors (Z) are not computed;
00083 *          = 'I': Q is initialized to the unit matrix and the matrix Z
00084 *                 of right Schur vectors of (H,T) is returned;
00085 *          = 'V': Z must contain a unitary matrix Z1 on entry and
00086 *                 the product Z1*Z is returned.
00087 *
00088 *  N       (input) INTEGER
00089 *          The order of the matrices H, T, Q, and Z.  N >= 0.
00090 *
00091 *  ILO     (input) INTEGER
00092 *  IHI     (input) INTEGER
00093 *          ILO and IHI mark the rows and columns of H which are in
00094 *          Hessenberg form.  It is assumed that A is already upper
00095 *          triangular in rows and columns 1:ILO-1 and IHI+1:N.
00096 *          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
00097 *
00098 *  H       (input/output) COMPLEX*16 array, dimension (LDH, N)
00099 *          On entry, the N-by-N upper Hessenberg matrix H.
00100 *          On exit, if JOB = 'S', H contains the upper triangular
00101 *          matrix S from the generalized Schur factorization.
00102 *          If JOB = 'E', the diagonal of H matches that of S, but
00103 *          the rest of H is unspecified.
00104 *
00105 *  LDH     (input) INTEGER
00106 *          The leading dimension of the array H.  LDH >= max( 1, N ).
00107 *
00108 *  T       (input/output) COMPLEX*16 array, dimension (LDT, N)
00109 *          On entry, the N-by-N upper triangular matrix T.
00110 *          On exit, if JOB = 'S', T contains the upper triangular
00111 *          matrix P from the generalized Schur factorization.
00112 *          If JOB = 'E', the diagonal of T matches that of P, but
00113 *          the rest of T is unspecified.
00114 *
00115 *  LDT     (input) INTEGER
00116 *          The leading dimension of the array T.  LDT >= max( 1, N ).
00117 *
00118 *  ALPHA   (output) COMPLEX*16 array, dimension (N)
00119 *          The complex scalars alpha that define the eigenvalues of
00120 *          GNEP.  ALPHA(i) = S(i,i) in the generalized Schur
00121 *          factorization.
00122 *
00123 *  BETA    (output) COMPLEX*16 array, dimension (N)
00124 *          The real non-negative scalars beta that define the
00125 *          eigenvalues of GNEP.  BETA(i) = P(i,i) in the generalized
00126 *          Schur factorization.
00127 *
00128 *          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
00129 *          represent the j-th eigenvalue of the matrix pair (A,B), in
00130 *          one of the forms lambda = alpha/beta or mu = beta/alpha.
00131 *          Since either lambda or mu may overflow, they should not,
00132 *          in general, be computed.
00133 *
00134 *  Q       (input/output) COMPLEX*16 array, dimension (LDQ, N)
00135 *          On entry, if COMPZ = 'V', the unitary matrix Q1 used in the
00136 *          reduction of (A,B) to generalized Hessenberg form.
00137 *          On exit, if COMPZ = 'I', the unitary matrix of left Schur
00138 *          vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
00139 *          left Schur vectors of (A,B).
00140 *          Not referenced if COMPZ = 'N'.
00141 *
00142 *  LDQ     (input) INTEGER
00143 *          The leading dimension of the array Q.  LDQ >= 1.
00144 *          If COMPQ='V' or 'I', then LDQ >= N.
00145 *
00146 *  Z       (input/output) COMPLEX*16 array, dimension (LDZ, N)
00147 *          On entry, if COMPZ = 'V', the unitary matrix Z1 used in the
00148 *          reduction of (A,B) to generalized Hessenberg form.
00149 *          On exit, if COMPZ = 'I', the unitary matrix of right Schur
00150 *          vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
00151 *          right Schur vectors of (A,B).
00152 *          Not referenced if COMPZ = 'N'.
00153 *
00154 *  LDZ     (input) INTEGER
00155 *          The leading dimension of the array Z.  LDZ >= 1.
00156 *          If COMPZ='V' or 'I', then LDZ >= N.
00157 *
00158 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
00159 *          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
00160 *
00161 *  LWORK   (input) INTEGER
00162 *          The dimension of the array WORK.  LWORK >= max(1,N).
00163 *
00164 *          If LWORK = -1, then a workspace query is assumed; the routine
00165 *          only calculates the optimal size of the WORK array, returns
00166 *          this value as the first entry of the WORK array, and no error
00167 *          message related to LWORK is issued by XERBLA.
00168 *
00169 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
00170 *
00171 *  INFO    (output) INTEGER
00172 *          = 0: successful exit
00173 *          < 0: if INFO = -i, the i-th argument had an illegal value
00174 *          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
00175 *                     in Schur form, but ALPHA(i) and BETA(i),
00176 *                     i=INFO+1,...,N should be correct.
00177 *          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
00178 *                     in Schur form, but ALPHA(i) and BETA(i),
00179 *                     i=INFO-N+1,...,N should be correct.
00180 *
00181 *  Further Details
00182 *  ===============
00183 *
00184 *  We assume that complex ABS works as long as its value is less than
00185 *  overflow.
00186 *
00187 *  =====================================================================
00188 *
00189 *     .. Parameters ..
00190       COMPLEX*16         CZERO, CONE
00191       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00192      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00193       DOUBLE PRECISION   ZERO, ONE
00194       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00195       DOUBLE PRECISION   HALF
00196       PARAMETER          ( HALF = 0.5D+0 )
00197 *     ..
00198 *     .. Local Scalars ..
00199       LOGICAL            ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
00200       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
00201      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
00202      $                   JR, MAXIT
00203       DOUBLE PRECISION   ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
00204      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
00205       COMPLEX*16         ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
00206      $                   CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
00207      $                   U12, X
00208 *     ..
00209 *     .. External Functions ..
00210       LOGICAL            LSAME
00211       DOUBLE PRECISION   DLAMCH, ZLANHS
00212       EXTERNAL           LSAME, DLAMCH, ZLANHS
00213 *     ..
00214 *     .. External Subroutines ..
00215       EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL
00216 *     ..
00217 *     .. Intrinsic Functions ..
00218       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN,
00219      $                   SQRT
00220 *     ..
00221 *     .. Statement Functions ..
00222       DOUBLE PRECISION   ABS1
00223 *     ..
00224 *     .. Statement Function definitions ..
00225       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
00226 *     ..
00227 *     .. Executable Statements ..
00228 *
00229 *     Decode JOB, COMPQ, COMPZ
00230 *
00231       IF( LSAME( JOB, 'E' ) ) THEN
00232          ILSCHR = .FALSE.
00233          ISCHUR = 1
00234       ELSE IF( LSAME( JOB, 'S' ) ) THEN
00235          ILSCHR = .TRUE.
00236          ISCHUR = 2
00237       ELSE
00238          ISCHUR = 0
00239       END IF
00240 *
00241       IF( LSAME( COMPQ, 'N' ) ) THEN
00242          ILQ = .FALSE.
00243          ICOMPQ = 1
00244       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
00245          ILQ = .TRUE.
00246          ICOMPQ = 2
00247       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
00248          ILQ = .TRUE.
00249          ICOMPQ = 3
00250       ELSE
00251          ICOMPQ = 0
00252       END IF
00253 *
00254       IF( LSAME( COMPZ, 'N' ) ) THEN
00255          ILZ = .FALSE.
00256          ICOMPZ = 1
00257       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00258          ILZ = .TRUE.
00259          ICOMPZ = 2
00260       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00261          ILZ = .TRUE.
00262          ICOMPZ = 3
00263       ELSE
00264          ICOMPZ = 0
00265       END IF
00266 *
00267 *     Check Argument Values
00268 *
00269       INFO = 0
00270       WORK( 1 ) = MAX( 1, N )
00271       LQUERY = ( LWORK.EQ.-1 )
00272       IF( ISCHUR.EQ.0 ) THEN
00273          INFO = -1
00274       ELSE IF( ICOMPQ.EQ.0 ) THEN
00275          INFO = -2
00276       ELSE IF( ICOMPZ.EQ.0 ) THEN
00277          INFO = -3
00278       ELSE IF( N.LT.0 ) THEN
00279          INFO = -4
00280       ELSE IF( ILO.LT.1 ) THEN
00281          INFO = -5
00282       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
00283          INFO = -6
00284       ELSE IF( LDH.LT.N ) THEN
00285          INFO = -8
00286       ELSE IF( LDT.LT.N ) THEN
00287          INFO = -10
00288       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
00289          INFO = -14
00290       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
00291          INFO = -16
00292       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00293          INFO = -18
00294       END IF
00295       IF( INFO.NE.0 ) THEN
00296          CALL XERBLA( 'ZHGEQZ', -INFO )
00297          RETURN
00298       ELSE IF( LQUERY ) THEN
00299          RETURN
00300       END IF
00301 *
00302 *     Quick return if possible
00303 *
00304 *     WORK( 1 ) = CMPLX( 1 )
00305       IF( N.LE.0 ) THEN
00306          WORK( 1 ) = DCMPLX( 1 )
00307          RETURN
00308       END IF
00309 *
00310 *     Initialize Q and Z
00311 *
00312       IF( ICOMPQ.EQ.3 )
00313      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
00314       IF( ICOMPZ.EQ.3 )
00315      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
00316 *
00317 *     Machine Constants
00318 *
00319       IN = IHI + 1 - ILO
00320       SAFMIN = DLAMCH( 'S' )
00321       ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
00322       ANORM = ZLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK )
00323       BNORM = ZLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK )
00324       ATOL = MAX( SAFMIN, ULP*ANORM )
00325       BTOL = MAX( SAFMIN, ULP*BNORM )
00326       ASCALE = ONE / MAX( SAFMIN, ANORM )
00327       BSCALE = ONE / MAX( SAFMIN, BNORM )
00328 *
00329 *
00330 *     Set Eigenvalues IHI+1:N
00331 *
00332       DO 10 J = IHI + 1, N
00333          ABSB = ABS( T( J, J ) )
00334          IF( ABSB.GT.SAFMIN ) THEN
00335             SIGNBC = DCONJG( T( J, J ) / ABSB )
00336             T( J, J ) = ABSB
00337             IF( ILSCHR ) THEN
00338                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
00339                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
00340             ELSE
00341                H( J, J ) = H( J, J )*SIGNBC
00342             END IF
00343             IF( ILZ )
00344      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
00345          ELSE
00346             T( J, J ) = CZERO
00347          END IF
00348          ALPHA( J ) = H( J, J )
00349          BETA( J ) = T( J, J )
00350    10 CONTINUE
00351 *
00352 *     If IHI < ILO, skip QZ steps
00353 *
00354       IF( IHI.LT.ILO )
00355      $   GO TO 190
00356 *
00357 *     MAIN QZ ITERATION LOOP
00358 *
00359 *     Initialize dynamic indices
00360 *
00361 *     Eigenvalues ILAST+1:N have been found.
00362 *        Column operations modify rows IFRSTM:whatever
00363 *        Row operations modify columns whatever:ILASTM
00364 *
00365 *     If only eigenvalues are being computed, then
00366 *        IFRSTM is the row of the last splitting row above row ILAST;
00367 *        this is always at least ILO.
00368 *     IITER counts iterations since the last eigenvalue was found,
00369 *        to tell when to use an extraordinary shift.
00370 *     MAXIT is the maximum number of QZ sweeps allowed.
00371 *
00372       ILAST = IHI
00373       IF( ILSCHR ) THEN
00374          IFRSTM = 1
00375          ILASTM = N
00376       ELSE
00377          IFRSTM = ILO
00378          ILASTM = IHI
00379       END IF
00380       IITER = 0
00381       ESHIFT = CZERO
00382       MAXIT = 30*( IHI-ILO+1 )
00383 *
00384       DO 170 JITER = 1, MAXIT
00385 *
00386 *        Check for too many iterations.
00387 *
00388          IF( JITER.GT.MAXIT )
00389      $      GO TO 180
00390 *
00391 *        Split the matrix if possible.
00392 *
00393 *        Two tests:
00394 *           1: H(j,j-1)=0  or  j=ILO
00395 *           2: T(j,j)=0
00396 *
00397 *        Special case: j=ILAST
00398 *
00399          IF( ILAST.EQ.ILO ) THEN
00400             GO TO 60
00401          ELSE
00402             IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
00403                H( ILAST, ILAST-1 ) = CZERO
00404                GO TO 60
00405             END IF
00406          END IF
00407 *
00408          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
00409             T( ILAST, ILAST ) = CZERO
00410             GO TO 50
00411          END IF
00412 *
00413 *        General case: j<ILAST
00414 *
00415          DO 40 J = ILAST - 1, ILO, -1
00416 *
00417 *           Test 1: for H(j,j-1)=0 or j=ILO
00418 *
00419             IF( J.EQ.ILO ) THEN
00420                ILAZRO = .TRUE.
00421             ELSE
00422                IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN
00423                   H( J, J-1 ) = CZERO
00424                   ILAZRO = .TRUE.
00425                ELSE
00426                   ILAZRO = .FALSE.
00427                END IF
00428             END IF
00429 *
00430 *           Test 2: for T(j,j)=0
00431 *
00432             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
00433                T( J, J ) = CZERO
00434 *
00435 *              Test 1a: Check for 2 consecutive small subdiagonals in A
00436 *
00437                ILAZR2 = .FALSE.
00438                IF( .NOT.ILAZRO ) THEN
00439                   IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1,
00440      $                J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) )
00441      $                ILAZR2 = .TRUE.
00442                END IF
00443 *
00444 *              If both tests pass (1 & 2), i.e., the leading diagonal
00445 *              element of B in the block is zero, split a 1x1 block off
00446 *              at the top. (I.e., at the J-th row/column) The leading
00447 *              diagonal element of the remainder can also be zero, so
00448 *              this may have to be done repeatedly.
00449 *
00450                IF( ILAZRO .OR. ILAZR2 ) THEN
00451                   DO 20 JCH = J, ILAST - 1
00452                      CTEMP = H( JCH, JCH )
00453                      CALL ZLARTG( CTEMP, H( JCH+1, JCH ), C, S,
00454      $                            H( JCH, JCH ) )
00455                      H( JCH+1, JCH ) = CZERO
00456                      CALL ZROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
00457      $                          H( JCH+1, JCH+1 ), LDH, C, S )
00458                      CALL ZROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
00459      $                          T( JCH+1, JCH+1 ), LDT, C, S )
00460                      IF( ILQ )
00461      $                  CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00462      $                             C, DCONJG( S ) )
00463                      IF( ILAZR2 )
00464      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
00465                      ILAZR2 = .FALSE.
00466                      IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
00467                         IF( JCH+1.GE.ILAST ) THEN
00468                            GO TO 60
00469                         ELSE
00470                            IFIRST = JCH + 1
00471                            GO TO 70
00472                         END IF
00473                      END IF
00474                      T( JCH+1, JCH+1 ) = CZERO
00475    20             CONTINUE
00476                   GO TO 50
00477                ELSE
00478 *
00479 *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
00480 *                 Then process as in the case T(ILAST,ILAST)=0
00481 *
00482                   DO 30 JCH = J, ILAST - 1
00483                      CTEMP = T( JCH, JCH+1 )
00484                      CALL ZLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S,
00485      $                            T( JCH, JCH+1 ) )
00486                      T( JCH+1, JCH+1 ) = CZERO
00487                      IF( JCH.LT.ILASTM-1 )
00488      $                  CALL ZROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
00489      $                             T( JCH+1, JCH+2 ), LDT, C, S )
00490                      CALL ZROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
00491      $                          H( JCH+1, JCH-1 ), LDH, C, S )
00492                      IF( ILQ )
00493      $                  CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00494      $                             C, DCONJG( S ) )
00495                      CTEMP = H( JCH+1, JCH )
00496                      CALL ZLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S,
00497      $                            H( JCH+1, JCH ) )
00498                      H( JCH+1, JCH-1 ) = CZERO
00499                      CALL ZROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
00500      $                          H( IFRSTM, JCH-1 ), 1, C, S )
00501                      CALL ZROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
00502      $                          T( IFRSTM, JCH-1 ), 1, C, S )
00503                      IF( ILZ )
00504      $                  CALL ZROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
00505      $                             C, S )
00506    30             CONTINUE
00507                   GO TO 50
00508                END IF
00509             ELSE IF( ILAZRO ) THEN
00510 *
00511 *              Only test 1 passed -- work on J:ILAST
00512 *
00513                IFIRST = J
00514                GO TO 70
00515             END IF
00516 *
00517 *           Neither test passed -- try next J
00518 *
00519    40    CONTINUE
00520 *
00521 *        (Drop-through is "impossible")
00522 *
00523          INFO = 2*N + 1
00524          GO TO 210
00525 *
00526 *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
00527 *        1x1 block.
00528 *
00529    50    CONTINUE
00530          CTEMP = H( ILAST, ILAST )
00531          CALL ZLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S,
00532      $                H( ILAST, ILAST ) )
00533          H( ILAST, ILAST-1 ) = CZERO
00534          CALL ZROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
00535      $              H( IFRSTM, ILAST-1 ), 1, C, S )
00536          CALL ZROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
00537      $              T( IFRSTM, ILAST-1 ), 1, C, S )
00538          IF( ILZ )
00539      $      CALL ZROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
00540 *
00541 *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
00542 *
00543    60    CONTINUE
00544          ABSB = ABS( T( ILAST, ILAST ) )
00545          IF( ABSB.GT.SAFMIN ) THEN
00546             SIGNBC = DCONJG( T( ILAST, ILAST ) / ABSB )
00547             T( ILAST, ILAST ) = ABSB
00548             IF( ILSCHR ) THEN
00549                CALL ZSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 )
00550                CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
00551      $                     1 )
00552             ELSE
00553                H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC
00554             END IF
00555             IF( ILZ )
00556      $         CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
00557          ELSE
00558             T( ILAST, ILAST ) = CZERO
00559          END IF
00560          ALPHA( ILAST ) = H( ILAST, ILAST )
00561          BETA( ILAST ) = T( ILAST, ILAST )
00562 *
00563 *        Go to next block -- exit if finished.
00564 *
00565          ILAST = ILAST - 1
00566          IF( ILAST.LT.ILO )
00567      $      GO TO 190
00568 *
00569 *        Reset counters
00570 *
00571          IITER = 0
00572          ESHIFT = CZERO
00573          IF( .NOT.ILSCHR ) THEN
00574             ILASTM = ILAST
00575             IF( IFRSTM.GT.ILAST )
00576      $         IFRSTM = ILO
00577          END IF
00578          GO TO 160
00579 *
00580 *        QZ step
00581 *
00582 *        This iteration only involves rows/columns IFIRST:ILAST.  We
00583 *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
00584 *
00585    70    CONTINUE
00586          IITER = IITER + 1
00587          IF( .NOT.ILSCHR ) THEN
00588             IFRSTM = IFIRST
00589          END IF
00590 *
00591 *        Compute the Shift.
00592 *
00593 *        At this point, IFIRST < ILAST, and the diagonal elements of
00594 *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
00595 *        magnitude)
00596 *
00597          IF( ( IITER / 10 )*10.NE.IITER ) THEN
00598 *
00599 *           The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
00600 *           the bottom-right 2x2 block of A inv(B) which is nearest to
00601 *           the bottom-right element.
00602 *
00603 *           We factor B as U*D, where U has unit diagonals, and
00604 *           compute (A*inv(D))*inv(U).
00605 *
00606             U12 = ( BSCALE*T( ILAST-1, ILAST ) ) /
00607      $            ( BSCALE*T( ILAST, ILAST ) )
00608             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
00609      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
00610             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
00611      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
00612             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
00613      $             ( BSCALE*T( ILAST, ILAST ) )
00614             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
00615      $             ( BSCALE*T( ILAST, ILAST ) )
00616             ABI22 = AD22 - U12*AD21
00617 *
00618             T1 = HALF*( AD11+ABI22 )
00619             RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
00620             TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) +
00621      $             DIMAG( T1-ABI22 )*DIMAG( RTDISC )
00622             IF( TEMP.LE.ZERO ) THEN
00623                SHIFT = T1 + RTDISC
00624             ELSE
00625                SHIFT = T1 - RTDISC
00626             END IF
00627          ELSE
00628 *
00629 *           Exceptional shift.  Chosen for no particularly good reason.
00630 *
00631             ESHIFT = ESHIFT + DCONJG( ( ASCALE*H( ILAST-1, ILAST ) ) /
00632      $               ( BSCALE*T( ILAST-1, ILAST-1 ) ) )
00633             SHIFT = ESHIFT
00634          END IF
00635 *
00636 *        Now check for two consecutive small subdiagonals.
00637 *
00638          DO 80 J = ILAST - 1, IFIRST + 1, -1
00639             ISTART = J
00640             CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) )
00641             TEMP = ABS1( CTEMP )
00642             TEMP2 = ASCALE*ABS1( H( J+1, J ) )
00643             TEMPR = MAX( TEMP, TEMP2 )
00644             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
00645                TEMP = TEMP / TEMPR
00646                TEMP2 = TEMP2 / TEMPR
00647             END IF
00648             IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL )
00649      $         GO TO 90
00650    80    CONTINUE
00651 *
00652          ISTART = IFIRST
00653          CTEMP = ASCALE*H( IFIRST, IFIRST ) -
00654      $           SHIFT*( BSCALE*T( IFIRST, IFIRST ) )
00655    90    CONTINUE
00656 *
00657 *        Do an implicit-shift QZ sweep.
00658 *
00659 *        Initial Q
00660 *
00661          CTEMP2 = ASCALE*H( ISTART+1, ISTART )
00662          CALL ZLARTG( CTEMP, CTEMP2, C, S, CTEMP3 )
00663 *
00664 *        Sweep
00665 *
00666          DO 150 J = ISTART, ILAST - 1
00667             IF( J.GT.ISTART ) THEN
00668                CTEMP = H( J, J-1 )
00669                CALL ZLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
00670                H( J+1, J-1 ) = CZERO
00671             END IF
00672 *
00673             DO 100 JC = J, ILASTM
00674                CTEMP = C*H( J, JC ) + S*H( J+1, JC )
00675                H( J+1, JC ) = -DCONJG( S )*H( J, JC ) + C*H( J+1, JC )
00676                H( J, JC ) = CTEMP
00677                CTEMP2 = C*T( J, JC ) + S*T( J+1, JC )
00678                T( J+1, JC ) = -DCONJG( S )*T( J, JC ) + C*T( J+1, JC )
00679                T( J, JC ) = CTEMP2
00680   100       CONTINUE
00681             IF( ILQ ) THEN
00682                DO 110 JR = 1, N
00683                   CTEMP = C*Q( JR, J ) + DCONJG( S )*Q( JR, J+1 )
00684                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
00685                   Q( JR, J ) = CTEMP
00686   110          CONTINUE
00687             END IF
00688 *
00689             CTEMP = T( J+1, J+1 )
00690             CALL ZLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
00691             T( J+1, J ) = CZERO
00692 *
00693             DO 120 JR = IFRSTM, MIN( J+2, ILAST )
00694                CTEMP = C*H( JR, J+1 ) + S*H( JR, J )
00695                H( JR, J ) = -DCONJG( S )*H( JR, J+1 ) + C*H( JR, J )
00696                H( JR, J+1 ) = CTEMP
00697   120       CONTINUE
00698             DO 130 JR = IFRSTM, J
00699                CTEMP = C*T( JR, J+1 ) + S*T( JR, J )
00700                T( JR, J ) = -DCONJG( S )*T( JR, J+1 ) + C*T( JR, J )
00701                T( JR, J+1 ) = CTEMP
00702   130       CONTINUE
00703             IF( ILZ ) THEN
00704                DO 140 JR = 1, N
00705                   CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
00706                   Z( JR, J ) = -DCONJG( S )*Z( JR, J+1 ) + C*Z( JR, J )
00707                   Z( JR, J+1 ) = CTEMP
00708   140          CONTINUE
00709             END IF
00710   150    CONTINUE
00711 *
00712   160    CONTINUE
00713 *
00714   170 CONTINUE
00715 *
00716 *     Drop-through = non-convergence
00717 *
00718   180 CONTINUE
00719       INFO = ILAST
00720       GO TO 210
00721 *
00722 *     Successful completion of all QZ steps
00723 *
00724   190 CONTINUE
00725 *
00726 *     Set Eigenvalues 1:ILO-1
00727 *
00728       DO 200 J = 1, ILO - 1
00729          ABSB = ABS( T( J, J ) )
00730          IF( ABSB.GT.SAFMIN ) THEN
00731             SIGNBC = DCONJG( T( J, J ) / ABSB )
00732             T( J, J ) = ABSB
00733             IF( ILSCHR ) THEN
00734                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
00735                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
00736             ELSE
00737                H( J, J ) = H( J, J )*SIGNBC
00738             END IF
00739             IF( ILZ )
00740      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
00741          ELSE
00742             T( J, J ) = CZERO
00743          END IF
00744          ALPHA( J ) = H( J, J )
00745          BETA( J ) = T( J, J )
00746   200 CONTINUE
00747 *
00748 *     Normal Termination
00749 *
00750       INFO = 0
00751 *
00752 *     Exit (other than argument error) -- return optimal workspace size
00753 *
00754   210 CONTINUE
00755       WORK( 1 ) = DCMPLX( N )
00756       RETURN
00757 *
00758 *     End of ZHGEQZ
00759 *
00760       END