LAPACK: SRC/chgeqz.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHGEQZ( 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       REAL               RWORK( * )
00016       COMPLEX            ALPHA( * ), BETA( * ), H( LDH, * ),
00017      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
00018      $                   Z( LDZ, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  CHGEQZ 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 CGGHRD.
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 CGGHRD 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 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 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 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 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 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 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 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) REAL 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            CZERO, CONE
00191       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00192      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00193       REAL               ZERO, ONE
00194       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00195       REAL               HALF
00196       PARAMETER          ( HALF = 0.5E+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       REAL               ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
00204      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
00205       COMPLEX            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       REAL               CLANHS, SLAMCH
00212       EXTERNAL           LSAME, CLANHS, SLAMCH
00213 *     ..
00214 *     .. External Subroutines ..
00215       EXTERNAL           CLARTG, CLASET, CROT, CSCAL, XERBLA
00216 *     ..
00217 *     .. Intrinsic Functions ..
00218       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT
00219 *     ..
00220 *     .. Statement Functions ..
00221       REAL               ABS1
00222 *     ..
00223 *     .. Statement Function definitions ..
00224       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
00225 *     ..
00226 *     .. Executable Statements ..
00227 *
00228 *     Decode JOB, COMPQ, COMPZ
00229 *
00230       IF( LSAME( JOB, 'E' ) ) THEN
00231          ILSCHR = .FALSE.
00232          ISCHUR = 1
00233       ELSE IF( LSAME( JOB, 'S' ) ) THEN
00234          ILSCHR = .TRUE.
00235          ISCHUR = 2
00236       ELSE
00237          ISCHUR = 0
00238       END IF
00239 *
00240       IF( LSAME( COMPQ, 'N' ) ) THEN
00241          ILQ = .FALSE.
00242          ICOMPQ = 1
00243       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
00244          ILQ = .TRUE.
00245          ICOMPQ = 2
00246       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
00247          ILQ = .TRUE.
00248          ICOMPQ = 3
00249       ELSE
00250          ICOMPQ = 0
00251       END IF
00252 *
00253       IF( LSAME( COMPZ, 'N' ) ) THEN
00254          ILZ = .FALSE.
00255          ICOMPZ = 1
00256       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00257          ILZ = .TRUE.
00258          ICOMPZ = 2
00259       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00260          ILZ = .TRUE.
00261          ICOMPZ = 3
00262       ELSE
00263          ICOMPZ = 0
00264       END IF
00265 *
00266 *     Check Argument Values
00267 *
00268       INFO = 0
00269       WORK( 1 ) = MAX( 1, N )
00270       LQUERY = ( LWORK.EQ.-1 )
00271       IF( ISCHUR.EQ.0 ) THEN
00272          INFO = -1
00273       ELSE IF( ICOMPQ.EQ.0 ) THEN
00274          INFO = -2
00275       ELSE IF( ICOMPZ.EQ.0 ) THEN
00276          INFO = -3
00277       ELSE IF( N.LT.0 ) THEN
00278          INFO = -4
00279       ELSE IF( ILO.LT.1 ) THEN
00280          INFO = -5
00281       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
00282          INFO = -6
00283       ELSE IF( LDH.LT.N ) THEN
00284          INFO = -8
00285       ELSE IF( LDT.LT.N ) THEN
00286          INFO = -10
00287       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
00288          INFO = -14
00289       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
00290          INFO = -16
00291       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00292          INFO = -18
00293       END IF
00294       IF( INFO.NE.0 ) THEN
00295          CALL XERBLA( 'CHGEQZ', -INFO )
00296          RETURN
00297       ELSE IF( LQUERY ) THEN
00298          RETURN
00299       END IF
00300 *
00301 *     Quick return if possible
00302 *
00303 c     WORK( 1 ) = CMPLX( 1 )
00304       IF( N.LE.0 ) THEN
00305          WORK( 1 ) = CMPLX( 1 )
00306          RETURN
00307       END IF
00308 *
00309 *     Initialize Q and Z
00310 *
00311       IF( ICOMPQ.EQ.3 )
00312      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
00313       IF( ICOMPZ.EQ.3 )
00314      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
00315 *
00316 *     Machine Constants
00317 *
00318       IN = IHI + 1 - ILO
00319       SAFMIN = SLAMCH( 'S' )
00320       ULP = SLAMCH( 'E' )*SLAMCH( 'B' )
00321       ANORM = CLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK )
00322       BNORM = CLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK )
00323       ATOL = MAX( SAFMIN, ULP*ANORM )
00324       BTOL = MAX( SAFMIN, ULP*BNORM )
00325       ASCALE = ONE / MAX( SAFMIN, ANORM )
00326       BSCALE = ONE / MAX( SAFMIN, BNORM )
00327 *
00328 *
00329 *     Set Eigenvalues IHI+1:N
00330 *
00331       DO 10 J = IHI + 1, N
00332          ABSB = ABS( T( J, J ) )
00333          IF( ABSB.GT.SAFMIN ) THEN
00334             SIGNBC = CONJG( T( J, J ) / ABSB )
00335             T( J, J ) = ABSB
00336             IF( ILSCHR ) THEN
00337                CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
00338                CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
00339             ELSE
00340                H( J, J ) = H( J, J )*SIGNBC
00341             END IF
00342             IF( ILZ )
00343      $         CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
00344          ELSE
00345             T( J, J ) = CZERO
00346          END IF
00347          ALPHA( J ) = H( J, J )
00348          BETA( J ) = T( J, J )
00349    10 CONTINUE
00350 *
00351 *     If IHI < ILO, skip QZ steps
00352 *
00353       IF( IHI.LT.ILO )
00354      $   GO TO 190
00355 *
00356 *     MAIN QZ ITERATION LOOP
00357 *
00358 *     Initialize dynamic indices
00359 *
00360 *     Eigenvalues ILAST+1:N have been found.
00361 *        Column operations modify rows IFRSTM:whatever
00362 *        Row operations modify columns whatever:ILASTM
00363 *
00364 *     If only eigenvalues are being computed, then
00365 *        IFRSTM is the row of the last splitting row above row ILAST;
00366 *        this is always at least ILO.
00367 *     IITER counts iterations since the last eigenvalue was found,
00368 *        to tell when to use an extraordinary shift.
00369 *     MAXIT is the maximum number of QZ sweeps allowed.
00370 *
00371       ILAST = IHI
00372       IF( ILSCHR ) THEN
00373          IFRSTM = 1
00374          ILASTM = N
00375       ELSE
00376          IFRSTM = ILO
00377          ILASTM = IHI
00378       END IF
00379       IITER = 0
00380       ESHIFT = CZERO
00381       MAXIT = 30*( IHI-ILO+1 )
00382 *
00383       DO 170 JITER = 1, MAXIT
00384 *
00385 *        Check for too many iterations.
00386 *
00387          IF( JITER.GT.MAXIT )
00388      $      GO TO 180
00389 *
00390 *        Split the matrix if possible.
00391 *
00392 *        Two tests:
00393 *           1: H(j,j-1)=0  or  j=ILO
00394 *           2: T(j,j)=0
00395 *
00396 *        Special case: j=ILAST
00397 *
00398          IF( ILAST.EQ.ILO ) THEN
00399             GO TO 60
00400          ELSE
00401             IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
00402                H( ILAST, ILAST-1 ) = CZERO
00403                GO TO 60
00404             END IF
00405          END IF
00406 *
00407          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
00408             T( ILAST, ILAST ) = CZERO
00409             GO TO 50
00410          END IF
00411 *
00412 *        General case: j<ILAST
00413 *
00414          DO 40 J = ILAST - 1, ILO, -1
00415 *
00416 *           Test 1: for H(j,j-1)=0 or j=ILO
00417 *
00418             IF( J.EQ.ILO ) THEN
00419                ILAZRO = .TRUE.
00420             ELSE
00421                IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN
00422                   H( J, J-1 ) = CZERO
00423                   ILAZRO = .TRUE.
00424                ELSE
00425                   ILAZRO = .FALSE.
00426                END IF
00427             END IF
00428 *
00429 *           Test 2: for T(j,j)=0
00430 *
00431             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
00432                T( J, J ) = CZERO
00433 *
00434 *              Test 1a: Check for 2 consecutive small subdiagonals in A
00435 *
00436                ILAZR2 = .FALSE.
00437                IF( .NOT.ILAZRO ) THEN
00438                   IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1,
00439      $                J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) )
00440      $                ILAZR2 = .TRUE.
00441                END IF
00442 *
00443 *              If both tests pass (1 & 2), i.e., the leading diagonal
00444 *              element of B in the block is zero, split a 1x1 block off
00445 *              at the top. (I.e., at the J-th row/column) The leading
00446 *              diagonal element of the remainder can also be zero, so
00447 *              this may have to be done repeatedly.
00448 *
00449                IF( ILAZRO .OR. ILAZR2 ) THEN
00450                   DO 20 JCH = J, ILAST - 1
00451                      CTEMP = H( JCH, JCH )
00452                      CALL CLARTG( CTEMP, H( JCH+1, JCH ), C, S,
00453      $                            H( JCH, JCH ) )
00454                      H( JCH+1, JCH ) = CZERO
00455                      CALL CROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
00456      $                          H( JCH+1, JCH+1 ), LDH, C, S )
00457                      CALL CROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
00458      $                          T( JCH+1, JCH+1 ), LDT, C, S )
00459                      IF( ILQ )
00460      $                  CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00461      $                             C, CONJG( S ) )
00462                      IF( ILAZR2 )
00463      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
00464                      ILAZR2 = .FALSE.
00465                      IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
00466                         IF( JCH+1.GE.ILAST ) THEN
00467                            GO TO 60
00468                         ELSE
00469                            IFIRST = JCH + 1
00470                            GO TO 70
00471                         END IF
00472                      END IF
00473                      T( JCH+1, JCH+1 ) = CZERO
00474    20             CONTINUE
00475                   GO TO 50
00476                ELSE
00477 *
00478 *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
00479 *                 Then process as in the case T(ILAST,ILAST)=0
00480 *
00481                   DO 30 JCH = J, ILAST - 1
00482                      CTEMP = T( JCH, JCH+1 )
00483                      CALL CLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S,
00484      $                            T( JCH, JCH+1 ) )
00485                      T( JCH+1, JCH+1 ) = CZERO
00486                      IF( JCH.LT.ILASTM-1 )
00487      $                  CALL CROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
00488      $                             T( JCH+1, JCH+2 ), LDT, C, S )
00489                      CALL CROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
00490      $                          H( JCH+1, JCH-1 ), LDH, C, S )
00491                      IF( ILQ )
00492      $                  CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00493      $                             C, CONJG( S ) )
00494                      CTEMP = H( JCH+1, JCH )
00495                      CALL CLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S,
00496      $                            H( JCH+1, JCH ) )
00497                      H( JCH+1, JCH-1 ) = CZERO
00498                      CALL CROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
00499      $                          H( IFRSTM, JCH-1 ), 1, C, S )
00500                      CALL CROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
00501      $                          T( IFRSTM, JCH-1 ), 1, C, S )
00502                      IF( ILZ )
00503      $                  CALL CROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
00504      $                             C, S )
00505    30             CONTINUE
00506                   GO TO 50
00507                END IF
00508             ELSE IF( ILAZRO ) THEN
00509 *
00510 *              Only test 1 passed -- work on J:ILAST
00511 *
00512                IFIRST = J
00513                GO TO 70
00514             END IF
00515 *
00516 *           Neither test passed -- try next J
00517 *
00518    40    CONTINUE
00519 *
00520 *        (Drop-through is "impossible")
00521 *
00522          INFO = 2*N + 1
00523          GO TO 210
00524 *
00525 *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
00526 *        1x1 block.
00527 *
00528    50    CONTINUE
00529          CTEMP = H( ILAST, ILAST )
00530          CALL CLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S,
00531      $                H( ILAST, ILAST ) )
00532          H( ILAST, ILAST-1 ) = CZERO
00533          CALL CROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
00534      $              H( IFRSTM, ILAST-1 ), 1, C, S )
00535          CALL CROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
00536      $              T( IFRSTM, ILAST-1 ), 1, C, S )
00537          IF( ILZ )
00538      $      CALL CROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
00539 *
00540 *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
00541 *
00542    60    CONTINUE
00543          ABSB = ABS( T( ILAST, ILAST ) )
00544          IF( ABSB.GT.SAFMIN ) THEN
00545             SIGNBC = CONJG( T( ILAST, ILAST ) / ABSB )
00546             T( ILAST, ILAST ) = ABSB
00547             IF( ILSCHR ) THEN
00548                CALL CSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 )
00549                CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
00550      $                     1 )
00551             ELSE
00552                H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC
00553             END IF
00554             IF( ILZ )
00555      $         CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
00556          ELSE
00557             T( ILAST, ILAST ) = CZERO
00558          END IF
00559          ALPHA( ILAST ) = H( ILAST, ILAST )
00560          BETA( ILAST ) = T( ILAST, ILAST )
00561 *
00562 *        Go to next block -- exit if finished.
00563 *
00564          ILAST = ILAST - 1
00565          IF( ILAST.LT.ILO )
00566      $      GO TO 190
00567 *
00568 *        Reset counters
00569 *
00570          IITER = 0
00571          ESHIFT = CZERO
00572          IF( .NOT.ILSCHR ) THEN
00573             ILASTM = ILAST
00574             IF( IFRSTM.GT.ILAST )
00575      $         IFRSTM = ILO
00576          END IF
00577          GO TO 160
00578 *
00579 *        QZ step
00580 *
00581 *        This iteration only involves rows/columns IFIRST:ILAST.  We
00582 *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
00583 *
00584    70    CONTINUE
00585          IITER = IITER + 1
00586          IF( .NOT.ILSCHR ) THEN
00587             IFRSTM = IFIRST
00588          END IF
00589 *
00590 *        Compute the Shift.
00591 *
00592 *        At this point, IFIRST < ILAST, and the diagonal elements of
00593 *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
00594 *        magnitude)
00595 *
00596          IF( ( IITER / 10 )*10.NE.IITER ) THEN
00597 *
00598 *           The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
00599 *           the bottom-right 2x2 block of A inv(B) which is nearest to
00600 *           the bottom-right element.
00601 *
00602 *           We factor B as U*D, where U has unit diagonals, and
00603 *           compute (A*inv(D))*inv(U).
00604 *
00605             U12 = ( BSCALE*T( ILAST-1, ILAST ) ) /
00606      $            ( BSCALE*T( ILAST, ILAST ) )
00607             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
00608      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
00609             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
00610      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
00611             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
00612      $             ( BSCALE*T( ILAST, ILAST ) )
00613             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
00614      $             ( BSCALE*T( ILAST, ILAST ) )
00615             ABI22 = AD22 - U12*AD21
00616 *
00617             T1 = HALF*( AD11+ABI22 )
00618             RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
00619             TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) +
00620      $             AIMAG( T1-ABI22 )*AIMAG( RTDISC )
00621             IF( TEMP.LE.ZERO ) THEN
00622                SHIFT = T1 + RTDISC
00623             ELSE
00624                SHIFT = T1 - RTDISC
00625             END IF
00626          ELSE
00627 *
00628 *           Exceptional shift.  Chosen for no particularly good reason.
00629 *
00630             ESHIFT = ESHIFT + CONJG( ( ASCALE*H( ILAST-1, ILAST ) ) /
00631      $               ( BSCALE*T( ILAST-1, ILAST-1 ) ) )
00632             SHIFT = ESHIFT
00633          END IF
00634 *
00635 *        Now check for two consecutive small subdiagonals.
00636 *
00637          DO 80 J = ILAST - 1, IFIRST + 1, -1
00638             ISTART = J
00639             CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) )
00640             TEMP = ABS1( CTEMP )
00641             TEMP2 = ASCALE*ABS1( H( J+1, J ) )
00642             TEMPR = MAX( TEMP, TEMP2 )
00643             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
00644                TEMP = TEMP / TEMPR
00645                TEMP2 = TEMP2 / TEMPR
00646             END IF
00647             IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL )
00648      $         GO TO 90
00649    80    CONTINUE
00650 *
00651          ISTART = IFIRST
00652          CTEMP = ASCALE*H( IFIRST, IFIRST ) -
00653      $           SHIFT*( BSCALE*T( IFIRST, IFIRST ) )
00654    90    CONTINUE
00655 *
00656 *        Do an implicit-shift QZ sweep.
00657 *
00658 *        Initial Q
00659 *
00660          CTEMP2 = ASCALE*H( ISTART+1, ISTART )
00661          CALL CLARTG( CTEMP, CTEMP2, C, S, CTEMP3 )
00662 *
00663 *        Sweep
00664 *
00665          DO 150 J = ISTART, ILAST - 1
00666             IF( J.GT.ISTART ) THEN
00667                CTEMP = H( J, J-1 )
00668                CALL CLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
00669                H( J+1, J-1 ) = CZERO
00670             END IF
00671 *
00672             DO 100 JC = J, ILASTM
00673                CTEMP = C*H( J, JC ) + S*H( J+1, JC )
00674                H( J+1, JC ) = -CONJG( S )*H( J, JC ) + C*H( J+1, JC )
00675                H( J, JC ) = CTEMP
00676                CTEMP2 = C*T( J, JC ) + S*T( J+1, JC )
00677                T( J+1, JC ) = -CONJG( S )*T( J, JC ) + C*T( J+1, JC )
00678                T( J, JC ) = CTEMP2
00679   100       CONTINUE
00680             IF( ILQ ) THEN
00681                DO 110 JR = 1, N
00682                   CTEMP = C*Q( JR, J ) + CONJG( S )*Q( JR, J+1 )
00683                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
00684                   Q( JR, J ) = CTEMP
00685   110          CONTINUE
00686             END IF
00687 *
00688             CTEMP = T( J+1, J+1 )
00689             CALL CLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
00690             T( J+1, J ) = CZERO
00691 *
00692             DO 120 JR = IFRSTM, MIN( J+2, ILAST )
00693                CTEMP = C*H( JR, J+1 ) + S*H( JR, J )
00694                H( JR, J ) = -CONJG( S )*H( JR, J+1 ) + C*H( JR, J )
00695                H( JR, J+1 ) = CTEMP
00696   120       CONTINUE
00697             DO 130 JR = IFRSTM, J
00698                CTEMP = C*T( JR, J+1 ) + S*T( JR, J )
00699                T( JR, J ) = -CONJG( S )*T( JR, J+1 ) + C*T( JR, J )
00700                T( JR, J+1 ) = CTEMP
00701   130       CONTINUE
00702             IF( ILZ ) THEN
00703                DO 140 JR = 1, N
00704                   CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
00705                   Z( JR, J ) = -CONJG( S )*Z( JR, J+1 ) + C*Z( JR, J )
00706                   Z( JR, J+1 ) = CTEMP
00707   140          CONTINUE
00708             END IF
00709   150    CONTINUE
00710 *
00711   160    CONTINUE
00712 *
00713   170 CONTINUE
00714 *
00715 *     Drop-through = non-convergence
00716 *
00717   180 CONTINUE
00718       INFO = ILAST
00719       GO TO 210
00720 *
00721 *     Successful completion of all QZ steps
00722 *
00723   190 CONTINUE
00724 *
00725 *     Set Eigenvalues 1:ILO-1
00726 *
00727       DO 200 J = 1, ILO - 1
00728          ABSB = ABS( T( J, J ) )
00729          IF( ABSB.GT.SAFMIN ) THEN
00730             SIGNBC = CONJG( T( J, J ) / ABSB )
00731             T( J, J ) = ABSB
00732             IF( ILSCHR ) THEN
00733                CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
00734                CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
00735             ELSE
00736                H( J, J ) = H( J, J )*SIGNBC
00737             END IF
00738             IF( ILZ )
00739      $         CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
00740          ELSE
00741             T( J, J ) = CZERO
00742          END IF
00743          ALPHA( J ) = H( J, J )
00744          BETA( J ) = T( J, J )
00745   200 CONTINUE
00746 *
00747 *     Normal Termination
00748 *
00749       INFO = 0
00750 *
00751 *     Exit (other than argument error) -- return optimal workspace size
00752 *
00753   210 CONTINUE
00754       WORK( 1 ) = CMPLX( N )
00755       RETURN
00756 *
00757 *     End of CHGEQZ
00758 *
00759       END