001:       SUBROUTINE CHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
002:      $                   ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
003:      $                   RWORK, INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          COMPQ, COMPZ, JOB
011:       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               RWORK( * )
015:       COMPLEX            ALPHA( * ), BETA( * ), H( LDH, * ),
016:      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
017:      $                   Z( LDZ, * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  CHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
024: *  where H is an upper Hessenberg matrix and T is upper triangular,
025: *  using the single-shift QZ method.
026: *  Matrix pairs of this type are produced by the reduction to
027: *  generalized upper Hessenberg form of a complex matrix pair (A,B):
028: *  
029: *     A = Q1*H*Z1**H,  B = Q1*T*Z1**H,
030: *  
031: *  as computed by CGGHRD.
032: *  
033: *  If JOB='S', then the Hessenberg-triangular pair (H,T) is
034: *  also reduced to generalized Schur form,
035: *  
036: *     H = Q*S*Z**H,  T = Q*P*Z**H,
037: *  
038: *  where Q and Z are unitary matrices and S and P are upper triangular.
039: *  
040: *  Optionally, the unitary matrix Q from the generalized Schur
041: *  factorization may be postmultiplied into an input matrix Q1, and the
042: *  unitary matrix Z may be postmultiplied into an input matrix Z1.
043: *  If Q1 and Z1 are the unitary matrices from CGGHRD that reduced
044: *  the matrix pair (A,B) to generalized Hessenberg form, then the output
045: *  matrices Q1*Q and Z1*Z are the unitary factors from the generalized
046: *  Schur factorization of (A,B):
047: *  
048: *     A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.
049: *  
050: *  To avoid overflow, eigenvalues of the matrix pair (H,T)
051: *  (equivalently, of (A,B)) are computed as a pair of complex values
052: *  (alpha,beta).  If beta is nonzero, lambda = alpha / beta is an
053: *  eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP)
054: *     A*x = lambda*B*x
055: *  and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
056: *  alternate form of the GNEP
057: *     mu*A*y = B*y.
058: *  The values of alpha and beta for the i-th eigenvalue can be read
059: *  directly from the generalized Schur form:  alpha = S(i,i),
060: *  beta = P(i,i).
061: *
062: *  Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
063: *       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
064: *       pp. 241--256.
065: *
066: *  Arguments
067: *  =========
068: *
069: *  JOB     (input) CHARACTER*1
070: *          = 'E': Compute eigenvalues only;
071: *          = 'S': Computer eigenvalues and the Schur form.
072: *
073: *  COMPQ   (input) CHARACTER*1
074: *          = 'N': Left Schur vectors (Q) are not computed;
075: *          = 'I': Q is initialized to the unit matrix and the matrix Q
076: *                 of left Schur vectors of (H,T) is returned;
077: *          = 'V': Q must contain a unitary matrix Q1 on entry and
078: *                 the product Q1*Q is returned.
079: *
080: *  COMPZ   (input) CHARACTER*1
081: *          = 'N': Right Schur vectors (Z) are not computed;
082: *          = 'I': Q is initialized to the unit matrix and the matrix Z
083: *                 of right Schur vectors of (H,T) is returned;
084: *          = 'V': Z must contain a unitary matrix Z1 on entry and
085: *                 the product Z1*Z is returned.
086: *
087: *  N       (input) INTEGER
088: *          The order of the matrices H, T, Q, and Z.  N >= 0.
089: *
090: *  ILO     (input) INTEGER
091: *  IHI     (input) INTEGER
092: *          ILO and IHI mark the rows and columns of H which are in
093: *          Hessenberg form.  It is assumed that A is already upper
094: *          triangular in rows and columns 1:ILO-1 and IHI+1:N.
095: *          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
096: *
097: *  H       (input/output) COMPLEX array, dimension (LDH, N)
098: *          On entry, the N-by-N upper Hessenberg matrix H.
099: *          On exit, if JOB = 'S', H contains the upper triangular
100: *          matrix S from the generalized Schur factorization.
101: *          If JOB = 'E', the diagonal of H matches that of S, but
102: *          the rest of H is unspecified.
103: *
104: *  LDH     (input) INTEGER
105: *          The leading dimension of the array H.  LDH >= max( 1, N ).
106: *
107: *  T       (input/output) COMPLEX array, dimension (LDT, N)
108: *          On entry, the N-by-N upper triangular matrix T.
109: *          On exit, if JOB = 'S', T contains the upper triangular
110: *          matrix P from the generalized Schur factorization.
111: *          If JOB = 'E', the diagonal of T matches that of P, but
112: *          the rest of T is unspecified.
113: *
114: *  LDT     (input) INTEGER
115: *          The leading dimension of the array T.  LDT >= max( 1, N ).
116: *
117: *  ALPHA   (output) COMPLEX array, dimension (N)
118: *          The complex scalars alpha that define the eigenvalues of
119: *          GNEP.  ALPHA(i) = S(i,i) in the generalized Schur
120: *          factorization.
121: *
122: *  BETA    (output) COMPLEX array, dimension (N)
123: *          The real non-negative scalars beta that define the
124: *          eigenvalues of GNEP.  BETA(i) = P(i,i) in the generalized
125: *          Schur factorization.
126: *
127: *          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
128: *          represent the j-th eigenvalue of the matrix pair (A,B), in
129: *          one of the forms lambda = alpha/beta or mu = beta/alpha.
130: *          Since either lambda or mu may overflow, they should not,
131: *          in general, be computed.
132: *
133: *  Q       (input/output) COMPLEX array, dimension (LDQ, N)
134: *          On entry, if COMPZ = 'V', the unitary matrix Q1 used in the
135: *          reduction of (A,B) to generalized Hessenberg form.
136: *          On exit, if COMPZ = 'I', the unitary matrix of left Schur
137: *          vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
138: *          left Schur vectors of (A,B).
139: *          Not referenced if COMPZ = 'N'.
140: *
141: *  LDQ     (input) INTEGER
142: *          The leading dimension of the array Q.  LDQ >= 1.
143: *          If COMPQ='V' or 'I', then LDQ >= N.
144: *
145: *  Z       (input/output) COMPLEX array, dimension (LDZ, N)
146: *          On entry, if COMPZ = 'V', the unitary matrix Z1 used in the
147: *          reduction of (A,B) to generalized Hessenberg form.
148: *          On exit, if COMPZ = 'I', the unitary matrix of right Schur
149: *          vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
150: *          right Schur vectors of (A,B).
151: *          Not referenced if COMPZ = 'N'.
152: *
153: *  LDZ     (input) INTEGER
154: *          The leading dimension of the array Z.  LDZ >= 1.
155: *          If COMPZ='V' or 'I', then LDZ >= N.
156: *
157: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
158: *          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
159: *
160: *  LWORK   (input) INTEGER
161: *          The dimension of the array WORK.  LWORK >= max(1,N).
162: *
163: *          If LWORK = -1, then a workspace query is assumed; the routine
164: *          only calculates the optimal size of the WORK array, returns
165: *          this value as the first entry of the WORK array, and no error
166: *          message related to LWORK is issued by XERBLA.
167: *
168: *  RWORK   (workspace) REAL array, dimension (N)
169: *
170: *  INFO    (output) INTEGER
171: *          = 0: successful exit
172: *          < 0: if INFO = -i, the i-th argument had an illegal value
173: *          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
174: *                     in Schur form, but ALPHA(i) and BETA(i),
175: *                     i=INFO+1,...,N should be correct.
176: *          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
177: *                     in Schur form, but ALPHA(i) and BETA(i),
178: *                     i=INFO-N+1,...,N should be correct.
179: *
180: *  Further Details
181: *  ===============
182: *
183: *  We assume that complex ABS works as long as its value is less than
184: *  overflow.
185: *
186: *  =====================================================================
187: *
188: *     .. Parameters ..
189:       COMPLEX            CZERO, CONE
190:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
191:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
192:       REAL               ZERO, ONE
193:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
194:       REAL               HALF
195:       PARAMETER          ( HALF = 0.5E+0 )
196: *     ..
197: *     .. Local Scalars ..
198:       LOGICAL            ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
199:       INTEGER            ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
200:      $                   ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
201:      $                   JR, MAXIT
202:       REAL               ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
203:      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
204:       COMPLEX            ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
205:      $                   CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,
206:      $                   U12, X
207: *     ..
208: *     .. External Functions ..
209:       LOGICAL            LSAME
210:       REAL               CLANHS, SLAMCH
211:       EXTERNAL           LSAME, CLANHS, SLAMCH
212: *     ..
213: *     .. External Subroutines ..
214:       EXTERNAL           CLARTG, CLASET, CROT, CSCAL, XERBLA
215: *     ..
216: *     .. Intrinsic Functions ..
217:       INTRINSIC          ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL, SQRT
218: *     ..
219: *     .. Statement Functions ..
220:       REAL               ABS1
221: *     ..
222: *     .. Statement Function definitions ..
223:       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
224: *     ..
225: *     .. Executable Statements ..
226: *
227: *     Decode JOB, COMPQ, COMPZ
228: *
229:       IF( LSAME( JOB, 'E' ) ) THEN
230:          ILSCHR = .FALSE.
231:          ISCHUR = 1
232:       ELSE IF( LSAME( JOB, 'S' ) ) THEN
233:          ILSCHR = .TRUE.
234:          ISCHUR = 2
235:       ELSE
236:          ISCHUR = 0
237:       END IF
238: *
239:       IF( LSAME( COMPQ, 'N' ) ) THEN
240:          ILQ = .FALSE.
241:          ICOMPQ = 1
242:       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
243:          ILQ = .TRUE.
244:          ICOMPQ = 2
245:       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
246:          ILQ = .TRUE.
247:          ICOMPQ = 3
248:       ELSE
249:          ICOMPQ = 0
250:       END IF
251: *
252:       IF( LSAME( COMPZ, 'N' ) ) THEN
253:          ILZ = .FALSE.
254:          ICOMPZ = 1
255:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
256:          ILZ = .TRUE.
257:          ICOMPZ = 2
258:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
259:          ILZ = .TRUE.
260:          ICOMPZ = 3
261:       ELSE
262:          ICOMPZ = 0
263:       END IF
264: *
265: *     Check Argument Values
266: *
267:       INFO = 0
268:       WORK( 1 ) = MAX( 1, N )
269:       LQUERY = ( LWORK.EQ.-1 )
270:       IF( ISCHUR.EQ.0 ) THEN
271:          INFO = -1
272:       ELSE IF( ICOMPQ.EQ.0 ) THEN
273:          INFO = -2
274:       ELSE IF( ICOMPZ.EQ.0 ) THEN
275:          INFO = -3
276:       ELSE IF( N.LT.0 ) THEN
277:          INFO = -4
278:       ELSE IF( ILO.LT.1 ) THEN
279:          INFO = -5
280:       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
281:          INFO = -6
282:       ELSE IF( LDH.LT.N ) THEN
283:          INFO = -8
284:       ELSE IF( LDT.LT.N ) THEN
285:          INFO = -10
286:       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
287:          INFO = -14
288:       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
289:          INFO = -16
290:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
291:          INFO = -18
292:       END IF
293:       IF( INFO.NE.0 ) THEN
294:          CALL XERBLA( 'CHGEQZ', -INFO )
295:          RETURN
296:       ELSE IF( LQUERY ) THEN
297:          RETURN
298:       END IF
299: *
300: *     Quick return if possible
301: *
302: c     WORK( 1 ) = CMPLX( 1 )
303:       IF( N.LE.0 ) THEN
304:          WORK( 1 ) = CMPLX( 1 )
305:          RETURN
306:       END IF
307: *
308: *     Initialize Q and Z
309: *
310:       IF( ICOMPQ.EQ.3 )
311:      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
312:       IF( ICOMPZ.EQ.3 )
313:      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
314: *
315: *     Machine Constants
316: *
317:       IN = IHI + 1 - ILO
318:       SAFMIN = SLAMCH( 'S' )
319:       ULP = SLAMCH( 'E' )*SLAMCH( 'B' )
320:       ANORM = CLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK )
321:       BNORM = CLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK )
322:       ATOL = MAX( SAFMIN, ULP*ANORM )
323:       BTOL = MAX( SAFMIN, ULP*BNORM )
324:       ASCALE = ONE / MAX( SAFMIN, ANORM )
325:       BSCALE = ONE / MAX( SAFMIN, BNORM )
326: *
327: *
328: *     Set Eigenvalues IHI+1:N
329: *
330:       DO 10 J = IHI + 1, N
331:          ABSB = ABS( T( J, J ) )
332:          IF( ABSB.GT.SAFMIN ) THEN
333:             SIGNBC = CONJG( T( J, J ) / ABSB )
334:             T( J, J ) = ABSB
335:             IF( ILSCHR ) THEN
336:                CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
337:                CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
338:             ELSE
339:                H( J, J ) = H( J, J )*SIGNBC
340:             END IF
341:             IF( ILZ )
342:      $         CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
343:          ELSE
344:             T( J, J ) = CZERO
345:          END IF
346:          ALPHA( J ) = H( J, J )
347:          BETA( J ) = T( J, J )
348:    10 CONTINUE
349: *
350: *     If IHI < ILO, skip QZ steps
351: *
352:       IF( IHI.LT.ILO )
353:      $   GO TO 190
354: *
355: *     MAIN QZ ITERATION LOOP
356: *
357: *     Initialize dynamic indices
358: *
359: *     Eigenvalues ILAST+1:N have been found.
360: *        Column operations modify rows IFRSTM:whatever
361: *        Row operations modify columns whatever:ILASTM
362: *
363: *     If only eigenvalues are being computed, then
364: *        IFRSTM is the row of the last splitting row above row ILAST;
365: *        this is always at least ILO.
366: *     IITER counts iterations since the last eigenvalue was found,
367: *        to tell when to use an extraordinary shift.
368: *     MAXIT is the maximum number of QZ sweeps allowed.
369: *
370:       ILAST = IHI
371:       IF( ILSCHR ) THEN
372:          IFRSTM = 1
373:          ILASTM = N
374:       ELSE
375:          IFRSTM = ILO
376:          ILASTM = IHI
377:       END IF
378:       IITER = 0
379:       ESHIFT = CZERO
380:       MAXIT = 30*( IHI-ILO+1 )
381: *
382:       DO 170 JITER = 1, MAXIT
383: *
384: *        Check for too many iterations.
385: *
386:          IF( JITER.GT.MAXIT )
387:      $      GO TO 180
388: *
389: *        Split the matrix if possible.
390: *
391: *        Two tests:
392: *           1: H(j,j-1)=0  or  j=ILO
393: *           2: T(j,j)=0
394: *
395: *        Special case: j=ILAST
396: *
397:          IF( ILAST.EQ.ILO ) THEN
398:             GO TO 60
399:          ELSE
400:             IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
401:                H( ILAST, ILAST-1 ) = CZERO
402:                GO TO 60
403:             END IF
404:          END IF
405: *
406:          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
407:             T( ILAST, ILAST ) = CZERO
408:             GO TO 50
409:          END IF
410: *
411: *        General case: j<ILAST
412: *
413:          DO 40 J = ILAST - 1, ILO, -1
414: *
415: *           Test 1: for H(j,j-1)=0 or j=ILO
416: *
417:             IF( J.EQ.ILO ) THEN
418:                ILAZRO = .TRUE.
419:             ELSE
420:                IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN
421:                   H( J, J-1 ) = CZERO
422:                   ILAZRO = .TRUE.
423:                ELSE
424:                   ILAZRO = .FALSE.
425:                END IF
426:             END IF
427: *
428: *           Test 2: for T(j,j)=0
429: *
430:             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
431:                T( J, J ) = CZERO
432: *
433: *              Test 1a: Check for 2 consecutive small subdiagonals in A
434: *
435:                ILAZR2 = .FALSE.
436:                IF( .NOT.ILAZRO ) THEN
437:                   IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1,
438:      $                J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) )
439:      $                ILAZR2 = .TRUE.
440:                END IF
441: *
442: *              If both tests pass (1 & 2), i.e., the leading diagonal
443: *              element of B in the block is zero, split a 1x1 block off
444: *              at the top. (I.e., at the J-th row/column) The leading
445: *              diagonal element of the remainder can also be zero, so
446: *              this may have to be done repeatedly.
447: *
448:                IF( ILAZRO .OR. ILAZR2 ) THEN
449:                   DO 20 JCH = J, ILAST - 1
450:                      CTEMP = H( JCH, JCH )
451:                      CALL CLARTG( CTEMP, H( JCH+1, JCH ), C, S,
452:      $                            H( JCH, JCH ) )
453:                      H( JCH+1, JCH ) = CZERO
454:                      CALL CROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
455:      $                          H( JCH+1, JCH+1 ), LDH, C, S )
456:                      CALL CROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
457:      $                          T( JCH+1, JCH+1 ), LDT, C, S )
458:                      IF( ILQ )
459:      $                  CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
460:      $                             C, CONJG( S ) )
461:                      IF( ILAZR2 )
462:      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
463:                      ILAZR2 = .FALSE.
464:                      IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
465:                         IF( JCH+1.GE.ILAST ) THEN
466:                            GO TO 60
467:                         ELSE
468:                            IFIRST = JCH + 1
469:                            GO TO 70
470:                         END IF
471:                      END IF
472:                      T( JCH+1, JCH+1 ) = CZERO
473:    20             CONTINUE
474:                   GO TO 50
475:                ELSE
476: *
477: *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
478: *                 Then process as in the case T(ILAST,ILAST)=0
479: *
480:                   DO 30 JCH = J, ILAST - 1
481:                      CTEMP = T( JCH, JCH+1 )
482:                      CALL CLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S,
483:      $                            T( JCH, JCH+1 ) )
484:                      T( JCH+1, JCH+1 ) = CZERO
485:                      IF( JCH.LT.ILASTM-1 )
486:      $                  CALL CROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
487:      $                             T( JCH+1, JCH+2 ), LDT, C, S )
488:                      CALL CROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
489:      $                          H( JCH+1, JCH-1 ), LDH, C, S )
490:                      IF( ILQ )
491:      $                  CALL CROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
492:      $                             C, CONJG( S ) )
493:                      CTEMP = H( JCH+1, JCH )
494:                      CALL CLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S,
495:      $                            H( JCH+1, JCH ) )
496:                      H( JCH+1, JCH-1 ) = CZERO
497:                      CALL CROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
498:      $                          H( IFRSTM, JCH-1 ), 1, C, S )
499:                      CALL CROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
500:      $                          T( IFRSTM, JCH-1 ), 1, C, S )
501:                      IF( ILZ )
502:      $                  CALL CROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
503:      $                             C, S )
504:    30             CONTINUE
505:                   GO TO 50
506:                END IF
507:             ELSE IF( ILAZRO ) THEN
508: *
509: *              Only test 1 passed -- work on J:ILAST
510: *
511:                IFIRST = J
512:                GO TO 70
513:             END IF
514: *
515: *           Neither test passed -- try next J
516: *
517:    40    CONTINUE
518: *
519: *        (Drop-through is "impossible")
520: *
521:          INFO = 2*N + 1
522:          GO TO 210
523: *
524: *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
525: *        1x1 block.
526: *
527:    50    CONTINUE
528:          CTEMP = H( ILAST, ILAST )
529:          CALL CLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S,
530:      $                H( ILAST, ILAST ) )
531:          H( ILAST, ILAST-1 ) = CZERO
532:          CALL CROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
533:      $              H( IFRSTM, ILAST-1 ), 1, C, S )
534:          CALL CROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
535:      $              T( IFRSTM, ILAST-1 ), 1, C, S )
536:          IF( ILZ )
537:      $      CALL CROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
538: *
539: *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
540: *
541:    60    CONTINUE
542:          ABSB = ABS( T( ILAST, ILAST ) )
543:          IF( ABSB.GT.SAFMIN ) THEN
544:             SIGNBC = CONJG( T( ILAST, ILAST ) / ABSB )
545:             T( ILAST, ILAST ) = ABSB
546:             IF( ILSCHR ) THEN
547:                CALL CSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 )
548:                CALL CSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
549:      $                     1 )
550:             ELSE
551:                H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC
552:             END IF
553:             IF( ILZ )
554:      $         CALL CSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
555:          ELSE
556:             T( ILAST, ILAST ) = CZERO
557:          END IF
558:          ALPHA( ILAST ) = H( ILAST, ILAST )
559:          BETA( ILAST ) = T( ILAST, ILAST )
560: *
561: *        Go to next block -- exit if finished.
562: *
563:          ILAST = ILAST - 1
564:          IF( ILAST.LT.ILO )
565:      $      GO TO 190
566: *
567: *        Reset counters
568: *
569:          IITER = 0
570:          ESHIFT = CZERO
571:          IF( .NOT.ILSCHR ) THEN
572:             ILASTM = ILAST
573:             IF( IFRSTM.GT.ILAST )
574:      $         IFRSTM = ILO
575:          END IF
576:          GO TO 160
577: *
578: *        QZ step
579: *
580: *        This iteration only involves rows/columns IFIRST:ILAST.  We
581: *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
582: *
583:    70    CONTINUE
584:          IITER = IITER + 1
585:          IF( .NOT.ILSCHR ) THEN
586:             IFRSTM = IFIRST
587:          END IF
588: *
589: *        Compute the Shift.
590: *
591: *        At this point, IFIRST < ILAST, and the diagonal elements of
592: *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
593: *        magnitude)
594: *
595:          IF( ( IITER / 10 )*10.NE.IITER ) THEN
596: *
597: *           The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
598: *           the bottom-right 2x2 block of A inv(B) which is nearest to
599: *           the bottom-right element.
600: *
601: *           We factor B as U*D, where U has unit diagonals, and
602: *           compute (A*inv(D))*inv(U).
603: *
604:             U12 = ( BSCALE*T( ILAST-1, ILAST ) ) /
605:      $            ( BSCALE*T( ILAST, ILAST ) )
606:             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
607:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
608:             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
609:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
610:             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
611:      $             ( BSCALE*T( ILAST, ILAST ) )
612:             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
613:      $             ( BSCALE*T( ILAST, ILAST ) )
614:             ABI22 = AD22 - U12*AD21
615: *
616:             T1 = HALF*( AD11+ABI22 )
617:             RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
618:             TEMP = REAL( T1-ABI22 )*REAL( RTDISC ) +
619:      $             AIMAG( T1-ABI22 )*AIMAG( RTDISC )
620:             IF( TEMP.LE.ZERO ) THEN
621:                SHIFT = T1 + RTDISC
622:             ELSE
623:                SHIFT = T1 - RTDISC
624:             END IF
625:          ELSE
626: *
627: *           Exceptional shift.  Chosen for no particularly good reason.
628: *
629:             ESHIFT = ESHIFT + CONJG( ( ASCALE*H( ILAST-1, ILAST ) ) /
630:      $               ( BSCALE*T( ILAST-1, ILAST-1 ) ) )
631:             SHIFT = ESHIFT
632:          END IF
633: *
634: *        Now check for two consecutive small subdiagonals.
635: *
636:          DO 80 J = ILAST - 1, IFIRST + 1, -1
637:             ISTART = J
638:             CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) )
639:             TEMP = ABS1( CTEMP )
640:             TEMP2 = ASCALE*ABS1( H( J+1, J ) )
641:             TEMPR = MAX( TEMP, TEMP2 )
642:             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
643:                TEMP = TEMP / TEMPR
644:                TEMP2 = TEMP2 / TEMPR
645:             END IF
646:             IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL )
647:      $         GO TO 90
648:    80    CONTINUE
649: *
650:          ISTART = IFIRST
651:          CTEMP = ASCALE*H( IFIRST, IFIRST ) -
652:      $           SHIFT*( BSCALE*T( IFIRST, IFIRST ) )
653:    90    CONTINUE
654: *
655: *        Do an implicit-shift QZ sweep.
656: *
657: *        Initial Q
658: *
659:          CTEMP2 = ASCALE*H( ISTART+1, ISTART )
660:          CALL CLARTG( CTEMP, CTEMP2, C, S, CTEMP3 )
661: *
662: *        Sweep
663: *
664:          DO 150 J = ISTART, ILAST - 1
665:             IF( J.GT.ISTART ) THEN
666:                CTEMP = H( J, J-1 )
667:                CALL CLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
668:                H( J+1, J-1 ) = CZERO
669:             END IF
670: *
671:             DO 100 JC = J, ILASTM
672:                CTEMP = C*H( J, JC ) + S*H( J+1, JC )
673:                H( J+1, JC ) = -CONJG( S )*H( J, JC ) + C*H( J+1, JC )
674:                H( J, JC ) = CTEMP
675:                CTEMP2 = C*T( J, JC ) + S*T( J+1, JC )
676:                T( J+1, JC ) = -CONJG( S )*T( J, JC ) + C*T( J+1, JC )
677:                T( J, JC ) = CTEMP2
678:   100       CONTINUE
679:             IF( ILQ ) THEN
680:                DO 110 JR = 1, N
681:                   CTEMP = C*Q( JR, J ) + CONJG( S )*Q( JR, J+1 )
682:                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
683:                   Q( JR, J ) = CTEMP
684:   110          CONTINUE
685:             END IF
686: *
687:             CTEMP = T( J+1, J+1 )
688:             CALL CLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
689:             T( J+1, J ) = CZERO
690: *
691:             DO 120 JR = IFRSTM, MIN( J+2, ILAST )
692:                CTEMP = C*H( JR, J+1 ) + S*H( JR, J )
693:                H( JR, J ) = -CONJG( S )*H( JR, J+1 ) + C*H( JR, J )
694:                H( JR, J+1 ) = CTEMP
695:   120       CONTINUE
696:             DO 130 JR = IFRSTM, J
697:                CTEMP = C*T( JR, J+1 ) + S*T( JR, J )
698:                T( JR, J ) = -CONJG( S )*T( JR, J+1 ) + C*T( JR, J )
699:                T( JR, J+1 ) = CTEMP
700:   130       CONTINUE
701:             IF( ILZ ) THEN
702:                DO 140 JR = 1, N
703:                   CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
704:                   Z( JR, J ) = -CONJG( S )*Z( JR, J+1 ) + C*Z( JR, J )
705:                   Z( JR, J+1 ) = CTEMP
706:   140          CONTINUE
707:             END IF
708:   150    CONTINUE
709: *
710:   160    CONTINUE
711: *
712:   170 CONTINUE
713: *
714: *     Drop-through = non-convergence
715: *
716:   180 CONTINUE
717:       INFO = ILAST
718:       GO TO 210
719: *
720: *     Successful completion of all QZ steps
721: *
722:   190 CONTINUE
723: *
724: *     Set Eigenvalues 1:ILO-1
725: *
726:       DO 200 J = 1, ILO - 1
727:          ABSB = ABS( T( J, J ) )
728:          IF( ABSB.GT.SAFMIN ) THEN
729:             SIGNBC = CONJG( T( J, J ) / ABSB )
730:             T( J, J ) = ABSB
731:             IF( ILSCHR ) THEN
732:                CALL CSCAL( J-1, SIGNBC, T( 1, J ), 1 )
733:                CALL CSCAL( J, SIGNBC, H( 1, J ), 1 )
734:             ELSE
735:                H( J, J ) = H( J, J )*SIGNBC
736:             END IF
737:             IF( ILZ )
738:      $         CALL CSCAL( N, SIGNBC, Z( 1, J ), 1 )
739:          ELSE
740:             T( J, J ) = CZERO
741:          END IF
742:          ALPHA( J ) = H( J, J )
743:          BETA( J ) = T( J, J )
744:   200 CONTINUE
745: *
746: *     Normal Termination
747: *
748:       INFO = 0
749: *
750: *     Exit (other than argument error) -- return optimal workspace size
751: *
752:   210 CONTINUE
753:       WORK( 1 ) = CMPLX( N )
754:       RETURN
755: *
756: *     End of CHGEQZ
757: *
758:       END
759: