001:       SUBROUTINE ZHGEQZ( 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:       DOUBLE PRECISION   RWORK( * )
015:       COMPLEX*16         ALPHA( * ), BETA( * ), H( LDH, * ),
016:      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
017:      $                   Z( LDZ, * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  ZHGEQZ 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 ZGGHRD.
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 ZGGHRD 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*16 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*16 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*16 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*16 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*16 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*16 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*16 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) DOUBLE PRECISION 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*16         CZERO, CONE
190:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
191:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
192:       DOUBLE PRECISION   ZERO, ONE
193:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
194:       DOUBLE PRECISION   HALF
195:       PARAMETER          ( HALF = 0.5D+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:       DOUBLE PRECISION   ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
203:      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
204:       COMPLEX*16         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:       DOUBLE PRECISION   DLAMCH, ZLANHS
211:       EXTERNAL           LSAME, DLAMCH, ZLANHS
212: *     ..
213: *     .. External Subroutines ..
214:       EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL
215: *     ..
216: *     .. Intrinsic Functions ..
217:       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN,
218:      $                   SQRT
219: *     ..
220: *     .. Statement Functions ..
221:       DOUBLE PRECISION   ABS1
222: *     ..
223: *     .. Statement Function definitions ..
224:       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
225: *     ..
226: *     .. Executable Statements ..
227: *
228: *     Decode JOB, COMPQ, COMPZ
229: *
230:       IF( LSAME( JOB, 'E' ) ) THEN
231:          ILSCHR = .FALSE.
232:          ISCHUR = 1
233:       ELSE IF( LSAME( JOB, 'S' ) ) THEN
234:          ILSCHR = .TRUE.
235:          ISCHUR = 2
236:       ELSE
237:          ISCHUR = 0
238:       END IF
239: *
240:       IF( LSAME( COMPQ, 'N' ) ) THEN
241:          ILQ = .FALSE.
242:          ICOMPQ = 1
243:       ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
244:          ILQ = .TRUE.
245:          ICOMPQ = 2
246:       ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
247:          ILQ = .TRUE.
248:          ICOMPQ = 3
249:       ELSE
250:          ICOMPQ = 0
251:       END IF
252: *
253:       IF( LSAME( COMPZ, 'N' ) ) THEN
254:          ILZ = .FALSE.
255:          ICOMPZ = 1
256:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
257:          ILZ = .TRUE.
258:          ICOMPZ = 2
259:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
260:          ILZ = .TRUE.
261:          ICOMPZ = 3
262:       ELSE
263:          ICOMPZ = 0
264:       END IF
265: *
266: *     Check Argument Values
267: *
268:       INFO = 0
269:       WORK( 1 ) = MAX( 1, N )
270:       LQUERY = ( LWORK.EQ.-1 )
271:       IF( ISCHUR.EQ.0 ) THEN
272:          INFO = -1
273:       ELSE IF( ICOMPQ.EQ.0 ) THEN
274:          INFO = -2
275:       ELSE IF( ICOMPZ.EQ.0 ) THEN
276:          INFO = -3
277:       ELSE IF( N.LT.0 ) THEN
278:          INFO = -4
279:       ELSE IF( ILO.LT.1 ) THEN
280:          INFO = -5
281:       ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
282:          INFO = -6
283:       ELSE IF( LDH.LT.N ) THEN
284:          INFO = -8
285:       ELSE IF( LDT.LT.N ) THEN
286:          INFO = -10
287:       ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
288:          INFO = -14
289:       ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
290:          INFO = -16
291:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
292:          INFO = -18
293:       END IF
294:       IF( INFO.NE.0 ) THEN
295:          CALL XERBLA( 'ZHGEQZ', -INFO )
296:          RETURN
297:       ELSE IF( LQUERY ) THEN
298:          RETURN
299:       END IF
300: *
301: *     Quick return if possible
302: *
303: *     WORK( 1 ) = CMPLX( 1 )
304:       IF( N.LE.0 ) THEN
305:          WORK( 1 ) = DCMPLX( 1 )
306:          RETURN
307:       END IF
308: *
309: *     Initialize Q and Z
310: *
311:       IF( ICOMPQ.EQ.3 )
312:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
313:       IF( ICOMPZ.EQ.3 )
314:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
315: *
316: *     Machine Constants
317: *
318:       IN = IHI + 1 - ILO
319:       SAFMIN = DLAMCH( 'S' )
320:       ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
321:       ANORM = ZLANHS( 'F', IN, H( ILO, ILO ), LDH, RWORK )
322:       BNORM = ZLANHS( 'F', IN, T( ILO, ILO ), LDT, RWORK )
323:       ATOL = MAX( SAFMIN, ULP*ANORM )
324:       BTOL = MAX( SAFMIN, ULP*BNORM )
325:       ASCALE = ONE / MAX( SAFMIN, ANORM )
326:       BSCALE = ONE / MAX( SAFMIN, BNORM )
327: *
328: *
329: *     Set Eigenvalues IHI+1:N
330: *
331:       DO 10 J = IHI + 1, N
332:          ABSB = ABS( T( J, J ) )
333:          IF( ABSB.GT.SAFMIN ) THEN
334:             SIGNBC = DCONJG( T( J, J ) / ABSB )
335:             T( J, J ) = ABSB
336:             IF( ILSCHR ) THEN
337:                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
338:                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
339:             ELSE
340:                H( J, J ) = H( J, J )*SIGNBC
341:             END IF
342:             IF( ILZ )
343:      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
344:          ELSE
345:             T( J, J ) = CZERO
346:          END IF
347:          ALPHA( J ) = H( J, J )
348:          BETA( J ) = T( J, J )
349:    10 CONTINUE
350: *
351: *     If IHI < ILO, skip QZ steps
352: *
353:       IF( IHI.LT.ILO )
354:      $   GO TO 190
355: *
356: *     MAIN QZ ITERATION LOOP
357: *
358: *     Initialize dynamic indices
359: *
360: *     Eigenvalues ILAST+1:N have been found.
361: *        Column operations modify rows IFRSTM:whatever
362: *        Row operations modify columns whatever:ILASTM
363: *
364: *     If only eigenvalues are being computed, then
365: *        IFRSTM is the row of the last splitting row above row ILAST;
366: *        this is always at least ILO.
367: *     IITER counts iterations since the last eigenvalue was found,
368: *        to tell when to use an extraordinary shift.
369: *     MAXIT is the maximum number of QZ sweeps allowed.
370: *
371:       ILAST = IHI
372:       IF( ILSCHR ) THEN
373:          IFRSTM = 1
374:          ILASTM = N
375:       ELSE
376:          IFRSTM = ILO
377:          ILASTM = IHI
378:       END IF
379:       IITER = 0
380:       ESHIFT = CZERO
381:       MAXIT = 30*( IHI-ILO+1 )
382: *
383:       DO 170 JITER = 1, MAXIT
384: *
385: *        Check for too many iterations.
386: *
387:          IF( JITER.GT.MAXIT )
388:      $      GO TO 180
389: *
390: *        Split the matrix if possible.
391: *
392: *        Two tests:
393: *           1: H(j,j-1)=0  or  j=ILO
394: *           2: T(j,j)=0
395: *
396: *        Special case: j=ILAST
397: *
398:          IF( ILAST.EQ.ILO ) THEN
399:             GO TO 60
400:          ELSE
401:             IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
402:                H( ILAST, ILAST-1 ) = CZERO
403:                GO TO 60
404:             END IF
405:          END IF
406: *
407:          IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
408:             T( ILAST, ILAST ) = CZERO
409:             GO TO 50
410:          END IF
411: *
412: *        General case: j<ILAST
413: *
414:          DO 40 J = ILAST - 1, ILO, -1
415: *
416: *           Test 1: for H(j,j-1)=0 or j=ILO
417: *
418:             IF( J.EQ.ILO ) THEN
419:                ILAZRO = .TRUE.
420:             ELSE
421:                IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN
422:                   H( J, J-1 ) = CZERO
423:                   ILAZRO = .TRUE.
424:                ELSE
425:                   ILAZRO = .FALSE.
426:                END IF
427:             END IF
428: *
429: *           Test 2: for T(j,j)=0
430: *
431:             IF( ABS( T( J, J ) ).LT.BTOL ) THEN
432:                T( J, J ) = CZERO
433: *
434: *              Test 1a: Check for 2 consecutive small subdiagonals in A
435: *
436:                ILAZR2 = .FALSE.
437:                IF( .NOT.ILAZRO ) THEN
438:                   IF( ABS1( H( J, J-1 ) )*( ASCALE*ABS1( H( J+1,
439:      $                J ) ) ).LE.ABS1( H( J, J ) )*( ASCALE*ATOL ) )
440:      $                ILAZR2 = .TRUE.
441:                END IF
442: *
443: *              If both tests pass (1 & 2), i.e., the leading diagonal
444: *              element of B in the block is zero, split a 1x1 block off
445: *              at the top. (I.e., at the J-th row/column) The leading
446: *              diagonal element of the remainder can also be zero, so
447: *              this may have to be done repeatedly.
448: *
449:                IF( ILAZRO .OR. ILAZR2 ) THEN
450:                   DO 20 JCH = J, ILAST - 1
451:                      CTEMP = H( JCH, JCH )
452:                      CALL ZLARTG( CTEMP, H( JCH+1, JCH ), C, S,
453:      $                            H( JCH, JCH ) )
454:                      H( JCH+1, JCH ) = CZERO
455:                      CALL ZROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
456:      $                          H( JCH+1, JCH+1 ), LDH, C, S )
457:                      CALL ZROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
458:      $                          T( JCH+1, JCH+1 ), LDT, C, S )
459:                      IF( ILQ )
460:      $                  CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
461:      $                             C, DCONJG( S ) )
462:                      IF( ILAZR2 )
463:      $                  H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
464:                      ILAZR2 = .FALSE.
465:                      IF( ABS1( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
466:                         IF( JCH+1.GE.ILAST ) THEN
467:                            GO TO 60
468:                         ELSE
469:                            IFIRST = JCH + 1
470:                            GO TO 70
471:                         END IF
472:                      END IF
473:                      T( JCH+1, JCH+1 ) = CZERO
474:    20             CONTINUE
475:                   GO TO 50
476:                ELSE
477: *
478: *                 Only test 2 passed -- chase the zero to T(ILAST,ILAST)
479: *                 Then process as in the case T(ILAST,ILAST)=0
480: *
481:                   DO 30 JCH = J, ILAST - 1
482:                      CTEMP = T( JCH, JCH+1 )
483:                      CALL ZLARTG( CTEMP, T( JCH+1, JCH+1 ), C, S,
484:      $                            T( JCH, JCH+1 ) )
485:                      T( JCH+1, JCH+1 ) = CZERO
486:                      IF( JCH.LT.ILASTM-1 )
487:      $                  CALL ZROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
488:      $                             T( JCH+1, JCH+2 ), LDT, C, S )
489:                      CALL ZROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
490:      $                          H( JCH+1, JCH-1 ), LDH, C, S )
491:                      IF( ILQ )
492:      $                  CALL ZROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
493:      $                             C, DCONJG( S ) )
494:                      CTEMP = H( JCH+1, JCH )
495:                      CALL ZLARTG( CTEMP, H( JCH+1, JCH-1 ), C, S,
496:      $                            H( JCH+1, JCH ) )
497:                      H( JCH+1, JCH-1 ) = CZERO
498:                      CALL ZROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
499:      $                          H( IFRSTM, JCH-1 ), 1, C, S )
500:                      CALL ZROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
501:      $                          T( IFRSTM, JCH-1 ), 1, C, S )
502:                      IF( ILZ )
503:      $                  CALL ZROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
504:      $                             C, S )
505:    30             CONTINUE
506:                   GO TO 50
507:                END IF
508:             ELSE IF( ILAZRO ) THEN
509: *
510: *              Only test 1 passed -- work on J:ILAST
511: *
512:                IFIRST = J
513:                GO TO 70
514:             END IF
515: *
516: *           Neither test passed -- try next J
517: *
518:    40    CONTINUE
519: *
520: *        (Drop-through is "impossible")
521: *
522:          INFO = 2*N + 1
523:          GO TO 210
524: *
525: *        T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
526: *        1x1 block.
527: *
528:    50    CONTINUE
529:          CTEMP = H( ILAST, ILAST )
530:          CALL ZLARTG( CTEMP, H( ILAST, ILAST-1 ), C, S,
531:      $                H( ILAST, ILAST ) )
532:          H( ILAST, ILAST-1 ) = CZERO
533:          CALL ZROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
534:      $              H( IFRSTM, ILAST-1 ), 1, C, S )
535:          CALL ZROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
536:      $              T( IFRSTM, ILAST-1 ), 1, C, S )
537:          IF( ILZ )
538:      $      CALL ZROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
539: *
540: *        H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
541: *
542:    60    CONTINUE
543:          ABSB = ABS( T( ILAST, ILAST ) )
544:          IF( ABSB.GT.SAFMIN ) THEN
545:             SIGNBC = DCONJG( T( ILAST, ILAST ) / ABSB )
546:             T( ILAST, ILAST ) = ABSB
547:             IF( ILSCHR ) THEN
548:                CALL ZSCAL( ILAST-IFRSTM, SIGNBC, T( IFRSTM, ILAST ), 1 )
549:                CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
550:      $                     1 )
551:             ELSE
552:                H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC
553:             END IF
554:             IF( ILZ )
555:      $         CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
556:          ELSE
557:             T( ILAST, ILAST ) = CZERO
558:          END IF
559:          ALPHA( ILAST ) = H( ILAST, ILAST )
560:          BETA( ILAST ) = T( ILAST, ILAST )
561: *
562: *        Go to next block -- exit if finished.
563: *
564:          ILAST = ILAST - 1
565:          IF( ILAST.LT.ILO )
566:      $      GO TO 190
567: *
568: *        Reset counters
569: *
570:          IITER = 0
571:          ESHIFT = CZERO
572:          IF( .NOT.ILSCHR ) THEN
573:             ILASTM = ILAST
574:             IF( IFRSTM.GT.ILAST )
575:      $         IFRSTM = ILO
576:          END IF
577:          GO TO 160
578: *
579: *        QZ step
580: *
581: *        This iteration only involves rows/columns IFIRST:ILAST.  We
582: *        assume IFIRST < ILAST, and that the diagonal of B is non-zero.
583: *
584:    70    CONTINUE
585:          IITER = IITER + 1
586:          IF( .NOT.ILSCHR ) THEN
587:             IFRSTM = IFIRST
588:          END IF
589: *
590: *        Compute the Shift.
591: *
592: *        At this point, IFIRST < ILAST, and the diagonal elements of
593: *        T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
594: *        magnitude)
595: *
596:          IF( ( IITER / 10 )*10.NE.IITER ) THEN
597: *
598: *           The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
599: *           the bottom-right 2x2 block of A inv(B) which is nearest to
600: *           the bottom-right element.
601: *
602: *           We factor B as U*D, where U has unit diagonals, and
603: *           compute (A*inv(D))*inv(U).
604: *
605:             U12 = ( BSCALE*T( ILAST-1, ILAST ) ) /
606:      $            ( BSCALE*T( ILAST, ILAST ) )
607:             AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
608:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
609:             AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
610:      $             ( BSCALE*T( ILAST-1, ILAST-1 ) )
611:             AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
612:      $             ( BSCALE*T( ILAST, ILAST ) )
613:             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
614:      $             ( BSCALE*T( ILAST, ILAST ) )
615:             ABI22 = AD22 - U12*AD21
616: *
617:             T1 = HALF*( AD11+ABI22 )
618:             RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )
619:             TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) +
620:      $             DIMAG( T1-ABI22 )*DIMAG( RTDISC )
621:             IF( TEMP.LE.ZERO ) THEN
622:                SHIFT = T1 + RTDISC
623:             ELSE
624:                SHIFT = T1 - RTDISC
625:             END IF
626:          ELSE
627: *
628: *           Exceptional shift.  Chosen for no particularly good reason.
629: *
630:             ESHIFT = ESHIFT + DCONJG( ( ASCALE*H( ILAST-1, ILAST ) ) /
631:      $               ( BSCALE*T( ILAST-1, ILAST-1 ) ) )
632:             SHIFT = ESHIFT
633:          END IF
634: *
635: *        Now check for two consecutive small subdiagonals.
636: *
637:          DO 80 J = ILAST - 1, IFIRST + 1, -1
638:             ISTART = J
639:             CTEMP = ASCALE*H( J, J ) - SHIFT*( BSCALE*T( J, J ) )
640:             TEMP = ABS1( CTEMP )
641:             TEMP2 = ASCALE*ABS1( H( J+1, J ) )
642:             TEMPR = MAX( TEMP, TEMP2 )
643:             IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
644:                TEMP = TEMP / TEMPR
645:                TEMP2 = TEMP2 / TEMPR
646:             END IF
647:             IF( ABS1( H( J, J-1 ) )*TEMP2.LE.TEMP*ATOL )
648:      $         GO TO 90
649:    80    CONTINUE
650: *
651:          ISTART = IFIRST
652:          CTEMP = ASCALE*H( IFIRST, IFIRST ) -
653:      $           SHIFT*( BSCALE*T( IFIRST, IFIRST ) )
654:    90    CONTINUE
655: *
656: *        Do an implicit-shift QZ sweep.
657: *
658: *        Initial Q
659: *
660:          CTEMP2 = ASCALE*H( ISTART+1, ISTART )
661:          CALL ZLARTG( CTEMP, CTEMP2, C, S, CTEMP3 )
662: *
663: *        Sweep
664: *
665:          DO 150 J = ISTART, ILAST - 1
666:             IF( J.GT.ISTART ) THEN
667:                CTEMP = H( J, J-1 )
668:                CALL ZLARTG( CTEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
669:                H( J+1, J-1 ) = CZERO
670:             END IF
671: *
672:             DO 100 JC = J, ILASTM
673:                CTEMP = C*H( J, JC ) + S*H( J+1, JC )
674:                H( J+1, JC ) = -DCONJG( S )*H( J, JC ) + C*H( J+1, JC )
675:                H( J, JC ) = CTEMP
676:                CTEMP2 = C*T( J, JC ) + S*T( J+1, JC )
677:                T( J+1, JC ) = -DCONJG( S )*T( J, JC ) + C*T( J+1, JC )
678:                T( J, JC ) = CTEMP2
679:   100       CONTINUE
680:             IF( ILQ ) THEN
681:                DO 110 JR = 1, N
682:                   CTEMP = C*Q( JR, J ) + DCONJG( S )*Q( JR, J+1 )
683:                   Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
684:                   Q( JR, J ) = CTEMP
685:   110          CONTINUE
686:             END IF
687: *
688:             CTEMP = T( J+1, J+1 )
689:             CALL ZLARTG( CTEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
690:             T( J+1, J ) = CZERO
691: *
692:             DO 120 JR = IFRSTM, MIN( J+2, ILAST )
693:                CTEMP = C*H( JR, J+1 ) + S*H( JR, J )
694:                H( JR, J ) = -DCONJG( S )*H( JR, J+1 ) + C*H( JR, J )
695:                H( JR, J+1 ) = CTEMP
696:   120       CONTINUE
697:             DO 130 JR = IFRSTM, J
698:                CTEMP = C*T( JR, J+1 ) + S*T( JR, J )
699:                T( JR, J ) = -DCONJG( S )*T( JR, J+1 ) + C*T( JR, J )
700:                T( JR, J+1 ) = CTEMP
701:   130       CONTINUE
702:             IF( ILZ ) THEN
703:                DO 140 JR = 1, N
704:                   CTEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
705:                   Z( JR, J ) = -DCONJG( S )*Z( JR, J+1 ) + C*Z( JR, J )
706:                   Z( JR, J+1 ) = CTEMP
707:   140          CONTINUE
708:             END IF
709:   150    CONTINUE
710: *
711:   160    CONTINUE
712: *
713:   170 CONTINUE
714: *
715: *     Drop-through = non-convergence
716: *
717:   180 CONTINUE
718:       INFO = ILAST
719:       GO TO 210
720: *
721: *     Successful completion of all QZ steps
722: *
723:   190 CONTINUE
724: *
725: *     Set Eigenvalues 1:ILO-1
726: *
727:       DO 200 J = 1, ILO - 1
728:          ABSB = ABS( T( J, J ) )
729:          IF( ABSB.GT.SAFMIN ) THEN
730:             SIGNBC = DCONJG( T( J, J ) / ABSB )
731:             T( J, J ) = ABSB
732:             IF( ILSCHR ) THEN
733:                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
734:                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
735:             ELSE
736:                H( J, J ) = H( J, J )*SIGNBC
737:             END IF
738:             IF( ILZ )
739:      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
740:          ELSE
741:             T( J, J ) = CZERO
742:          END IF
743:          ALPHA( J ) = H( J, J )
744:          BETA( J ) = T( J, J )
745:   200 CONTINUE
746: *
747: *     Normal Termination
748: *
749:       INFO = 0
750: *
751: *     Exit (other than argument error) -- return optimal workspace size
752: *
753:   210 CONTINUE
754:       WORK( 1 ) = DCMPLX( N )
755:       RETURN
756: *
757: *     End of ZHGEQZ
758: *
759:       END
760: