0001:       SUBROUTINE SHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
0002:      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
0003:      $                   LWORK, INFO )
0004: *
0005: *  -- LAPACK routine (version 3.2) --
0006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
0007: *     November 2006
0008: *
0009: *     .. Scalar Arguments ..
0010:       CHARACTER          COMPQ, COMPZ, JOB
0011:       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
0012: *     ..
0013: *     .. Array Arguments ..
0014:       REAL               ALPHAI( * ), ALPHAR( * ), BETA( * ),
0015:      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
0016:      $                   WORK( * ), Z( LDZ, * )
0017: *     ..
0018: *
0019: *  Purpose
0020: *  =======
0021: *
0022: *  SHGEQZ computes the eigenvalues of a real matrix pair (H,T),
0023: *  where H is an upper Hessenberg matrix and T is upper triangular,
0024: *  using the double-shift QZ method.
0025: *  Matrix pairs of this type are produced by the reduction to
0026: *  generalized upper Hessenberg form of a real matrix pair (A,B):
0027: *
0028: *     A = Q1*H*Z1**T,  B = Q1*T*Z1**T,
0029: *
0030: *  as computed by SGGHRD.
0031: *
0032: *  If JOB='S', then the Hessenberg-triangular pair (H,T) is
0033: *  also reduced to generalized Schur form,
0034: *  
0035: *     H = Q*S*Z**T,  T = Q*P*Z**T,
0036: *  
0037: *  where Q and Z are orthogonal matrices, P is an upper triangular
0038: *  matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
0039: *  diagonal blocks.
0040: *
0041: *  The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
0042: *  (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
0043: *  eigenvalues.
0044: *
0045: *  Additionally, the 2-by-2 upper triangular diagonal blocks of P
0046: *  corresponding to 2-by-2 blocks of S are reduced to positive diagonal
0047: *  form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
0048: *  P(j,j) > 0, and P(j+1,j+1) > 0.
0049: *
0050: *  Optionally, the orthogonal matrix Q from the generalized Schur
0051: *  factorization may be postmultiplied into an input matrix Q1, and the
0052: *  orthogonal matrix Z may be postmultiplied into an input matrix Z1.
0053: *  If Q1 and Z1 are the orthogonal matrices from SGGHRD that reduced
0054: *  the matrix pair (A,B) to generalized upper Hessenberg form, then the
0055: *  output matrices Q1*Q and Z1*Z are the orthogonal factors from the
0056: *  generalized Schur factorization of (A,B):
0057: *
0058: *     A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.
0059: *  
0060: *  To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
0061: *  of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
0062: *  complex and beta real.
0063: *  If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
0064: *  generalized nonsymmetric eigenvalue problem (GNEP)
0065: *     A*x = lambda*B*x
0066: *  and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
0067: *  alternate form of the GNEP
0068: *     mu*A*y = B*y.
0069: *  Real eigenvalues can be read directly from the generalized Schur
0070: *  form: 
0071: *    alpha = S(i,i), beta = P(i,i).
0072: *
0073: *  Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
0074: *       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
0075: *       pp. 241--256.
0076: *
0077: *  Arguments
0078: *  =========
0079: *
0080: *  JOB     (input) CHARACTER*1
0081: *          = 'E': Compute eigenvalues only;
0082: *          = 'S': Compute eigenvalues and the Schur form. 
0083: *
0084: *  COMPQ   (input) CHARACTER*1
0085: *          = 'N': Left Schur vectors (Q) are not computed;
0086: *          = 'I': Q is initialized to the unit matrix and the matrix Q
0087: *                 of left Schur vectors of (H,T) is returned;
0088: *          = 'V': Q must contain an orthogonal matrix Q1 on entry and
0089: *                 the product Q1*Q is returned.
0090: *
0091: *  COMPZ   (input) CHARACTER*1
0092: *          = 'N': Right Schur vectors (Z) are not computed;
0093: *          = 'I': Z is initialized to the unit matrix and the matrix Z
0094: *                 of right Schur vectors of (H,T) is returned;
0095: *          = 'V': Z must contain an orthogonal matrix Z1 on entry and
0096: *                 the product Z1*Z is returned.
0097: *
0098: *  N       (input) INTEGER
0099: *          The order of the matrices H, T, Q, and Z.  N >= 0.
0100: *
0101: *  ILO     (input) INTEGER
0102: *  IHI     (input) INTEGER
0103: *          ILO and IHI mark the rows and columns of H which are in
0104: *          Hessenberg form.  It is assumed that A is already upper
0105: *          triangular in rows and columns 1:ILO-1 and IHI+1:N.
0106: *          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
0107: *
0108: *  H       (input/output) REAL array, dimension (LDH, N)
0109: *          On entry, the N-by-N upper Hessenberg matrix H.
0110: *          On exit, if JOB = 'S', H contains the upper quasi-triangular
0111: *          matrix S from the generalized Schur factorization;
0112: *          2-by-2 diagonal blocks (corresponding to complex conjugate
0113: *          pairs of eigenvalues) are returned in standard form, with
0114: *          H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0.
0115: *          If JOB = 'E', the diagonal blocks of H match those of S, but
0116: *          the rest of H is unspecified.
0117: *
0118: *  LDH     (input) INTEGER
0119: *          The leading dimension of the array H.  LDH >= max( 1, N ).
0120: *
0121: *  T       (input/output) REAL array, dimension (LDT, N)
0122: *          On entry, the N-by-N upper triangular matrix T.
0123: *          On exit, if JOB = 'S', T contains the upper triangular
0124: *          matrix P from the generalized Schur factorization;
0125: *          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
0126: *          are reduced to positive diagonal form, i.e., if H(j+1,j) is
0127: *          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
0128: *          T(j+1,j+1) > 0.
0129: *          If JOB = 'E', the diagonal blocks of T match those of P, but
0130: *          the rest of T is unspecified.
0131: *
0132: *  LDT     (input) INTEGER
0133: *          The leading dimension of the array T.  LDT >= max( 1, N ).
0134: *
0135: *  ALPHAR  (output) REAL array, dimension (N)
0136: *          The real parts of each scalar alpha defining an eigenvalue
0137: *          of GNEP.
0138: *
0139: *  ALPHAI  (output) REAL array, dimension (N)
0140: *          The imaginary parts of each scalar alpha defining an
0141: *          eigenvalue of GNEP.
0142: *          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
0143: *          positive, then the j-th and (j+1)-st eigenvalues are a
0144: *          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
0145: *
0146: *  BETA    (output) REAL array, dimension (N)
0147: *          The scalars beta that define the eigenvalues of GNEP.
0148: *          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
0149: *          beta = BETA(j) represent the j-th eigenvalue of the matrix
0150: *          pair (A,B), in one of the forms lambda = alpha/beta or
0151: *          mu = beta/alpha.  Since either lambda or mu may overflow,
0152: *          they should not, in general, be computed.
0153: *
0154: *  Q       (input/output) REAL array, dimension (LDQ, N)
0155: *          On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in
0156: *          the reduction of (A,B) to generalized Hessenberg form.
0157: *          On exit, if COMPZ = 'I', the orthogonal matrix of left Schur
0158: *          vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix
0159: *          of left Schur vectors of (A,B).
0160: *          Not referenced if COMPZ = 'N'.
0161: *
0162: *  LDQ     (input) INTEGER
0163: *          The leading dimension of the array Q.  LDQ >= 1.
0164: *          If COMPQ='V' or 'I', then LDQ >= N.
0165: *
0166: *  Z       (input/output) REAL array, dimension (LDZ, N)
0167: *          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
0168: *          the reduction of (A,B) to generalized Hessenberg form.
0169: *          On exit, if COMPZ = 'I', the orthogonal matrix of
0170: *          right Schur vectors of (H,T), and if COMPZ = 'V', the
0171: *          orthogonal matrix of right Schur vectors of (A,B).
0172: *          Not referenced if COMPZ = 'N'.
0173: *
0174: *  LDZ     (input) INTEGER
0175: *          The leading dimension of the array Z.  LDZ >= 1.
0176: *          If COMPZ='V' or 'I', then LDZ >= N.
0177: *
0178: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
0179: *          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
0180: *
0181: *  LWORK   (input) INTEGER
0182: *          The dimension of the array WORK.  LWORK >= max(1,N).
0183: *
0184: *          If LWORK = -1, then a workspace query is assumed; the routine
0185: *          only calculates the optimal size of the WORK array, returns
0186: *          this value as the first entry of the WORK array, and no error
0187: *          message related to LWORK is issued by XERBLA.
0188: *
0189: *  INFO    (output) INTEGER
0190: *          = 0: successful exit
0191: *          < 0: if INFO = -i, the i-th argument had an illegal value
0192: *          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
0193: *                     in Schur form, but ALPHAR(i), ALPHAI(i), and
0194: *                     BETA(i), i=INFO+1,...,N should be correct.
0195: *          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
0196: *                     in Schur form, but ALPHAR(i), ALPHAI(i), and
0197: *                     BETA(i), i=INFO-N+1,...,N should be correct.
0198: *
0199: *  Further Details
0200: *  ===============
0201: *
0202: *  Iteration counters:
0203: *
0204: *  JITER  -- counts iterations.
0205: *  IITER  -- counts iterations run since ILAST was last
0206: *            changed.  This is therefore reset only when a 1-by-1 or
0207: *            2-by-2 block deflates off the bottom.
0208: *
0209: *  =====================================================================
0210: *
0211: *     .. Parameters ..
0212: *    $                     SAFETY = 1.0E+0 )
0213:       REAL               HALF, ZERO, ONE, SAFETY
0214:       PARAMETER          ( HALF = 0.5E+0, ZERO = 0.0E+0, ONE = 1.0E+0,
0215:      $                   SAFETY = 1.0E+2 )
0216: *     ..
0217: *     .. Local Scalars ..
0218:       LOGICAL            ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
0219:      $                   LQUERY
0220:       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
0221:      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
0222:      $                   JR, MAXIT
0223:       REAL               A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
0224:      $                   AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
0225:      $                   AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
0226:      $                   B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
0227:      $                   BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
0228:      $                   CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
0229:      $                   SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
0230:      $                   TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
0231:      $                   U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
0232:      $                   WR2
0233: *     ..
0234: *     .. Local Arrays ..
0235:       REAL               V( 3 )
0236: *     ..
0237: *     .. External Functions ..
0238:       LOGICAL            LSAME
0239:       REAL               SLAMCH, SLANHS, SLAPY2, SLAPY3
0240:       EXTERNAL           LSAME, SLAMCH, SLANHS, SLAPY2, SLAPY3
0241: *     ..
0242: *     .. External Subroutines ..
0243:       EXTERNAL           SLAG2, SLARFG, SLARTG, SLASET, SLASV2, SROT,
0244:      $                   XERBLA
0245: *     ..
0246: *     .. Intrinsic Functions ..
0247:       INTRINSIC          ABS, MAX, MIN, REAL, SQRT
0248: *     ..
0249: *     .. Executable Statements ..
0250: *
0251: *     Decode JOB, COMPQ, COMPZ
0252: *
0253:       IF( LSAME( JOB, 'E' ) ) THEN
0254:          ILSCHR = .FALSE.
0255:          ISCHUR = 1
0256:       ELSE IF( LSAME( JOB, 'S' ) ) THEN
0257:          ILSCHR = .TRUE.
0258:          ISCHUR = 2
0259:       ELSE
0260:          ISCHUR = 0
0261:       END IF
0262: *
0263:       IF( LSAME( COMPQ, 'N' ) ) THEN
0264:          ILQ = .FALSE.
0265:          ICOMPQ = 1
0266:       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
0267:          ILQ = .TRUE.
0268:          ICOMPQ = 2
0269:       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
0270:          ILQ = .TRUE.
0271:          ICOMPQ = 3
0272:       ELSE
0273:          ICOMPQ = 0
0274:       END IF
0275: *
0276:       IF( LSAME( COMPZ, 'N' ) ) THEN
0277:          ILZ = .FALSE.
0278:          ICOMPZ = 1
0279:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
0280:          ILZ = .TRUE.
0281:          ICOMPZ = 2
0282:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
0283:          ILZ = .TRUE.
0284:          ICOMPZ = 3
0285:       ELSE
0286:          ICOMPZ = 0
0287:       END IF
0288: *
0289: *     Check Argument Values
0290: *
0291:       INFO = 0
0292:       WORK( 1 ) = MAX( 1, N )
0293:       LQUERY = ( LWORK.EQ.-1 )
0294:       IF( ISCHUR.EQ.0 ) THEN
0295:          INFO = -1
0296:       ELSE IF( ICOMPQ.EQ.0 ) THEN
0297:          INFO = -2
0298:       ELSE IF( ICOMPZ.EQ.0 ) THEN
0299:          INFO = -3
0300:       ELSE IF( N.LT.0 ) THEN
0301:          INFO = -4
0302:       ELSE IF( ILO.LT.1 ) THEN
0303:          INFO = -5
0304:       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
0305:          INFO = -6
0306:       ELSE IF( LDH.LT.N ) THEN
0307:          INFO = -8
0308:       ELSE IF( LDT.LT.N ) THEN
0309:          INFO = -10
0310:       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
0311:          INFO = -15
0312:       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
0313:          INFO = -17
0314:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
0315:          INFO = -19
0316:       END IF
0317:       IF( INFO.NE.0 ) THEN
0318:          CALL XERBLA( 'SHGEQZ', -INFO )
0319:          RETURN
0320:       ELSE IF( LQUERY ) THEN
0321:          RETURN
0322:       END IF
0323: *
0324: *     Quick return if possible
0325: *
0326:       IF( N.LE.0 ) THEN
0327:          WORK( 1 ) = REAL( 1 )
0328:          RETURN
0329:       END IF
0330: *
0331: *     Initialize Q and Z
0332: *
0333:       IF( ICOMPQ.EQ.3 )
0334:      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
0335:       IF( ICOMPZ.EQ.3 )
0336:      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
0337: *
0338: *     Machine Constants
0339: *
0340:       IN = IHI + 1 - ILO
0341:       SAFMIN = SLAMCH( 'S' )
0342:       SAFMAX = ONE / SAFMIN
0343:       ULP = SLAMCH( 'E' )*SLAMCH( 'B' )
0344:       ANORM = SLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK )
0345:       BNORM = SLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK )
0346:       ATOL = MAX( SAFMIN, ULP*ANORM )
0347:       BTOL = MAX( SAFMIN, ULP*BNORM )
0348:       ASCALE = ONE / MAX( SAFMIN, ANORM )
0349:       BSCALE = ONE / MAX( SAFMIN, BNORM )
0350: *
0351: *     Set Eigenvalues IHI+1:N
0352: *
0353:       DO 30 J = IHI + 1, N
0354:          IF( T( J, J ).LT.ZERO ) THEN
0355:             IF( ILSCHR ) THEN
0356:                DO 10 JR = 1, J
0357:                   H( JR, J ) = -H( JR, J )
0358:                   T( JR, J ) = -T( JR, J )
0359:    10          CONTINUE
0360:             ELSE
0361:                H( J, J ) = -H( J, J )
0362:                T( J, J ) = -T( J, J )
0363:             END IF
0364:             IF( ILZ ) THEN
0365:                DO 20 JR = 1, N
0366:                   Z( JR, J ) = -Z( JR, J )
0367:    20          CONTINUE
0368:             END IF
0369:          END IF
0370:          ALPHAR( J ) = H( J, J )
0371:          ALPHAI( J ) = ZERO
0372:          BETA( J ) = T( J, J )
0373:    30 CONTINUE
0374: *
0375: *     If IHI < ILO, skip QZ steps
0376: *
0377:       IF( IHI.LT.ILO )
0378:      $   GO TO 380
0379: *
0380: *     MAIN QZ ITERATION LOOP
0381: *
0382: *     Initialize dynamic indices
0383: *
0384: *     Eigenvalues ILAST+1:N have been found.
0385: *        Column operations modify rows IFRSTM:whatever.
0386: *        Row operations modify columns whatever:ILASTM.
0387: *
0388: *     If only eigenvalues are being computed, then
0389: *        IFRSTM is the row of the last splitting row above row ILAST;
0390: *        this is always at least ILO.
0391: *     IITER counts iterations since the last eigenvalue was found,
0392: *        to tell when to use an extraordinary shift.
0393: *     MAXIT is the maximum number of QZ sweeps allowed.
0394: *
0395:       ILAST = IHI
0396:       IF( ILSCHR ) THEN
0397:          IFRSTM = 1
0398:          ILASTM = N
0399:       ELSE
0400:          IFRSTM = ILO
0401:          ILASTM = IHI
0402:       END IF
0403:       IITER = 0
0404:       ESHIFT = ZERO
0405:       MAXIT = 30*( IHI-ILO+1 )
0406: *
0407:       DO 360 JITER = 1, MAXIT
0408: *
0409: *        Split the matrix if possible.
0410: *
0411: *        Two tests:
0412: *           1: H(j,j-1)=0  or  j=ILO
0413: *           2: T(j,j)=0
0414: *
0415:          IF( ILAST.EQ.ILO ) THEN
0416: *
0417: *           Special case: j=ILAST
0418: *
0419:             GO TO 80
0420:          ELSE
0421:             IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
0422:                H( ILAST, ILAST-1 ) = ZERO
0423:                GO TO 80
0424:             END IF
0425:          END IF
0426: *
0427:          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
0428:             T( ILAST, ILAST ) = ZERO
0429:             GO TO 70
0430:          END IF
0431: *
0432: *        General case: j<ILAST
0433: *
0434:          DO 60 J = ILAST - 1, ILO, -1
0435: *
0436: *           Test 1: for H(j,j-1)=0 or j=ILO
0437: *
0438:             IF( J.EQ.ILO ) THEN
0439:                ILAZRO = .TRUE.
0440:             ELSE
0441:                IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN
0442:                   H( J, J-1 ) = ZERO
0443:                   ILAZRO = .TRUE.
0444:                ELSE
0445:                   ILAZRO = .FALSE.
0446:                END IF
0447:             END IF
0448: *
0449: *           Test 2: for T(j,j)=0
0450: *
0451:             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
0452:                T( J, J ) = ZERO
0453: *
0454: *              Test 1a: Check for 2 consecutive small subdiagonals in A
0455: *
0456:                ILAZR2 = .FALSE.
0457:                IF( .NOT.ILAZRO ) THEN
0458:                   TEMP = ABS( H( J, J-1 ) )
0459:                   TEMP2 = ABS( H( J, J ) )
0460:                   TEMPR = MAX( TEMP, TEMP2 )
0461:                   IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
0462:                      TEMP = TEMP / TEMPR
0463:                      TEMP2 = TEMP2 / TEMPR
0464:                   END IF
0465:                   IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2*
0466:      $                ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
0467:                END IF
0468: *
0469: *              If both tests pass (1 & 2), i.e., the leading diagonal
0470: *              element of B in the block is zero, split a 1x1 block off
0471: *              at the top. (I.e., at the J-th row/column) The leading
0472: *              diagonal element of the remainder can also be zero, so
0473: *              this may have to be done repeatedly.
0474: *
0475:                IF( ILAZRO .OR. ILAZR2 ) THEN
0476:                   DO 40 JCH = J, ILAST - 1
0477:                      TEMP = H( JCH, JCH )
0478:                      CALL SLARTG( TEMP, H( JCH+1, JCH ), C, S,
0479:      $                            H( JCH, JCH ) )
0480:                      H( JCH+1, JCH ) = ZERO
0481:                      CALL SROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
0482:      $                          H( JCH+1, JCH+1 ), LDH, C, S )
0483:                      CALL SROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
0484:      $                          T( JCH+1, JCH+1 ), LDT, C, S )
0485:                      IF( ILQ )
0486:      $                  CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
0487:      $                             C, S )
0488:                      IF( ILAZR2 )
0489:      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
0490:                      ILAZR2 = .FALSE.
0491:                      IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
0492:                         IF( JCH+1.GE.ILAST ) THEN
0493:                            GO TO 80
0494:                         ELSE
0495:                            IFIRST = JCH + 1
0496:                            GO TO 110
0497:                         END IF
0498:                      END IF
0499:                      T( JCH+1, JCH+1 ) = ZERO
0500:    40             CONTINUE
0501:                   GO TO 70
0502:                ELSE
0503: *
0504: *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
0505: *                 Then process as in the case T(ILAST,ILAST)=0
0506: *
0507:                   DO 50 JCH = J, ILAST - 1
0508:                      TEMP = T( JCH, JCH+1 )
0509:                      CALL SLARTG( TEMP, T( JCH+1, JCH+1 ), C, S,
0510:      $                            T( JCH, JCH+1 ) )
0511:                      T( JCH+1, JCH+1 ) = ZERO
0512:                      IF( JCH.LT.ILASTM-1 )
0513:      $                  CALL SROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
0514:      $                             T( JCH+1, JCH+2 ), LDT, C, S )
0515:                      CALL SROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
0516:      $                          H( JCH+1, JCH-1 ), LDH, C, S )
0517:                      IF( ILQ )
0518:      $                  CALL SROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
0519:      $                             C, S )
0520:                      TEMP = H( JCH+1, JCH )
0521:                      CALL SLARTG( TEMP, H( JCH+1, JCH-1 ), C, S,
0522:      $                            H( JCH+1, JCH ) )
0523:                      H( JCH+1, JCH-1 ) = ZERO
0524:                      CALL SROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
0525:      $                          H( IFRSTM, JCH-1 ), 1, C, S )
0526:                      CALL SROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
0527:      $                          T( IFRSTM, JCH-1 ), 1, C, S )
0528:                      IF( ILZ )
0529:      $                  CALL SROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
0530:      $                             C, S )
0531:    50             CONTINUE
0532:                   GO TO 70
0533:                END IF
0534:             ELSE IF( ILAZRO ) THEN
0535: *
0536: *              Only test 1 passed -- work on J:ILAST
0537: *
0538:                IFIRST = J
0539:                GO TO 110
0540:             END IF
0541: *
0542: *           Neither test passed -- try next J
0543: *
0544:    60    CONTINUE
0545: *
0546: *        (Drop-through is "impossible")
0547: *
0548:          INFO = N + 1
0549:          GO TO 420
0550: *
0551: *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
0552: *        1x1 block.
0553: *
0554:    70    CONTINUE
0555:          TEMP = H( ILAST, ILAST )
0556:          CALL SLARTG( TEMP, H( ILAST, ILAST-1 ), C, S,
0557:      $                H( ILAST, ILAST ) )
0558:          H( ILAST, ILAST-1 ) = ZERO
0559:          CALL SROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
0560:      $              H( IFRSTM, ILAST-1 ), 1, C, S )
0561:          CALL SROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
0562:      $              T( IFRSTM, ILAST-1 ), 1, C, S )
0563:          IF( ILZ )
0564:      $      CALL SROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
0565: *
0566: *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
0567: *                              and BETA
0568: *
0569:    80    CONTINUE
0570:          IF( T( ILAST, ILAST ).LT.ZERO ) THEN
0571:             IF( ILSCHR ) THEN
0572:                DO 90 J = IFRSTM, ILAST
0573:                   H( J, ILAST ) = -H( J, ILAST )
0574:                   T( J, ILAST ) = -T( J, ILAST )
0575:    90          CONTINUE
0576:             ELSE
0577:                H( ILAST, ILAST ) = -H( ILAST, ILAST )
0578:                T( ILAST, ILAST ) = -T( ILAST, ILAST )
0579:             END IF
0580:             IF( ILZ ) THEN
0581:                DO 100 J = 1, N
0582:                   Z( J, ILAST ) = -Z( J, ILAST )
0583:   100          CONTINUE
0584:             END IF
0585:          END IF
0586:          ALPHAR( ILAST ) = H( ILAST, ILAST )
0587:          ALPHAI( ILAST ) = ZERO
0588:          BETA( ILAST ) = T( ILAST, ILAST )
0589: *
0590: *        Go to next block -- exit if finished.
0591: *
0592:          ILAST = ILAST - 1
0593:          IF( ILAST.LT.ILO )
0594:      $      GO TO 380
0595: *
0596: *        Reset counters
0597: *
0598:          IITER = 0
0599:          ESHIFT = ZERO
0600:          IF( .NOT.ILSCHR ) THEN
0601:             ILASTM = ILAST
0602:             IF( IFRSTM.GT.ILAST )
0603:      $         IFRSTM = ILO
0604:          END IF
0605:          GO TO 350
0606: *
0607: *        QZ step
0608: *
0609: *        This iteration only involves rows/columns IFIRST:ILAST. We
0610: *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
0611: *
0612:   110    CONTINUE
0613:          IITER = IITER + 1
0614:          IF( .NOT.ILSCHR ) THEN
0615:             IFRSTM = IFIRST
0616:          END IF
0617: *
0618: *        Compute single shifts.
0619: *
0620: *        At this point, IFIRST < ILAST, and the diagonal elements of
0621: *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
0622: *        magnitude)
0623: *
0624:          IF( ( IITER / 10 )*10.EQ.IITER ) THEN
0625: *
0626: *           Exceptional shift.  Chosen for no particularly good reason.
0627: *           (Single shift only.)
0628: *
0629:             IF( ( REAL( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
0630:      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
0631:                ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /
0632:      $                  T( ILAST-1, ILAST-1 )
0633:             ELSE
0634:                ESHIFT = ESHIFT + ONE / ( SAFMIN*REAL( MAXIT ) )
0635:             END IF
0636:             S1 = ONE
0637:             WR = ESHIFT
0638: *
0639:          ELSE
0640: *
0641: *           Shifts based on the generalized eigenvalues of the
0642: *           bottom-right 2x2 block of A and B. The first eigenvalue
0643: *           returned by SLAG2 is the Wilkinson shift (AEP p.512),
0644: *
0645:             CALL SLAG2( H( ILAST-1, ILAST-1 ), LDH,
0646:      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
0647:      $                  S2, WR, WR2, WI )
0648: *
0649:             TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
0650:             IF( WI.NE.ZERO )
0651:      $         GO TO 200
0652:          END IF
0653: *
0654: *        Fiddle with shift to avoid overflow
0655: *
0656:          TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
0657:          IF( S1.GT.TEMP ) THEN
0658:             SCALE = TEMP / S1
0659:          ELSE
0660:             SCALE = ONE
0661:          END IF
0662: *
0663:          TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
0664:          IF( ABS( WR ).GT.TEMP )
0665:      $      SCALE = MIN( SCALE, TEMP / ABS( WR ) )
0666:          S1 = SCALE*S1
0667:          WR = SCALE*WR
0668: *
0669: *        Now check for two consecutive small subdiagonals.
0670: *
0671:          DO 120 J = ILAST - 1, IFIRST + 1, -1
0672:             ISTART = J
0673:             TEMP = ABS( S1*H( J, J-1 ) )
0674:             TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) )
0675:             TEMPR = MAX( TEMP, TEMP2 )
0676:             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
0677:                TEMP = TEMP / TEMPR
0678:                TEMP2 = TEMP2 / TEMPR
0679:             END IF
0680:             IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
0681:      $          TEMP2 )GO TO 130
0682:   120    CONTINUE
0683: *
0684:          ISTART = IFIRST
0685:   130    CONTINUE
0686: *
0687: *        Do an implicit single-shift QZ sweep.
0688: *
0689: *        Initial Q
0690: *
0691:          TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART )
0692:          TEMP2 = S1*H( ISTART+1, ISTART )
0693:          CALL SLARTG( TEMP, TEMP2, C, S, TEMPR )
0694: *
0695: *        Sweep
0696: *
0697:          DO 190 J = ISTART, ILAST - 1
0698:             IF( J.GT.ISTART ) THEN
0699:                TEMP = H( J, J-1 )
0700:                CALL SLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
0701:                H( J+1, J-1 ) = ZERO
0702:             END IF
0703: *
0704:             DO 140 JC = J, ILASTM
0705:                TEMP = C*H( J, JC ) + S*H( J+1, JC )
0706:                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
0707:                H( J, JC ) = TEMP
0708:                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
0709:                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
0710:                T( J, JC ) = TEMP2
0711:   140       CONTINUE
0712:             IF( ILQ ) THEN
0713:                DO 150 JR = 1, N
0714:                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
0715:                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
0716:                   Q( JR, J ) = TEMP
0717:   150          CONTINUE
0718:             END IF
0719: *
0720:             TEMP = T( J+1, J+1 )
0721:             CALL SLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
0722:             T( J+1, J ) = ZERO
0723: *
0724:             DO 160 JR = IFRSTM, MIN( J+2, ILAST )
0725:                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
0726:                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
0727:                H( JR, J+1 ) = TEMP
0728:   160       CONTINUE
0729:             DO 170 JR = IFRSTM, J
0730:                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
0731:                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
0732:                T( JR, J+1 ) = TEMP
0733:   170       CONTINUE
0734:             IF( ILZ ) THEN
0735:                DO 180 JR = 1, N
0736:                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
0737:                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
0738:                   Z( JR, J+1 ) = TEMP
0739:   180          CONTINUE
0740:             END IF
0741:   190    CONTINUE
0742: *
0743:          GO TO 350
0744: *
0745: *        Use Francis double-shift
0746: *
0747: *        Note: the Francis double-shift should work with real shifts,
0748: *              but only if the block is at least 3x3.
0749: *              This code may break if this point is reached with
0750: *              a 2x2 block with real eigenvalues.
0751: *
0752:   200    CONTINUE
0753:          IF( IFIRST+1.EQ.ILAST ) THEN
0754: *
0755: *           Special case -- 2x2 block with complex eigenvectors
0756: *
0757: *           Step 1: Standardize, that is, rotate so that
0758: *
0759: *                       ( B11  0  )
0760: *                   B = (         )  with B11 non-negative.
0761: *                       (  0  B22 )
0762: *
0763:             CALL SLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ),
0764:      $                   T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
0765: *
0766:             IF( B11.LT.ZERO ) THEN
0767:                CR = -CR
0768:                SR = -SR
0769:                B11 = -B11
0770:                B22 = -B22
0771:             END IF
0772: *
0773:             CALL SROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,
0774:      $                 H( ILAST, ILAST-1 ), LDH, CL, SL )
0775:             CALL SROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
0776:      $                 H( IFRSTM, ILAST ), 1, CR, SR )
0777: *
0778:             IF( ILAST.LT.ILASTM )
0779:      $         CALL SROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
0780:      $                    T( ILAST, ILAST+1 ), LDH, CL, SL )
0781:             IF( IFRSTM.LT.ILAST-1 )
0782:      $         CALL SROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
0783:      $                    T( IFRSTM, ILAST ), 1, CR, SR )
0784: *
0785:             IF( ILQ )
0786:      $         CALL SROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
0787:      $                    SL )
0788:             IF( ILZ )
0789:      $         CALL SROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
0790:      $                    SR )
0791: *
0792:             T( ILAST-1, ILAST-1 ) = B11
0793:             T( ILAST-1, ILAST ) = ZERO
0794:             T( ILAST, ILAST-1 ) = ZERO
0795:             T( ILAST, ILAST ) = B22
0796: *
0797: *           If B22 is negative, negate column ILAST
0798: *
0799:             IF( B22.LT.ZERO ) THEN
0800:                DO 210 J = IFRSTM, ILAST
0801:                   H( J, ILAST ) = -H( J, ILAST )
0802:                   T( J, ILAST ) = -T( J, ILAST )
0803:   210          CONTINUE
0804: *
0805:                IF( ILZ ) THEN
0806:                   DO 220 J = 1, N
0807:                      Z( J, ILAST ) = -Z( J, ILAST )
0808:   220             CONTINUE
0809:                END IF
0810:             END IF
0811: *
0812: *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
0813: *
0814: *           Recompute shift
0815: *
0816:             CALL SLAG2( H( ILAST-1, ILAST-1 ), LDH,
0817:      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
0818:      $                  TEMP, WR, TEMP2, WI )
0819: *
0820: *           If standardization has perturbed the shift onto real line,
0821: *           do another (real single-shift) QR step.
0822: *
0823:             IF( WI.EQ.ZERO )
0824:      $         GO TO 350
0825:             S1INV = ONE / S1
0826: *
0827: *           Do EISPACK (QZVAL) computation of alpha and beta
0828: *
0829:             A11 = H( ILAST-1, ILAST-1 )
0830:             A21 = H( ILAST, ILAST-1 )
0831:             A12 = H( ILAST-1, ILAST )
0832:             A22 = H( ILAST, ILAST )
0833: *
0834: *           Compute complex Givens rotation on right
0835: *           (Assume some element of C = (sA - wB) > unfl )
0836: *                            __
0837: *           (sA - wB) ( CZ   -SZ )
0838: *                     ( SZ    CZ )
0839: *
0840:             C11R = S1*A11 - WR*B11
0841:             C11I = -WI*B11
0842:             C12 = S1*A12
0843:             C21 = S1*A21
0844:             C22R = S1*A22 - WR*B22
0845:             C22I = -WI*B22
0846: *
0847:             IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
0848:      $          ABS( C22R )+ABS( C22I ) ) THEN
0849:                T1 = SLAPY3( C12, C11R, C11I )
0850:                CZ = C12 / T1
0851:                SZR = -C11R / T1
0852:                SZI = -C11I / T1
0853:             ELSE
0854:                CZ = SLAPY2( C22R, C22I )
0855:                IF( CZ.LE.SAFMIN ) THEN
0856:                   CZ = ZERO
0857:                   SZR = ONE
0858:                   SZI = ZERO
0859:                ELSE
0860:                   TEMPR = C22R / CZ
0861:                   TEMPI = C22I / CZ
0862:                   T1 = SLAPY2( CZ, C21 )
0863:                   CZ = CZ / T1
0864:                   SZR = -C21*TEMPR / T1
0865:                   SZI = C21*TEMPI / T1
0866:                END IF
0867:             END IF
0868: *
0869: *           Compute Givens rotation on left
0870: *
0871: *           (  CQ   SQ )
0872: *           (  __      )  A or B
0873: *           ( -SQ   CQ )
0874: *
0875:             AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
0876:             BN = ABS( B11 ) + ABS( B22 )
0877:             WABS = ABS( WR ) + ABS( WI )
0878:             IF( S1*AN.GT.WABS*BN ) THEN
0879:                CQ = CZ*B11
0880:                SQR = SZR*B22
0881:                SQI = -SZI*B22
0882:             ELSE
0883:                A1R = CZ*A11 + SZR*A12
0884:                A1I = SZI*A12
0885:                A2R = CZ*A21 + SZR*A22
0886:                A2I = SZI*A22
0887:                CQ = SLAPY2( A1R, A1I )
0888:                IF( CQ.LE.SAFMIN ) THEN
0889:                   CQ = ZERO
0890:                   SQR = ONE
0891:                   SQI = ZERO
0892:                ELSE
0893:                   TEMPR = A1R / CQ
0894:                   TEMPI = A1I / CQ
0895:                   SQR = TEMPR*A2R + TEMPI*A2I
0896:                   SQI = TEMPI*A2R - TEMPR*A2I
0897:                END IF
0898:             END IF
0899:             T1 = SLAPY3( CQ, SQR, SQI )
0900:             CQ = CQ / T1
0901:             SQR = SQR / T1
0902:             SQI = SQI / T1
0903: *
0904: *           Compute diagonal elements of QBZ
0905: *
0906:             TEMPR = SQR*SZR - SQI*SZI
0907:             TEMPI = SQR*SZI + SQI*SZR
0908:             B1R = CQ*CZ*B11 + TEMPR*B22
0909:             B1I = TEMPI*B22
0910:             B1A = SLAPY2( B1R, B1I )
0911:             B2R = CQ*CZ*B22 + TEMPR*B11
0912:             B2I = -TEMPI*B11
0913:             B2A = SLAPY2( B2R, B2I )
0914: *
0915: *           Normalize so beta > 0, and Im( alpha1 ) > 0
0916: *
0917:             BETA( ILAST-1 ) = B1A
0918:             BETA( ILAST ) = B2A
0919:             ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
0920:             ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
0921:             ALPHAR( ILAST ) = ( WR*B2A )*S1INV
0922:             ALPHAI( ILAST ) = -( WI*B2A )*S1INV
0923: *
0924: *           Step 3: Go to next block -- exit if finished.
0925: *
0926:             ILAST = IFIRST - 1
0927:             IF( ILAST.LT.ILO )
0928:      $         GO TO 380
0929: *
0930: *           Reset counters
0931: *
0932:             IITER = 0
0933:             ESHIFT = ZERO
0934:             IF( .NOT.ILSCHR ) THEN
0935:                ILASTM = ILAST
0936:                IF( IFRSTM.GT.ILAST )
0937:      $            IFRSTM = ILO
0938:             END IF
0939:             GO TO 350
0940:          ELSE
0941: *
0942: *           Usual case: 3x3 or larger block, using Francis implicit
0943: *                       double-shift
0944: *
0945: *                                    2
0946: *           Eigenvalue equation is  w  - c w + d = 0,
0947: *
0948: *                                         -1 2        -1
0949: *           so compute 1st column of  (A B  )  - c A B   + d
0950: *           using the formula in QZIT (from EISPACK)
0951: *
0952: *           We assume that the block is at least 3x3
0953: *
0954:             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
0955:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
0956:             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
0957:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
0958:             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
0959:      $             ( BSCALE*T( ILAST, ILAST ) )
0960:             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
0961:      $             ( BSCALE*T( ILAST, ILAST ) )
0962:             U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST )
0963:             AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) /
0964:      $              ( BSCALE*T( IFIRST, IFIRST ) )
0965:             AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) /
0966:      $              ( BSCALE*T( IFIRST, IFIRST ) )
0967:             AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) /
0968:      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
0969:             AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) /
0970:      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
0971:             AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) /
0972:      $              ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
0973:             U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 )
0974: *
0975:             V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 +
0976:      $               AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L
0977:             V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )-
0978:      $               ( AD22-AD11L )+AD21*U12 )*AD21L
0979:             V( 3 ) = AD32L*AD21L
0980: *
0981:             ISTART = IFIRST
0982: *
0983:             CALL SLARFG( 3, V( 1 ), V( 2 ), 1, TAU )
0984:             V( 1 ) = ONE
0985: *
0986: *           Sweep
0987: *
0988:             DO 290 J = ISTART, ILAST - 2
0989: *
0990: *              All but last elements: use 3x3 Householder transforms.
0991: *
0992: *              Zero (j-1)st column of A
0993: *
0994:                IF( J.GT.ISTART ) THEN
0995:                   V( 1 ) = H( J, J-1 )
0996:                   V( 2 ) = H( J+1, J-1 )
0997:                   V( 3 ) = H( J+2, J-1 )
0998: *
0999:                   CALL SLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU )
1000:                   V( 1 ) = ONE
1001:                   H( J+1, J-1 ) = ZERO
1002:                   H( J+2, J-1 ) = ZERO
1003:                END IF
1004: *
1005:                DO 230 JC = J, ILASTM
1006:                   TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )*
1007:      $                   H( J+2, JC ) )
1008:                   H( J, JC ) = H( J, JC ) - TEMP
1009:                   H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 )
1010:                   H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 )
1011:                   TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )*
1012:      $                    T( J+2, JC ) )
1013:                   T( J, JC ) = T( J, JC ) - TEMP2
1014:                   T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 )
1015:                   T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 )
1016:   230          CONTINUE
1017:                IF( ILQ ) THEN
1018:                   DO 240 JR = 1, N
1019:                      TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )*
1020:      $                      Q( JR, J+2 ) )
1021:                      Q( JR, J ) = Q( JR, J ) - TEMP
1022:                      Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 )
1023:                      Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 )
1024:   240             CONTINUE
1025:                END IF
1026: *
1027: *              Zero j-th column of B (see SLAGBC for details)
1028: *
1029: *              Swap rows to pivot
1030: *
1031:                ILPIVT = .FALSE.
1032:                TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) )
1033:                TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) )
1034:                IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN
1035:                   SCALE = ZERO
1036:                   U1 = ONE
1037:                   U2 = ZERO
1038:                   GO TO 250
1039:                ELSE IF( TEMP.GE.TEMP2 ) THEN
1040:                   W11 = T( J+1, J+1 )
1041:                   W21 = T( J+2, J+1 )
1042:                   W12 = T( J+1, J+2 )
1043:                   W22 = T( J+2, J+2 )
1044:                   U1 = T( J+1, J )
1045:                   U2 = T( J+2, J )
1046:                ELSE
1047:                   W21 = T( J+1, J+1 )
1048:                   W11 = T( J+2, J+1 )
1049:                   W22 = T( J+1, J+2 )
1050:                   W12 = T( J+2, J+2 )
1051:                   U2 = T( J+1, J )
1052:                   U1 = T( J+2, J )
1053:                END IF
1054: *
1055: *              Swap columns if nec.
1056: *
1057:                IF( ABS( W12 ).GT.ABS( W11 ) ) THEN
1058:                   ILPIVT = .TRUE.
1059:                   TEMP = W12
1060:                   TEMP2 = W22
1061:                   W12 = W11
1062:                   W22 = W21
1063:                   W11 = TEMP
1064:                   W21 = TEMP2
1065:                END IF
1066: *
1067: *              LU-factor
1068: *
1069:                TEMP = W21 / W11
1070:                U2 = U2 - TEMP*U1
1071:                W22 = W22 - TEMP*W12
1072:                W21 = ZERO
1073: *
1074: *              Compute SCALE
1075: *
1076:                SCALE = ONE
1077:                IF( ABS( W22 ).LT.SAFMIN ) THEN
1078:                   SCALE = ZERO
1079:                   U2 = ONE
1080:                   U1 = -W12 / W11
1081:                   GO TO 250
1082:                END IF
1083:                IF( ABS( W22 ).LT.ABS( U2 ) )
1084:      $            SCALE = ABS( W22 / U2 )
1085:                IF( ABS( W11 ).LT.ABS( U1 ) )
1086:      $            SCALE = MIN( SCALE, ABS( W11 / U1 ) )
1087: *
1088: *              Solve
1089: *
1090:                U2 = ( SCALE*U2 ) / W22
1091:                U1 = ( SCALE*U1-W12*U2 ) / W11
1092: *
1093:   250          CONTINUE
1094:                IF( ILPIVT ) THEN
1095:                   TEMP = U2
1096:                   U2 = U1
1097:                   U1 = TEMP
1098:                END IF
1099: *
1100: *              Compute Householder Vector
1101: *
1102:                T1 = SQRT( SCALE**2+U1**2+U2**2 )
1103:                TAU = ONE + SCALE / T1
1104:                VS = -ONE / ( SCALE+T1 )
1105:                V( 1 ) = ONE
1106:                V( 2 ) = VS*U1
1107:                V( 3 ) = VS*U2
1108: *
1109: *              Apply transformations from the right.
1110: *
1111:                DO 260 JR = IFRSTM, MIN( J+3, ILAST )
1112:                   TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )*
1113:      $                   H( JR, J+2 ) )
1114:                   H( JR, J ) = H( JR, J ) - TEMP
1115:                   H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 )
1116:                   H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 )
1117:   260          CONTINUE
1118:                DO 270 JR = IFRSTM, J + 2
1119:                   TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )*
1120:      $                   T( JR, J+2 ) )
1121:                   T( JR, J ) = T( JR, J ) - TEMP
1122:                   T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 )
1123:                   T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 )
1124:   270          CONTINUE
1125:                IF( ILZ ) THEN
1126:                   DO 280 JR = 1, N
1127:                      TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )*
1128:      $                      Z( JR, J+2 ) )
1129:                      Z( JR, J ) = Z( JR, J ) - TEMP
1130:                      Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 )
1131:                      Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 )
1132:   280             CONTINUE
1133:                END IF
1134:                T( J+1, J ) = ZERO
1135:                T( J+2, J ) = ZERO
1136:   290       CONTINUE
1137: *
1138: *           Last elements: Use Givens rotations
1139: *
1140: *           Rotations from the left
1141: *
1142:             J = ILAST - 1
1143:             TEMP = H( J, J-1 )
1144:             CALL SLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
1145:             H( J+1, J-1 ) = ZERO
1146: *
1147:             DO 300 JC = J, ILASTM
1148:                TEMP = C*H( J, JC ) + S*H( J+1, JC )
1149:                H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
1150:                H( J, JC ) = TEMP
1151:                TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
1152:                T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
1153:                T( J, JC ) = TEMP2
1154:   300       CONTINUE
1155:             IF( ILQ ) THEN
1156:                DO 310 JR = 1, N
1157:                   TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
1158:                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
1159:                   Q( JR, J ) = TEMP
1160:   310          CONTINUE
1161:             END IF
1162: *
1163: *           Rotations from the right.
1164: *
1165:             TEMP = T( J+1, J+1 )
1166:             CALL SLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
1167:             T( J+1, J ) = ZERO
1168: *
1169:             DO 320 JR = IFRSTM, ILAST
1170:                TEMP = C*H( JR, J+1 ) + S*H( JR, J )
1171:                H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
1172:                H( JR, J+1 ) = TEMP
1173:   320       CONTINUE
1174:             DO 330 JR = IFRSTM, ILAST - 1
1175:                TEMP = C*T( JR, J+1 ) + S*T( JR, J )
1176:                T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
1177:                T( JR, J+1 ) = TEMP
1178:   330       CONTINUE
1179:             IF( ILZ ) THEN
1180:                DO 340 JR = 1, N
1181:                   TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
1182:                   Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
1183:                   Z( JR, J+1 ) = TEMP
1184:   340          CONTINUE
1185:             END IF
1186: *
1187: *           End of Double-Shift code
1188: *
1189:          END IF
1190: *
1191:          GO TO 350
1192: *
1193: *        End of iteration loop
1194: *
1195:   350    CONTINUE
1196:   360 CONTINUE
1197: *
1198: *     Drop-through = non-convergence
1199: *
1200:       INFO = ILAST
1201:       GO TO 420
1202: *
1203: *     Successful completion of all QZ steps
1204: *
1205:   380 CONTINUE
1206: *
1207: *     Set Eigenvalues 1:ILO-1
1208: *
1209:       DO 410 J = 1, ILO - 1
1210:          IF( T( J, J ).LT.ZERO ) THEN
1211:             IF( ILSCHR ) THEN
1212:                DO 390 JR = 1, J
1213:                   H( JR, J ) = -H( JR, J )
1214:                   T( JR, J ) = -T( JR, J )
1215:   390          CONTINUE
1216:             ELSE
1217:                H( J, J ) = -H( J, J )
1218:                T( J, J ) = -T( J, J )
1219:             END IF
1220:             IF( ILZ ) THEN
1221:                DO 400 JR = 1, N
1222:                   Z( JR, J ) = -Z( JR, J )
1223:   400          CONTINUE
1224:             END IF
1225:          END IF
1226:          ALPHAR( J ) = H( J, J )
1227:          ALPHAI( J ) = ZERO
1228:          BETA( J ) = T( J, J )
1229:   410 CONTINUE
1230: *
1231: *     Normal Termination
1232: *
1233:       INFO = 0
1234: *
1235: *     Exit (other than argument error) -- return optimal workspace size
1236: *
1237:   420 CONTINUE
1238:       WORK( 1 ) = REAL( N )
1239:       RETURN
1240: *
1241: *     End of SHGEQZ
1242: *
1243:       END
1244: