LAPACK: SRC/shgeqz.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
00002      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
00003      $                   LWORK, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2.1)                                  --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *  -- April 2009                                                      --
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               ALPHAI( * ), ALPHAR( * ), BETA( * ),
00016      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
00017      $                   WORK( * ), Z( LDZ, * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  SHGEQZ computes the eigenvalues of a real matrix pair (H,T),
00024 *  where H is an upper Hessenberg matrix and T is upper triangular,
00025 *  using the double-shift QZ method.
00026 *  Matrix pairs of this type are produced by the reduction to
00027 *  generalized upper Hessenberg form of a real matrix pair (A,B):
00028 *
00029 *     A = Q1*H*Z1**T,  B = Q1*T*Z1**T,
00030 *
00031 *  as computed by SGGHRD.
00032 *
00033 *  If JOB='S', then the Hessenberg-triangular pair (H,T) is
00034 *  also reduced to generalized Schur form,
00035 *  
00036 *     H = Q*S*Z**T,  T = Q*P*Z**T,
00037 *  
00038 *  where Q and Z are orthogonal matrices, P is an upper triangular
00039 *  matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
00040 *  diagonal blocks.
00041 *
00042 *  The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
00043 *  (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
00044 *  eigenvalues.
00045 *
00046 *  Additionally, the 2-by-2 upper triangular diagonal blocks of P
00047 *  corresponding to 2-by-2 blocks of S are reduced to positive diagonal
00048 *  form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
00049 *  P(j,j) > 0, and P(j+1,j+1) > 0.
00050 *
00051 *  Optionally, the orthogonal matrix Q from the generalized Schur
00052 *  factorization may be postmultiplied into an input matrix Q1, and the
00053 *  orthogonal matrix Z may be postmultiplied into an input matrix Z1.
00054 *  If Q1 and Z1 are the orthogonal matrices from SGGHRD that reduced
00055 *  the matrix pair (A,B) to generalized upper Hessenberg form, then the
00056 *  output matrices Q1*Q and Z1*Z are the orthogonal factors from the
00057 *  generalized Schur factorization of (A,B):
00058 *
00059 *     A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.
00060 *  
00061 *  To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
00062 *  of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
00063 *  complex and beta real.
00064 *  If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
00065 *  generalized nonsymmetric eigenvalue problem (GNEP)
00066 *     A*x = lambda*B*x
00067 *  and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
00068 *  alternate form of the GNEP
00069 *     mu*A*y = B*y.
00070 *  Real eigenvalues can be read directly from the generalized Schur
00071 *  form: 
00072 *    alpha = S(i,i), beta = P(i,i).
00073 *
00074 *  Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
00075 *       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
00076 *       pp. 241--256.
00077 *
00078 *  Arguments
00079 *  =========
00080 *
00081 *  JOB     (input) CHARACTER*1
00082 *          = 'E': Compute eigenvalues only;
00083 *          = 'S': Compute eigenvalues and the Schur form. 
00084 *
00085 *  COMPQ   (input) CHARACTER*1
00086 *          = 'N': Left Schur vectors (Q) are not computed;
00087 *          = 'I': Q is initialized to the unit matrix and the matrix Q
00088 *                 of left Schur vectors of (H,T) is returned;
00089 *          = 'V': Q must contain an orthogonal matrix Q1 on entry and
00090 *                 the product Q1*Q is returned.
00091 *
00092 *  COMPZ   (input) CHARACTER*1
00093 *          = 'N': Right Schur vectors (Z) are not computed;
00094 *          = 'I': Z is initialized to the unit matrix and the matrix Z
00095 *                 of right Schur vectors of (H,T) is returned;
00096 *          = 'V': Z must contain an orthogonal matrix Z1 on entry and
00097 *                 the product Z1*Z is returned.
00098 *
00099 *  N       (input) INTEGER
00100 *          The order of the matrices H, T, Q, and Z.  N >= 0.
00101 *
00102 *  ILO     (input) INTEGER
00103 *  IHI     (input) INTEGER
00104 *          ILO and IHI mark the rows and columns of H which are in
00105 *          Hessenberg form.  It is assumed that A is already upper
00106 *          triangular in rows and columns 1:ILO-1 and IHI+1:N.
00107 *          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
00108 *
00109 *  H       (input/output) REAL array, dimension (LDH, N)
00110 *          On entry, the N-by-N upper Hessenberg matrix H.
00111 *          On exit, if JOB = 'S', H contains the upper quasi-triangular
00112 *          matrix S from the generalized Schur factorization;
00113 *          2-by-2 diagonal blocks (corresponding to complex conjugate
00114 *          pairs of eigenvalues) are returned in standard form, with
00115 *          H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0.
00116 *          If JOB = 'E', the diagonal blocks of H match those of S, but
00117 *          the rest of H is unspecified.
00118 *
00119 *  LDH     (input) INTEGER
00120 *          The leading dimension of the array H.  LDH >= max( 1, N ).
00121 *
00122 *  T       (input/output) REAL array, dimension (LDT, N)
00123 *          On entry, the N-by-N upper triangular matrix T.
00124 *          On exit, if JOB = 'S', T contains the upper triangular
00125 *          matrix P from the generalized Schur factorization;
00126 *          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
00127 *          are reduced to positive diagonal form, i.e., if H(j+1,j) is
00128 *          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
00129 *          T(j+1,j+1) > 0.
00130 *          If JOB = 'E', the diagonal blocks of T match those of P, but
00131 *          the rest of T is unspecified.
00132 *
00133 *  LDT     (input) INTEGER
00134 *          The leading dimension of the array T.  LDT >= max( 1, N ).
00135 *
00136 *  ALPHAR  (output) REAL array, dimension (N)
00137 *          The real parts of each scalar alpha defining an eigenvalue
00138 *          of GNEP.
00139 *
00140 *  ALPHAI  (output) REAL array, dimension (N)
00141 *          The imaginary parts of each scalar alpha defining an
00142 *          eigenvalue of GNEP.
00143 *          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
00144 *          positive, then the j-th and (j+1)-st eigenvalues are a
00145 *          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
00146 *
00147 *  BETA    (output) REAL array, dimension (N)
00148 *          The scalars beta that define the eigenvalues of GNEP.
00149 *          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
00150 *          beta = BETA(j) represent the j-th eigenvalue of the matrix
00151 *          pair (A,B), in one of the forms lambda = alpha/beta or
00152 *          mu = beta/alpha.  Since either lambda or mu may overflow,
00153 *          they should not, in general, be computed.
00154 *
00155 *  Q       (input/output) REAL array, dimension (LDQ, N)
00156 *          On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in
00157 *          the reduction of (A,B) to generalized Hessenberg form.
00158 *          On exit, if COMPZ = 'I', the orthogonal matrix of left Schur
00159 *          vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix
00160 *          of left Schur vectors of (A,B).
00161 *          Not referenced if COMPZ = 'N'.
00162 *
00163 *  LDQ     (input) INTEGER
00164 *          The leading dimension of the array Q.  LDQ >= 1.
00165 *          If COMPQ='V' or 'I', then LDQ >= N.
00166 *
00167 *  Z       (input/output) REAL array, dimension (LDZ, N)
00168 *          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
00169 *          the reduction of (A,B) to generalized Hessenberg form.
00170 *          On exit, if COMPZ = 'I', the orthogonal matrix of
00171 *          right Schur vectors of (H,T), and if COMPZ = 'V', the
00172 *          orthogonal matrix of right Schur vectors of (A,B).
00173 *          Not referenced if COMPZ = 'N'.
00174 *
00175 *  LDZ     (input) INTEGER
00176 *          The leading dimension of the array Z.  LDZ >= 1.
00177 *          If COMPZ='V' or 'I', then LDZ >= N.
00178 *
00179 *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
00180 *          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
00181 *
00182 *  LWORK   (input) INTEGER
00183 *          The dimension of the array WORK.  LWORK >= max(1,N).
00184 *
00185 *          If LWORK = -1, then a workspace query is assumed; the routine
00186 *          only calculates the optimal size of the WORK array, returns
00187 *          this value as the first entry of the WORK array, and no error
00188 *          message related to LWORK is issued by XERBLA.
00189 *
00190 *  INFO    (output) INTEGER
00191 *          = 0: successful exit
00192 *          < 0: if INFO = -i, the i-th argument had an illegal value
00193 *          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
00194 *                     in Schur form, but ALPHAR(i), ALPHAI(i), and
00195 *                     BETA(i), i=INFO+1,...,N should be correct.
00196 *          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
00197 *                     in Schur form, but ALPHAR(i), ALPHAI(i), and
00198 *                     BETA(i), i=INFO-N+1,...,N should be correct.
00199 *
00200 *  Further Details
00201 *  ===============
00202 *
00203 *  Iteration counters:
00204 *
00205 *  JITER  -- counts iterations.
00206 *  IITER  -- counts iterations run since ILAST was last
00207 *            changed.  This is therefore reset only when a 1-by-1 or
00208 *            2-by-2 block deflates off the bottom.
00209 *
00210 *  =====================================================================
00211 *
00212 *     .. Parameters ..
00213 *    $                     SAFETY = 1.0E+0 )
00214       REAL               HALF, ZERO, ONE, SAFETY
00215       PARAMETER          ( HALF = 0.5E+0, ZERO = 0.0E+0, ONE = 1.0E+0,
00216      $                   SAFETY = 1.0E+2 )
00217 *     ..
00218 *     .. Local Scalars ..
00219       LOGICAL            ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
00220      $                   LQUERY
00221       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
00222      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
00223      $                   JR, MAXIT
00224       REAL               A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
00225      $                   AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
00226      $                   AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
00227      $                   B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
00228      $                   BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
00229      $                   CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
00230      $                   SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
00231      $                   TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
00232      $                   U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
00233      $                   WR2
00234 *     ..
00235 *     .. Local Arrays ..
00236       REAL               V( 3 )
00237 *     ..
00238 *     .. External Functions ..
00239       LOGICAL            LSAME
00240       REAL               SLAMCH, SLANHS, SLAPY2, SLAPY3
00241       EXTERNAL           LSAME, SLAMCH, SLANHS, SLAPY2, SLAPY3
00242 *     ..
00243 *     .. External Subroutines ..
00244       EXTERNAL           SLAG2, SLARFG, SLARTG, SLASET, SLASV2, SROT,
00245      $                   XERBLA
00246 *     ..
00247 *     .. Intrinsic Functions ..
00248       INTRINSIC          ABS, MAX, MIN, REAL, SQRT
00249 *     ..
00250 *     .. Executable Statements ..
00251 *
00252 *     Decode JOB, COMPQ, COMPZ
00253 *
00254       IF( LSAME( JOB, 'E' ) ) THEN
00255          ILSCHR = .FALSE.
00256          ISCHUR = 1
00257       ELSE IF( LSAME( JOB, 'S' ) ) THEN
00258          ILSCHR = .TRUE.
00259          ISCHUR = 2
00260       ELSE
00261          ISCHUR = 0
00262       END IF
00263 *
00264       IF( LSAME( COMPQ, 'N' ) ) THEN
00265          ILQ = .FALSE.
00266          ICOMPQ = 1
00267       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
00268          ILQ = .TRUE.
00269          ICOMPQ = 2
00270       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
00271          ILQ = .TRUE.
00272          ICOMPQ = 3
00273       ELSE
00274          ICOMPQ = 0
00275       END IF
00276 *
00277       IF( LSAME( COMPZ, 'N' ) ) THEN
00278          ILZ = .FALSE.
00279          ICOMPZ = 1
00280       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
00281          ILZ = .TRUE.
00282          ICOMPZ = 2
00283       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
00284          ILZ = .TRUE.
00285          ICOMPZ = 3
00286       ELSE
00287          ICOMPZ = 0
00288       END IF
00289 *
00290 *     Check Argument Values
00291 *
00292       INFO = 0
00293       WORK( 1 ) = MAX( 1, N )
00294       LQUERY = ( LWORK.EQ.-1 )
00295       IF( ISCHUR.EQ.0 ) THEN
00296          INFO = -1
00297       ELSE IF( ICOMPQ.EQ.0 ) THEN
00298          INFO = -2
00299       ELSE IF( ICOMPZ.EQ.0 ) THEN
00300          INFO = -3
00301       ELSE IF( N.LT.0 ) THEN
00302          INFO = -4
00303       ELSE IF( ILO.LT.1 ) THEN
00304          INFO = -5
00305       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
00306          INFO = -6
00307       ELSE IF( LDH.LT.N ) THEN
00308          INFO = -8
00309       ELSE IF( LDT.LT.N ) THEN
00310          INFO = -10
00311       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
00312          INFO = -15
00313       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
00314          INFO = -17
00315       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00316          INFO = -19
00317       END IF
00318       IF( INFO.NE.0 ) THEN
00319          CALL XERBLA( 'SHGEQZ', -INFO )
00320          RETURN
00321       ELSE IF( LQUERY ) THEN
00322          RETURN
00323       END IF
00324 *
00325 *     Quick return if possible
00326 *
00327       IF( N.LE.0 ) THEN
00328          WORK( 1 ) = REAL( 1 )
00329          RETURN
00330       END IF
00331 *
00332 *     Initialize Q and Z
00333 *
00334       IF( ICOMPQ.EQ.3 )
00335      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
00336       IF( ICOMPZ.EQ.3 )
00337      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
00338 *
00339 *     Machine Constants
00340 *
00341       IN = IHI + 1 - ILO
00342       SAFMIN = SLAMCH( 'S' )
00343       SAFMAX = ONE / SAFMIN
00344       ULP = SLAMCH( 'E' )*SLAMCH( 'B' )
00345       ANORM = SLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK )
00346       BNORM = SLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK )
00347       ATOL = MAX( SAFMIN, ULP*ANORM )
00348       BTOL = MAX( SAFMIN, ULP*BNORM )
00349       ASCALE = ONE / MAX( SAFMIN, ANORM )
00350       BSCALE = ONE / MAX( SAFMIN, BNORM )
00351 *
00352 *     Set Eigenvalues IHI+1:N
00353 *
00354       DO 30 J = IHI + 1, N
00355          IF( T( J, J ).LT.ZERO ) THEN
00356             IF( ILSCHR ) THEN
00357                DO 10 JR = 1, J
00358                   H( JR, J ) = -H( JR, J )
00359                   T( JR, J ) = -T( JR, J )
00360    10          CONTINUE
00361             ELSE
00362                H( J, J ) = -H( J, J )
00363                T( J, J ) = -T( J, J )
00364             END IF
00365             IF( ILZ ) THEN
00366                DO 20 JR = 1, N
00367                   Z( JR, J ) = -Z( JR, J )
00368    20          CONTINUE
00369             END IF
00370          END IF
00371          ALPHAR( J ) = H( J, J )
00372          ALPHAI( J ) = ZERO
00373          BETA( J ) = T( J, J )
00374    30 CONTINUE
00375 *
00376 *     If IHI < ILO, skip QZ steps
00377 *
00378       IF( IHI.LT.ILO )
00379      $   GO TO 380
00380 *
00381 *     MAIN QZ ITERATION LOOP
00382 *
00383 *     Initialize dynamic indices
00384 *
00385 *     Eigenvalues ILAST+1:N have been found.
00386 *        Column operations modify rows IFRSTM:whatever.
00387 *        Row operations modify columns whatever:ILASTM.
00388 *
00389 *     If only eigenvalues are being computed, then
00390 *        IFRSTM is the row of the last splitting row above row ILAST;
00391 *        this is always at least ILO.
00392 *     IITER counts iterations since the last eigenvalue was found,
00393 *        to tell when to use an extraordinary shift.
00394 *     MAXIT is the maximum number of QZ sweeps allowed.
00395 *
00396       ILAST = IHI
00397       IF( ILSCHR ) THEN
00398          IFRSTM = 1
00399          ILASTM = N
00400       ELSE
00401          IFRSTM = ILO
00402          ILASTM = IHI
00403       END IF
00404       IITER = 0
00405       ESHIFT = ZERO
00406       MAXIT = 30*( IHI-ILO+1 )
00407 *
00408       DO 360 JITER = 1, MAXIT
00409 *
00410 *        Split the matrix if possible.
00411 *
00412 *        Two tests:
00413 *           1: H(j,j-1)=0  or  j=ILO
00414 *           2: T(j,j)=0
00415 *
00416          IF( ILAST.EQ.ILO ) THEN
00417 *
00418 *           Special case: j=ILAST
00419 *
00420             GO TO 80
00421          ELSE
00422             IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
00423                H( ILAST, ILAST-1 ) = ZERO
00424                GO TO 80
00425             END IF
00426          END IF
00427 *
00428          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
00429             T( ILAST, ILAST ) = ZERO
00430             GO TO 70
00431          END IF
00432 *
00433 *        General case: j<ILAST
00434 *
00435          DO 60 J = ILAST - 1, ILO, -1
00436 *
00437 *           Test 1: for H(j,j-1)=0 or j=ILO
00438 *
00439             IF( J.EQ.ILO ) THEN
00440                ILAZRO = .TRUE.
00441             ELSE
00442                IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN
00443                   H( J, J-1 ) = ZERO
00444                   ILAZRO = .TRUE.
00445                ELSE
00446                   ILAZRO = .FALSE.
00447                END IF
00448             END IF
00449 *
00450 *           Test 2: for T(j,j)=0
00451 *
00452             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
00453                T( J, J ) = ZERO
00454 *
00455 *              Test 1a: Check for 2 consecutive small subdiagonals in A
00456 *
00457                ILAZR2 = .FALSE.
00458                IF( .NOT.ILAZRO ) THEN
00459                   TEMP = ABS( H( J, J-1 ) )
00460                   TEMP2 = ABS( H( J, J ) )
00461                   TEMPR = MAX( TEMP, TEMP2 )
00462                   IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
00463                      TEMP = TEMP / TEMPR
00464                      TEMP2 = TEMP2 / TEMPR
00465                   END IF
00466                   IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2*
00467      $                ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
00468                END IF
00469 *
00470 *              If both tests pass (1 & 2), i.e., the leading diagonal
00471 *              element of B in the block is zero, split a 1x1 block off
00472 *              at the top. (I.e., at the J-th row/column) The leading
00473 *              diagonal element of the remainder can also be zero, so
00474 *              this may have to be done repeatedly.
00475 *
00476                IF( ILAZRO .OR. ILAZR2 ) THEN
00477                   DO 40 JCH = J, ILAST - 1
00478                      TEMP = H( JCH, JCH )
00479                      CALL SLARTG( TEMP, H( JCH+1, JCH ), C, S,
00480      $                            H( JCH, JCH ) )
00481                      H( JCH+1, JCH ) = ZERO
00482                      CALL SROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
00483      $                          H( JCH+1, JCH+1 ), LDH, C, S )
00484                      CALL SROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
00485      $                          T( JCH+1, JCH+1 ), LDT, C, S )
00486                      IF( ILQ )
00487      $                  CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00488      $                             C, S )
00489                      IF( ILAZR2 )
00490      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
00491                      ILAZR2 = .FALSE.
00492                      IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
00493                         IF( JCH+1.GE.ILAST ) THEN
00494                            GO TO 80
00495                         ELSE
00496                            IFIRST = JCH + 1
00497                            GO TO 110
00498                         END IF
00499                      END IF
00500                      T( JCH+1, JCH+1 ) = ZERO
00501    40             CONTINUE
00502                   GO TO 70
00503                ELSE
00504 *
00505 *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
00506 *                 Then process as in the case T(ILAST,ILAST)=0
00507 *
00508                   DO 50 JCH = J, ILAST - 1
00509                      TEMP = T( JCH, JCH+1 )
00510                      CALL SLARTG( TEMP, T( JCH+1, JCH+1 ), C, S,
00511      $                            T( JCH, JCH+1 ) )
00512                      T( JCH+1, JCH+1 ) = ZERO
00513                      IF( JCH.LT.ILASTM-1 )
00514      $                  CALL SROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
00515      $                             T( JCH+1, JCH+2 ), LDT, C, S )
00516                      CALL SROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
00517      $                          H( JCH+1, JCH-1 ), LDH, C, S )
00518                      IF( ILQ )
00519      $                  CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
00520      $                             C, S )
00521                      TEMP = H( JCH+1, JCH )
00522                      CALL SLARTG( TEMP, H( JCH+1, JCH-1 ), C, S,
00523      $                            H( JCH+1, JCH ) )
00524                      H( JCH+1, JCH-1 ) = ZERO
00525                      CALL SROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
00526      $                          H( IFRSTM, JCH-1 ), 1, C, S )
00527                      CALL SROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
00528      $                          T( IFRSTM, JCH-1 ), 1, C, S )
00529                      IF( ILZ )
00530      $                  CALL SROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
00531      $                             C, S )
00532    50             CONTINUE
00533                   GO TO 70
00534                END IF
00535             ELSE IF( ILAZRO ) THEN
00536 *
00537 *              Only test 1 passed -- work on J:ILAST
00538 *
00539                IFIRST = J
00540                GO TO 110
00541             END IF
00542 *
00543 *           Neither test passed -- try next J
00544 *
00545    60    CONTINUE
00546 *
00547 *        (Drop-through is "impossible")
00548 *
00549          INFO = N + 1
00550          GO TO 420
00551 *
00552 *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
00553 *        1x1 block.
00554 *
00555    70    CONTINUE
00556          TEMP = H( ILAST, ILAST )
00557          CALL SLARTG( TEMP, H( ILAST, ILAST-1 ), C, S,
00558      $                H( ILAST, ILAST ) )
00559          H( ILAST, ILAST-1 ) = ZERO
00560          CALL SROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
00561      $              H( IFRSTM, ILAST-1 ), 1, C, S )
00562          CALL SROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
00563      $              T( IFRSTM, ILAST-1 ), 1, C, S )
00564          IF( ILZ )
00565      $      CALL SROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
00566 *
00567 *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
00568 *                              and BETA
00569 *
00570    80    CONTINUE
00571          IF( T( ILAST, ILAST ).LT.ZERO ) THEN
00572             IF( ILSCHR ) THEN
00573                DO 90 J = IFRSTM, ILAST
00574                   H( J, ILAST ) = -H( J, ILAST )
00575                   T( J, ILAST ) = -T( J, ILAST )
00576    90          CONTINUE
00577             ELSE
00578                H( ILAST, ILAST ) = -H( ILAST, ILAST )
00579                T( ILAST, ILAST ) = -T( ILAST, ILAST )
00580             END IF
00581             IF( ILZ ) THEN
00582                DO 100 J = 1, N
00583                   Z( J, ILAST ) = -Z( J, ILAST )
00584   100          CONTINUE
00585             END IF
00586          END IF
00587          ALPHAR( ILAST ) = H( ILAST, ILAST )
00588          ALPHAI( ILAST ) = ZERO
00589          BETA( ILAST ) = T( ILAST, ILAST )
00590 *
00591 *        Go to next block -- exit if finished.
00592 *
00593          ILAST = ILAST - 1
00594          IF( ILAST.LT.ILO )
00595      $      GO TO 380
00596 *
00597 *        Reset counters
00598 *
00599          IITER = 0
00600          ESHIFT = ZERO
00601          IF( .NOT.ILSCHR ) THEN
00602             ILASTM = ILAST
00603             IF( IFRSTM.GT.ILAST )
00604      $         IFRSTM = ILO
00605          END IF
00606          GO TO 350
00607 *
00608 *        QZ step
00609 *
00610 *        This iteration only involves rows/columns IFIRST:ILAST. We
00611 *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
00612 *
00613   110    CONTINUE
00614          IITER = IITER + 1
00615          IF( .NOT.ILSCHR ) THEN
00616             IFRSTM = IFIRST
00617          END IF
00618 *
00619 *        Compute single shifts.
00620 *
00621 *        At this point, IFIRST < ILAST, and the diagonal elements of
00622 *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
00623 *        magnitude)
00624 *
00625          IF( ( IITER / 10 )*10.EQ.IITER ) THEN
00626 *
00627 *           Exceptional shift.  Chosen for no particularly good reason.
00628 *           (Single shift only.)
00629 *
00630             IF( ( REAL( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
00631      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
00632                ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /
00633      $                  T( ILAST-1, ILAST-1 )
00634             ELSE
00635                ESHIFT = ESHIFT + ONE / ( SAFMIN*REAL( MAXIT ) )
00636             END IF
00637             S1 = ONE
00638             WR = ESHIFT
00639 *
00640          ELSE
00641 *
00642 *           Shifts based on the generalized eigenvalues of the
00643 *           bottom-right 2x2 block of A and B. The first eigenvalue
00644 *           returned by SLAG2 is the Wilkinson shift (AEP p.512),
00645 *
00646             CALL SLAG2( H( ILAST-1, ILAST-1 ), LDH,
00647      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
00648      $                  S2, WR, WR2, WI )
00649 *
00650             TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
00651             IF( WI.NE.ZERO )
00652      $         GO TO 200
00653          END IF
00654 *
00655 *        Fiddle with shift to avoid overflow
00656 *
00657          TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
00658          IF( S1.GT.TEMP ) THEN
00659             SCALE = TEMP / S1
00660          ELSE
00661             SCALE = ONE
00662          END IF
00663 *
00664          TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
00665          IF( ABS( WR ).GT.TEMP )
00666      $      SCALE = MIN( SCALE, TEMP / ABS( WR ) )
00667          S1 = SCALE*S1
00668          WR = SCALE*WR
00669 *
00670 *        Now check for two consecutive small subdiagonals.
00671 *
00672          DO 120 J = ILAST - 1, IFIRST + 1, -1
00673             ISTART = J
00674             TEMP = ABS( S1*H( J, J-1 ) )
00675             TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) )
00676             TEMPR = MAX( TEMP, TEMP2 )
00677             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
00678                TEMP = TEMP / TEMPR
00679                TEMP2 = TEMP2 / TEMPR
00680             END IF
00681             IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
00682      $          TEMP2 )GO TO 130
00683   120    CONTINUE
00684 *
00685          ISTART = IFIRST
00686   130    CONTINUE
00687 *
00688 *        Do an implicit single-shift QZ sweep.
00689 *
00690 *        Initial Q
00691 *
00692          TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART )
00693          TEMP2 = S1*H( ISTART+1, ISTART )
00694          CALL SLARTG( TEMP, TEMP2, C, S, TEMPR )
00695 *
00696 *        Sweep
00697 *
00698          DO 190 J = ISTART, ILAST - 1
00699             IF( J.GT.ISTART ) THEN
00700                TEMP = H( J, J-1 )
00701                CALL SLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
00702                H( J+1, J-1 ) = ZERO
00703             END IF
00704 *
00705             DO 140 JC = J, ILASTM
00706                TEMP = C*H( J, JC ) + S*H( J+1, JC )
00707                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
00708                H( J, JC ) = TEMP
00709                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
00710                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
00711                T( J, JC ) = TEMP2
00712   140       CONTINUE
00713             IF( ILQ ) THEN
00714                DO 150 JR = 1, N
00715                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
00716                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
00717                   Q( JR, J ) = TEMP
00718   150          CONTINUE
00719             END IF
00720 *
00721             TEMP = T( J+1, J+1 )
00722             CALL SLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
00723             T( J+1, J ) = ZERO
00724 *
00725             DO 160 JR = IFRSTM, MIN( J+2, ILAST )
00726                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
00727                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
00728                H( JR, J+1 ) = TEMP
00729   160       CONTINUE
00730             DO 170 JR = IFRSTM, J
00731                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
00732                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
00733                T( JR, J+1 ) = TEMP
00734   170       CONTINUE
00735             IF( ILZ ) THEN
00736                DO 180 JR = 1, N
00737                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
00738                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
00739                   Z( JR, J+1 ) = TEMP
00740   180          CONTINUE
00741             END IF
00742   190    CONTINUE
00743 *
00744          GO TO 350
00745 *
00746 *        Use Francis double-shift
00747 *
00748 *        Note: the Francis double-shift should work with real shifts,
00749 *              but only if the block is at least 3x3.
00750 *              This code may break if this point is reached with
00751 *              a 2x2 block with real eigenvalues.
00752 *
00753   200    CONTINUE
00754          IF( IFIRST+1.EQ.ILAST ) THEN
00755 *
00756 *           Special case -- 2x2 block with complex eigenvectors
00757 *
00758 *           Step 1: Standardize, that is, rotate so that
00759 *
00760 *                       ( B11  0  )
00761 *                   B = (         )  with B11 non-negative.
00762 *                       (  0  B22 )
00763 *
00764             CALL SLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ),
00765      $                   T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
00766 *
00767             IF( B11.LT.ZERO ) THEN
00768                CR = -CR
00769                SR = -SR
00770                B11 = -B11
00771                B22 = -B22
00772             END IF
00773 *
00774             CALL SROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,
00775      $                 H( ILAST, ILAST-1 ), LDH, CL, SL )
00776             CALL SROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
00777      $                 H( IFRSTM, ILAST ), 1, CR, SR )
00778 *
00779             IF( ILAST.LT.ILASTM )
00780      $         CALL SROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
00781      $                    T( ILAST, ILAST+1 ), LDT, CL, SL )
00782             IF( IFRSTM.LT.ILAST-1 )
00783      $         CALL SROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
00784      $                    T( IFRSTM, ILAST ), 1, CR, SR )
00785 *
00786             IF( ILQ )
00787      $         CALL SROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
00788      $                    SL )
00789             IF( ILZ )
00790      $         CALL SROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
00791      $                    SR )
00792 *
00793             T( ILAST-1, ILAST-1 ) = B11
00794             T( ILAST-1, ILAST ) = ZERO
00795             T( ILAST, ILAST-1 ) = ZERO
00796             T( ILAST, ILAST ) = B22
00797 *
00798 *           If B22 is negative, negate column ILAST
00799 *
00800             IF( B22.LT.ZERO ) THEN
00801                DO 210 J = IFRSTM, ILAST
00802                   H( J, ILAST ) = -H( J, ILAST )
00803                   T( J, ILAST ) = -T( J, ILAST )
00804   210          CONTINUE
00805 *
00806                IF( ILZ ) THEN
00807                   DO 220 J = 1, N
00808                      Z( J, ILAST ) = -Z( J, ILAST )
00809   220             CONTINUE
00810                END IF
00811             END IF
00812 *
00813 *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
00814 *
00815 *           Recompute shift
00816 *
00817             CALL SLAG2( H( ILAST-1, ILAST-1 ), LDH,
00818      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
00819      $                  TEMP, WR, TEMP2, WI )
00820 *
00821 *           If standardization has perturbed the shift onto real line,
00822 *           do another (real single-shift) QR step.
00823 *
00824             IF( WI.EQ.ZERO )
00825      $         GO TO 350
00826             S1INV = ONE / S1
00827 *
00828 *           Do EISPACK (QZVAL) computation of alpha and beta
00829 *
00830             A11 = H( ILAST-1, ILAST-1 )
00831             A21 = H( ILAST, ILAST-1 )
00832             A12 = H( ILAST-1, ILAST )
00833             A22 = H( ILAST, ILAST )
00834 *
00835 *           Compute complex Givens rotation on right
00836 *           (Assume some element of C = (sA - wB) > unfl )
00837 *                            __
00838 *           (sA - wB) ( CZ   -SZ )
00839 *                     ( SZ    CZ )
00840 *
00841             C11R = S1*A11 - WR*B11
00842             C11I = -WI*B11
00843             C12 = S1*A12
00844             C21 = S1*A21
00845             C22R = S1*A22 - WR*B22
00846             C22I = -WI*B22
00847 *
00848             IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
00849      $          ABS( C22R )+ABS( C22I ) ) THEN
00850                T1 = SLAPY3( C12, C11R, C11I )
00851                CZ = C12 / T1
00852                SZR = -C11R / T1
00853                SZI = -C11I / T1
00854             ELSE
00855                CZ = SLAPY2( C22R, C22I )
00856                IF( CZ.LE.SAFMIN ) THEN
00857                   CZ = ZERO
00858                   SZR = ONE
00859                   SZI = ZERO
00860                ELSE
00861                   TEMPR = C22R / CZ
00862                   TEMPI = C22I / CZ
00863                   T1 = SLAPY2( CZ, C21 )
00864                   CZ = CZ / T1
00865                   SZR = -C21*TEMPR / T1
00866                   SZI = C21*TEMPI / T1
00867                END IF
00868             END IF
00869 *
00870 *           Compute Givens rotation on left
00871 *
00872 *           (  CQ   SQ )
00873 *           (  __      )  A or B
00874 *           ( -SQ   CQ )
00875 *
00876             AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
00877             BN = ABS( B11 ) + ABS( B22 )
00878             WABS = ABS( WR ) + ABS( WI )
00879             IF( S1*AN.GT.WABS*BN ) THEN
00880                CQ = CZ*B11
00881                SQR = SZR*B22
00882                SQI = -SZI*B22
00883             ELSE
00884                A1R = CZ*A11 + SZR*A12
00885                A1I = SZI*A12
00886                A2R = CZ*A21 + SZR*A22
00887                A2I = SZI*A22
00888                CQ = SLAPY2( A1R, A1I )
00889                IF( CQ.LE.SAFMIN ) THEN
00890                   CQ = ZERO
00891                   SQR = ONE
00892                   SQI = ZERO
00893                ELSE
00894                   TEMPR = A1R / CQ
00895                   TEMPI = A1I / CQ
00896                   SQR = TEMPR*A2R + TEMPI*A2I
00897                   SQI = TEMPI*A2R - TEMPR*A2I
00898                END IF
00899             END IF
00900             T1 = SLAPY3( CQ, SQR, SQI )
00901             CQ = CQ / T1
00902             SQR = SQR / T1
00903             SQI = SQI / T1
00904 *
00905 *           Compute diagonal elements of QBZ
00906 *
00907             TEMPR = SQR*SZR - SQI*SZI
00908             TEMPI = SQR*SZI + SQI*SZR
00909             B1R = CQ*CZ*B11 + TEMPR*B22
00910             B1I = TEMPI*B22
00911             B1A = SLAPY2( B1R, B1I )
00912             B2R = CQ*CZ*B22 + TEMPR*B11
00913             B2I = -TEMPI*B11
00914             B2A = SLAPY2( B2R, B2I )
00915 *
00916 *           Normalize so beta > 0, and Im( alpha1 ) > 0
00917 *
00918             BETA( ILAST-1 ) = B1A
00919             BETA( ILAST ) = B2A
00920             ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
00921             ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
00922             ALPHAR( ILAST ) = ( WR*B2A )*S1INV
00923             ALPHAI( ILAST ) = -( WI*B2A )*S1INV
00924 *
00925 *           Step 3: Go to next block -- exit if finished.
00926 *
00927             ILAST = IFIRST - 1
00928             IF( ILAST.LT.ILO )
00929      $         GO TO 380
00930 *
00931 *           Reset counters
00932 *
00933             IITER = 0
00934             ESHIFT = ZERO
00935             IF( .NOT.ILSCHR ) THEN
00936                ILASTM = ILAST
00937                IF( IFRSTM.GT.ILAST )
00938      $            IFRSTM = ILO
00939             END IF
00940             GO TO 350
00941          ELSE
00942 *
00943 *           Usual case: 3x3 or larger block, using Francis implicit
00944 *                       double-shift
00945 *
00946 *                                    2
00947 *           Eigenvalue equation is  w  - c w + d = 0,
00948 *
00949 *                                         -1 2        -1
00950 *           so compute 1st column of  (A B  )  - c A B   + d
00951 *           using the formula in QZIT (from EISPACK)
00952 *
00953 *           We assume that the block is at least 3x3
00954 *
00955             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
00956      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
00957             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
00958      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
00959             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
00960      $             ( BSCALE*T( ILAST, ILAST ) )
00961             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
00962      $             ( BSCALE*T( ILAST, ILAST ) )
00963             U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST )
00964             AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) /
00965      $              ( BSCALE*T( IFIRST, IFIRST ) )
00966             AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) /
00967      $              ( BSCALE*T( IFIRST, IFIRST ) )
00968             AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) /
00969      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
00970             AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) /
00971      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
00972             AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) /
00973      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
00974             U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 )
00975 *
00976             V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 +
00977      $               AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L
00978             V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )-
00979      $               ( AD22-AD11L )+AD21*U12 )*AD21L
00980             V( 3 ) = AD32L*AD21L
00981 *
00982             ISTART = IFIRST
00983 *
00984             CALL SLARFG( 3, V( 1 ), V( 2 ), 1, TAU )
00985             V( 1 ) = ONE
00986 *
00987 *           Sweep
00988 *
00989             DO 290 J = ISTART, ILAST - 2
00990 *
00991 *              All but last elements: use 3x3 Householder transforms.
00992 *
00993 *              Zero (j-1)st column of A
00994 *
00995                IF( J.GT.ISTART ) THEN
00996                   V( 1 ) = H( J, J-1 )
00997                   V( 2 ) = H( J+1, J-1 )
00998                   V( 3 ) = H( J+2, J-1 )
00999 *
01000                   CALL SLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU )
01001                   V( 1 ) = ONE
01002                   H( J+1, J-1 ) = ZERO
01003                   H( J+2, J-1 ) = ZERO
01004                END IF
01005 *
01006                DO 230 JC = J, ILASTM
01007                   TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )*
01008      $                   H( J+2, JC ) )
01009                   H( J, JC ) = H( J, JC ) - TEMP
01010                   H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 )
01011                   H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 )
01012                   TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )*
01013      $                    T( J+2, JC ) )
01014                   T( J, JC ) = T( J, JC ) - TEMP2
01015                   T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 )
01016                   T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 )
01017   230          CONTINUE
01018                IF( ILQ ) THEN
01019                   DO 240 JR = 1, N
01020                      TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )*
01021      $                      Q( JR, J+2 ) )
01022                      Q( JR, J ) = Q( JR, J ) - TEMP
01023                      Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 )
01024                      Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 )
01025   240             CONTINUE
01026                END IF
01027 *
01028 *              Zero j-th column of B (see SLAGBC for details)
01029 *
01030 *              Swap rows to pivot
01031 *
01032                ILPIVT = .FALSE.
01033                TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) )
01034                TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) )
01035                IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN
01036                   SCALE = ZERO
01037                   U1 = ONE
01038                   U2 = ZERO
01039                   GO TO 250
01040                ELSE IF( TEMP.GE.TEMP2 ) THEN
01041                   W11 = T( J+1, J+1 )
01042                   W21 = T( J+2, J+1 )
01043                   W12 = T( J+1, J+2 )
01044                   W22 = T( J+2, J+2 )
01045                   U1 = T( J+1, J )
01046                   U2 = T( J+2, J )
01047                ELSE
01048                   W21 = T( J+1, J+1 )
01049                   W11 = T( J+2, J+1 )
01050                   W22 = T( J+1, J+2 )
01051                   W12 = T( J+2, J+2 )
01052                   U2 = T( J+1, J )
01053                   U1 = T( J+2, J )
01054                END IF
01055 *
01056 *              Swap columns if nec.
01057 *
01058                IF( ABS( W12 ).GT.ABS( W11 ) ) THEN
01059                   ILPIVT = .TRUE.
01060                   TEMP = W12
01061                   TEMP2 = W22
01062                   W12 = W11
01063                   W22 = W21
01064                   W11 = TEMP
01065                   W21 = TEMP2
01066                END IF
01067 *
01068 *              LU-factor
01069 *
01070                TEMP = W21 / W11
01071                U2 = U2 - TEMP*U1
01072                W22 = W22 - TEMP*W12
01073                W21 = ZERO
01074 *
01075 *              Compute SCALE
01076 *
01077                SCALE = ONE
01078                IF( ABS( W22 ).LT.SAFMIN ) THEN
01079                   SCALE = ZERO
01080                   U2 = ONE
01081                   U1 = -W12 / W11
01082                   GO TO 250
01083                END IF
01084                IF( ABS( W22 ).LT.ABS( U2 ) )
01085      $            SCALE = ABS( W22 / U2 )
01086                IF( ABS( W11 ).LT.ABS( U1 ) )
01087      $            SCALE = MIN( SCALE, ABS( W11 / U1 ) )
01088 *
01089 *              Solve
01090 *
01091                U2 = ( SCALE*U2 ) / W22
01092                U1 = ( SCALE*U1-W12*U2 ) / W11
01093 *
01094   250          CONTINUE
01095                IF( ILPIVT ) THEN
01096                   TEMP = U2
01097                   U2 = U1
01098                   U1 = TEMP
01099                END IF
01100 *
01101 *              Compute Householder Vector
01102 *
01103                T1 = SQRT( SCALE**2+U1**2+U2**2 )
01104                TAU = ONE + SCALE / T1
01105                VS = -ONE / ( SCALE+T1 )
01106                V( 1 ) = ONE
01107                V( 2 ) = VS*U1
01108                V( 3 ) = VS*U2
01109 *
01110 *              Apply transformations from the right.
01111 *
01112                DO 260 JR = IFRSTM, MIN( J+3, ILAST )
01113                   TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )*
01114      $                   H( JR, J+2 ) )
01115                   H( JR, J ) = H( JR, J ) - TEMP
01116                   H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 )
01117                   H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 )
01118   260          CONTINUE
01119                DO 270 JR = IFRSTM, J + 2
01120                   TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )*
01121      $                   T( JR, J+2 ) )
01122                   T( JR, J ) = T( JR, J ) - TEMP
01123                   T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 )
01124                   T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 )
01125   270          CONTINUE
01126                IF( ILZ ) THEN
01127                   DO 280 JR = 1, N
01128                      TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )*
01129      $                      Z( JR, J+2 ) )
01130                      Z( JR, J ) = Z( JR, J ) - TEMP
01131                      Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 )
01132                      Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 )
01133   280             CONTINUE
01134                END IF
01135                T( J+1, J ) = ZERO
01136                T( J+2, J ) = ZERO
01137   290       CONTINUE
01138 *
01139 *           Last elements: Use Givens rotations
01140 *
01141 *           Rotations from the left
01142 *
01143             J = ILAST - 1
01144             TEMP = H( J, J-1 )
01145             CALL SLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
01146             H( J+1, J-1 ) = ZERO
01147 *
01148             DO 300 JC = J, ILASTM
01149                TEMP = C*H( J, JC ) + S*H( J+1, JC )
01150                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
01151                H( J, JC ) = TEMP
01152                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
01153                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
01154                T( J, JC ) = TEMP2
01155   300       CONTINUE
01156             IF( ILQ ) THEN
01157                DO 310 JR = 1, N
01158                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
01159                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
01160                   Q( JR, J ) = TEMP
01161   310          CONTINUE
01162             END IF
01163 *
01164 *           Rotations from the right.
01165 *
01166             TEMP = T( J+1, J+1 )
01167             CALL SLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
01168             T( J+1, J ) = ZERO
01169 *
01170             DO 320 JR = IFRSTM, ILAST
01171                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
01172                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
01173                H( JR, J+1 ) = TEMP
01174   320       CONTINUE
01175             DO 330 JR = IFRSTM, ILAST - 1
01176                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
01177                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
01178                T( JR, J+1 ) = TEMP
01179   330       CONTINUE
01180             IF( ILZ ) THEN
01181                DO 340 JR = 1, N
01182                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
01183                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
01184                   Z( JR, J+1 ) = TEMP
01185   340          CONTINUE
01186             END IF
01187 *
01188 *           End of Double-Shift code
01189 *
01190          END IF
01191 *
01192          GO TO 350
01193 *
01194 *        End of iteration loop
01195 *
01196   350    CONTINUE
01197   360 CONTINUE
01198 *
01199 *     Drop-through = non-convergence
01200 *
01201       INFO = ILAST
01202       GO TO 420
01203 *
01204 *     Successful completion of all QZ steps
01205 *
01206   380 CONTINUE
01207 *
01208 *     Set Eigenvalues 1:ILO-1
01209 *
01210       DO 410 J = 1, ILO - 1
01211          IF( T( J, J ).LT.ZERO ) THEN
01212             IF( ILSCHR ) THEN
01213                DO 390 JR = 1, J
01214                   H( JR, J ) = -H( JR, J )
01215                   T( JR, J ) = -T( JR, J )
01216   390          CONTINUE
01217             ELSE
01218                H( J, J ) = -H( J, J )
01219                T( J, J ) = -T( J, J )
01220             END IF
01221             IF( ILZ ) THEN
01222                DO 400 JR = 1, N
01223                   Z( JR, J ) = -Z( JR, J )
01224   400          CONTINUE
01225             END IF
01226          END IF
01227          ALPHAR( J ) = H( J, J )
01228          ALPHAI( J ) = ZERO
01229          BETA( J ) = T( J, J )
01230   410 CONTINUE
01231 *
01232 *     Normal Termination
01233 *
01234       INFO = 0
01235 *
01236 *     Exit (other than argument error) -- return optimal workspace size
01237 *
01238   420 CONTINUE
01239       WORK( 1 ) = REAL( N )
01240       RETURN
01241 *
01242 *     End of SHGEQZ
01243 *
01244       END