LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
zhgeqz.f
Go to the documentation of this file.
1 *> \brief \b ZHGEQZ
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZHGEQZ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhgeqz.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhgeqz.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhgeqz.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
22 * ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
23 * RWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER COMPQ, COMPZ, JOB
27 * INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
28 * ..
29 * .. Array Arguments ..
30 * DOUBLE PRECISION RWORK( * )
31 * COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ),
32 * $ Q( LDQ, * ), T( LDT, * ), WORK( * ),
33 * $ Z( LDZ, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZHGEQZ computes the eigenvalues of a complex matrix pair (H,T),
43 *> where H is an upper Hessenberg matrix and T is upper triangular,
44 *> using the single-shift QZ method.
45 *> Matrix pairs of this type are produced by the reduction to
46 *> generalized upper Hessenberg form of a complex matrix pair (A,B):
47 *>
48 *> A = Q1*H*Z1**H, B = Q1*T*Z1**H,
49 *>
50 *> as computed by ZGGHRD.
51 *>
52 *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
53 *> also reduced to generalized Schur form,
54 *>
55 *> H = Q*S*Z**H, T = Q*P*Z**H,
56 *>
57 *> where Q and Z are unitary matrices and S and P are upper triangular.
58 *>
59 *> Optionally, the unitary matrix Q from the generalized Schur
60 *> factorization may be postmultiplied into an input matrix Q1, and the
61 *> unitary matrix Z may be postmultiplied into an input matrix Z1.
62 *> If Q1 and Z1 are the unitary matrices from ZGGHRD that reduced
63 *> the matrix pair (A,B) to generalized Hessenberg form, then the output
64 *> matrices Q1*Q and Z1*Z are the unitary factors from the generalized
65 *> Schur factorization of (A,B):
66 *>
67 *> A = (Q1*Q)*S*(Z1*Z)**H, B = (Q1*Q)*P*(Z1*Z)**H.
68 *>
69 *> To avoid overflow, eigenvalues of the matrix pair (H,T)
70 *> (equivalently, of (A,B)) are computed as a pair of complex values
71 *> (alpha,beta). If beta is nonzero, lambda = alpha / beta is an
72 *> eigenvalue of the generalized nonsymmetric eigenvalue problem (GNEP)
73 *> A*x = lambda*B*x
74 *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
75 *> alternate form of the GNEP
76 *> mu*A*y = B*y.
77 *> The values of alpha and beta for the i-th eigenvalue can be read
78 *> directly from the generalized Schur form: alpha = S(i,i),
79 *> beta = P(i,i).
80 *>
81 *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
82 *> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
83 *> pp. 241--256.
84 *> \endverbatim
85 *
86 * Arguments:
87 * ==========
88 *
89 *> \param[in] JOB
90 *> \verbatim
91 *> JOB is CHARACTER*1
92 *> = 'E': Compute eigenvalues only;
93 *> = 'S': Computer eigenvalues and the Schur form.
94 *> \endverbatim
95 *>
96 *> \param[in] COMPQ
97 *> \verbatim
98 *> COMPQ is CHARACTER*1
99 *> = 'N': Left Schur vectors (Q) are not computed;
100 *> = 'I': Q is initialized to the unit matrix and the matrix Q
101 *> of left Schur vectors of (H,T) is returned;
102 *> = 'V': Q must contain a unitary matrix Q1 on entry and
103 *> the product Q1*Q is returned.
104 *> \endverbatim
105 *>
106 *> \param[in] COMPZ
107 *> \verbatim
108 *> COMPZ is CHARACTER*1
109 *> = 'N': Right Schur vectors (Z) are not computed;
110 *> = 'I': Q is initialized to the unit matrix and the matrix Z
111 *> of right Schur vectors of (H,T) is returned;
112 *> = 'V': Z must contain a unitary matrix Z1 on entry and
113 *> the product Z1*Z is returned.
114 *> \endverbatim
115 *>
116 *> \param[in] N
117 *> \verbatim
118 *> N is INTEGER
119 *> The order of the matrices H, T, Q, and Z. N >= 0.
120 *> \endverbatim
121 *>
122 *> \param[in] ILO
123 *> \verbatim
124 *> ILO is INTEGER
125 *> \endverbatim
126 *>
127 *> \param[in] IHI
128 *> \verbatim
129 *> IHI is INTEGER
130 *> ILO and IHI mark the rows and columns of H which are in
131 *> Hessenberg form. It is assumed that A is already upper
132 *> triangular in rows and columns 1:ILO-1 and IHI+1:N.
133 *> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
134 *> \endverbatim
135 *>
136 *> \param[in,out] H
137 *> \verbatim
138 *> H is COMPLEX*16 array, dimension (LDH, N)
139 *> On entry, the N-by-N upper Hessenberg matrix H.
140 *> On exit, if JOB = 'S', H contains the upper triangular
141 *> matrix S from the generalized Schur factorization.
142 *> If JOB = 'E', the diagonal of H matches that of S, but
143 *> the rest of H is unspecified.
144 *> \endverbatim
145 *>
146 *> \param[in] LDH
147 *> \verbatim
148 *> LDH is INTEGER
149 *> The leading dimension of the array H. LDH >= max( 1, N ).
150 *> \endverbatim
151 *>
152 *> \param[in,out] T
153 *> \verbatim
154 *> T is COMPLEX*16 array, dimension (LDT, N)
155 *> On entry, the N-by-N upper triangular matrix T.
156 *> On exit, if JOB = 'S', T contains the upper triangular
157 *> matrix P from the generalized Schur factorization.
158 *> If JOB = 'E', the diagonal of T matches that of P, but
159 *> the rest of T is unspecified.
160 *> \endverbatim
161 *>
162 *> \param[in] LDT
163 *> \verbatim
164 *> LDT is INTEGER
165 *> The leading dimension of the array T. LDT >= max( 1, N ).
166 *> \endverbatim
167 *>
168 *> \param[out] ALPHA
169 *> \verbatim
170 *> ALPHA is COMPLEX*16 array, dimension (N)
171 *> The complex scalars alpha that define the eigenvalues of
172 *> GNEP. ALPHA(i) = S(i,i) in the generalized Schur
173 *> factorization.
174 *> \endverbatim
175 *>
176 *> \param[out] BETA
177 *> \verbatim
178 *> BETA is COMPLEX*16 array, dimension (N)
179 *> The real non-negative scalars beta that define the
180 *> eigenvalues of GNEP. BETA(i) = P(i,i) in the generalized
181 *> Schur factorization.
182 *>
183 *> Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
184 *> represent the j-th eigenvalue of the matrix pair (A,B), in
185 *> one of the forms lambda = alpha/beta or mu = beta/alpha.
186 *> Since either lambda or mu may overflow, they should not,
187 *> in general, be computed.
188 *> \endverbatim
189 *>
190 *> \param[in,out] Q
191 *> \verbatim
192 *> Q is COMPLEX*16 array, dimension (LDQ, N)
193 *> On entry, if COMPZ = 'V', the unitary matrix Q1 used in the
194 *> reduction of (A,B) to generalized Hessenberg form.
195 *> On exit, if COMPZ = 'I', the unitary matrix of left Schur
196 *> vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
197 *> left Schur vectors of (A,B).
198 *> Not referenced if COMPZ = 'N'.
199 *> \endverbatim
200 *>
201 *> \param[in] LDQ
202 *> \verbatim
203 *> LDQ is INTEGER
204 *> The leading dimension of the array Q. LDQ >= 1.
205 *> If COMPQ='V' or 'I', then LDQ >= N.
206 *> \endverbatim
207 *>
208 *> \param[in,out] Z
209 *> \verbatim
210 *> Z is COMPLEX*16 array, dimension (LDZ, N)
211 *> On entry, if COMPZ = 'V', the unitary matrix Z1 used in the
212 *> reduction of (A,B) to generalized Hessenberg form.
213 *> On exit, if COMPZ = 'I', the unitary matrix of right Schur
214 *> vectors of (H,T), and if COMPZ = 'V', the unitary matrix of
215 *> right Schur vectors of (A,B).
216 *> Not referenced if COMPZ = 'N'.
217 *> \endverbatim
218 *>
219 *> \param[in] LDZ
220 *> \verbatim
221 *> LDZ is INTEGER
222 *> The leading dimension of the array Z. LDZ >= 1.
223 *> If COMPZ='V' or 'I', then LDZ >= N.
224 *> \endverbatim
225 *>
226 *> \param[out] WORK
227 *> \verbatim
228 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
229 *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
230 *> \endverbatim
231 *>
232 *> \param[in] LWORK
233 *> \verbatim
234 *> LWORK is INTEGER
235 *> The dimension of the array WORK. LWORK >= max(1,N).
236 *>
237 *> If LWORK = -1, then a workspace query is assumed; the routine
238 *> only calculates the optimal size of the WORK array, returns
239 *> this value as the first entry of the WORK array, and no error
240 *> message related to LWORK is issued by XERBLA.
241 *> \endverbatim
242 *>
243 *> \param[out] RWORK
244 *> \verbatim
245 *> RWORK is DOUBLE PRECISION array, dimension (N)
246 *> \endverbatim
247 *>
248 *> \param[out] INFO
249 *> \verbatim
250 *> INFO is INTEGER
251 *> = 0: successful exit
252 *> < 0: if INFO = -i, the i-th argument had an illegal value
253 *> = 1,...,N: the QZ iteration did not converge. (H,T) is not
254 *> in Schur form, but ALPHA(i) and BETA(i),
255 *> i=INFO+1,...,N should be correct.
256 *> = N+1,...,2*N: the shift calculation failed. (H,T) is not
257 *> in Schur form, but ALPHA(i) and BETA(i),
258 *> i=INFO-N+1,...,N should be correct.
259 *> \endverbatim
260 *
261 * Authors:
262 * ========
263 *
264 *> \author Univ. of Tennessee
265 *> \author Univ. of California Berkeley
266 *> \author Univ. of Colorado Denver
267 *> \author NAG Ltd.
268 *
269 *> \date April 2012
270 *
271 *> \ingroup complex16GEcomputational
272 *
273 *> \par Further Details:
274 * =====================
275 *>
276 *> \verbatim
277 *>
278 *> We assume that complex ABS works as long as its value is less than
279 *> overflow.
280 *> \endverbatim
281 *>
282 * =====================================================================
283  SUBROUTINE zhgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
284  $ alpha, beta, q, ldq, z, ldz, work, lwork,
285  $ rwork, info )
286 *
287 * -- LAPACK computational routine (version 3.4.1) --
288 * -- LAPACK is a software package provided by Univ. of Tennessee, --
289 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290 * April 2012
291 *
292 * .. Scalar Arguments ..
293  CHARACTER compq, compz, job
294  INTEGER ihi, ilo, info, ldh, ldq, ldt, ldz, lwork, n
295 * ..
296 * .. Array Arguments ..
297  DOUBLE PRECISION rwork( * )
298  COMPLEX*16 alpha( * ), beta( * ), h( ldh, * ),
299  $ q( ldq, * ), t( ldt, * ), work( * ),
300  $ z( ldz, * )
301 * ..
302 *
303 * =====================================================================
304 *
305 * .. Parameters ..
306  COMPLEX*16 czero, cone
307  parameter( czero = ( 0.0d+0, 0.0d+0 ),
308  $ cone = ( 1.0d+0, 0.0d+0 ) )
309  DOUBLE PRECISION zero, one
310  parameter( zero = 0.0d+0, one = 1.0d+0 )
311  DOUBLE PRECISION half
312  parameter( half = 0.5d+0 )
313 * ..
314 * .. Local Scalars ..
315  LOGICAL ilazr2, ilazro, ilq, ilschr, ilz, lquery
316  INTEGER icompq, icompz, ifirst, ifrstm, iiter, ilast,
317  $ ilastm, in, ischur, istart, j, jc, jch, jiter,
318  $ jr, maxit
319  DOUBLE PRECISION absb, anorm, ascale, atol, bnorm, bscale, btol,
320  $ c, safmin, temp, temp2, tempr, ulp
321  COMPLEX*16 abi22, ad11, ad12, ad21, ad22, ctemp, ctemp2,
322  $ ctemp3, eshift, rtdisc, s, shift, signbc, t1,
323  $ u12, x
324 * ..
325 * .. External Functions ..
326  LOGICAL lsame
327  DOUBLE PRECISION dlamch, zlanhs
328  EXTERNAL lsame, dlamch, zlanhs
329 * ..
330 * .. External Subroutines ..
331  EXTERNAL xerbla, zlartg, zlaset, zrot, zscal
332 * ..
333 * .. Intrinsic Functions ..
334  INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min,
335  $ sqrt
336 * ..
337 * .. Statement Functions ..
338  DOUBLE PRECISION abs1
339 * ..
340 * .. Statement Function definitions ..
341  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
342 * ..
343 * .. Executable Statements ..
344 *
345 * Decode JOB, COMPQ, COMPZ
346 *
347  IF( lsame( job, 'E' ) ) THEN
348  ilschr = .false.
349  ischur = 1
350  ELSE IF( lsame( job, 'S' ) ) THEN
351  ilschr = .true.
352  ischur = 2
353  ELSE
354  ischur = 0
355  END IF
356 *
357  IF( lsame( compq, 'N' ) ) THEN
358  ilq = .false.
359  icompq = 1
360  ELSE IF( lsame( compq, 'V' ) ) THEN
361  ilq = .true.
362  icompq = 2
363  ELSE IF( lsame( compq, 'I' ) ) THEN
364  ilq = .true.
365  icompq = 3
366  ELSE
367  icompq = 0
368  END IF
369 *
370  IF( lsame( compz, 'N' ) ) THEN
371  ilz = .false.
372  icompz = 1
373  ELSE IF( lsame( compz, 'V' ) ) THEN
374  ilz = .true.
375  icompz = 2
376  ELSE IF( lsame( compz, 'I' ) ) THEN
377  ilz = .true.
378  icompz = 3
379  ELSE
380  icompz = 0
381  END IF
382 *
383 * Check Argument Values
384 *
385  info = 0
386  work( 1 ) = max( 1, n )
387  lquery = ( lwork.EQ.-1 )
388  IF( ischur.EQ.0 ) THEN
389  info = -1
390  ELSE IF( icompq.EQ.0 ) THEN
391  info = -2
392  ELSE IF( icompz.EQ.0 ) THEN
393  info = -3
394  ELSE IF( n.LT.0 ) THEN
395  info = -4
396  ELSE IF( ilo.LT.1 ) THEN
397  info = -5
398  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
399  info = -6
400  ELSE IF( ldh.LT.n ) THEN
401  info = -8
402  ELSE IF( ldt.LT.n ) THEN
403  info = -10
404  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
405  info = -14
406  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
407  info = -16
408  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
409  info = -18
410  END IF
411  IF( info.NE.0 ) THEN
412  CALL xerbla( 'ZHGEQZ', -info )
413  RETURN
414  ELSE IF( lquery ) THEN
415  RETURN
416  END IF
417 *
418 * Quick return if possible
419 *
420 * WORK( 1 ) = CMPLX( 1 )
421  IF( n.LE.0 ) THEN
422  work( 1 ) = dcmplx( 1 )
423  RETURN
424  END IF
425 *
426 * Initialize Q and Z
427 *
428  IF( icompq.EQ.3 )
429  $ CALL zlaset( 'Full', n, n, czero, cone, q, ldq )
430  IF( icompz.EQ.3 )
431  $ CALL zlaset( 'Full', n, n, czero, cone, z, ldz )
432 *
433 * Machine Constants
434 *
435  in = ihi + 1 - ilo
436  safmin = dlamch( 'S' )
437  ulp = dlamch( 'E' )*dlamch( 'B' )
438  anorm = zlanhs( 'F', in, h( ilo, ilo ), ldh, rwork )
439  bnorm = zlanhs( 'F', in, t( ilo, ilo ), ldt, rwork )
440  atol = max( safmin, ulp*anorm )
441  btol = max( safmin, ulp*bnorm )
442  ascale = one / max( safmin, anorm )
443  bscale = one / max( safmin, bnorm )
444 *
445 *
446 * Set Eigenvalues IHI+1:N
447 *
448  DO 10 j = ihi + 1, n
449  absb = abs( t( j, j ) )
450  IF( absb.GT.safmin ) THEN
451  signbc = dconjg( t( j, j ) / absb )
452  t( j, j ) = absb
453  IF( ilschr ) THEN
454  CALL zscal( j-1, signbc, t( 1, j ), 1 )
455  CALL zscal( j, signbc, h( 1, j ), 1 )
456  ELSE
457  h( j, j ) = h( j, j )*signbc
458  END IF
459  IF( ilz )
460  $ CALL zscal( n, signbc, z( 1, j ), 1 )
461  ELSE
462  t( j, j ) = czero
463  END IF
464  alpha( j ) = h( j, j )
465  beta( j ) = t( j, j )
466  10 CONTINUE
467 *
468 * If IHI < ILO, skip QZ steps
469 *
470  IF( ihi.LT.ilo )
471  $ go to 190
472 *
473 * MAIN QZ ITERATION LOOP
474 *
475 * Initialize dynamic indices
476 *
477 * Eigenvalues ILAST+1:N have been found.
478 * Column operations modify rows IFRSTM:whatever
479 * Row operations modify columns whatever:ILASTM
480 *
481 * If only eigenvalues are being computed, then
482 * IFRSTM is the row of the last splitting row above row ILAST;
483 * this is always at least ILO.
484 * IITER counts iterations since the last eigenvalue was found,
485 * to tell when to use an extraordinary shift.
486 * MAXIT is the maximum number of QZ sweeps allowed.
487 *
488  ilast = ihi
489  IF( ilschr ) THEN
490  ifrstm = 1
491  ilastm = n
492  ELSE
493  ifrstm = ilo
494  ilastm = ihi
495  END IF
496  iiter = 0
497  eshift = czero
498  maxit = 30*( ihi-ilo+1 )
499 *
500  DO 170 jiter = 1, maxit
501 *
502 * Check for too many iterations.
503 *
504  IF( jiter.GT.maxit )
505  $ go to 180
506 *
507 * Split the matrix if possible.
508 *
509 * Two tests:
510 * 1: H(j,j-1)=0 or j=ILO
511 * 2: T(j,j)=0
512 *
513 * Special case: j=ILAST
514 *
515  IF( ilast.EQ.ilo ) THEN
516  go to 60
517  ELSE
518  IF( abs1( h( ilast, ilast-1 ) ).LE.atol ) THEN
519  h( ilast, ilast-1 ) = czero
520  go to 60
521  END IF
522  END IF
523 *
524  IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
525  t( ilast, ilast ) = czero
526  go to 50
527  END IF
528 *
529 * General case: j<ILAST
530 *
531  DO 40 j = ilast - 1, ilo, -1
532 *
533 * Test 1: for H(j,j-1)=0 or j=ILO
534 *
535  IF( j.EQ.ilo ) THEN
536  ilazro = .true.
537  ELSE
538  IF( abs1( h( j, j-1 ) ).LE.atol ) THEN
539  h( j, j-1 ) = czero
540  ilazro = .true.
541  ELSE
542  ilazro = .false.
543  END IF
544  END IF
545 *
546 * Test 2: for T(j,j)=0
547 *
548  IF( abs( t( j, j ) ).LT.btol ) THEN
549  t( j, j ) = czero
550 *
551 * Test 1a: Check for 2 consecutive small subdiagonals in A
552 *
553  ilazr2 = .false.
554  IF( .NOT.ilazro ) THEN
555  IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
556  $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
557  $ ilazr2 = .true.
558  END IF
559 *
560 * If both tests pass (1 & 2), i.e., the leading diagonal
561 * element of B in the block is zero, split a 1x1 block off
562 * at the top. (I.e., at the J-th row/column) The leading
563 * diagonal element of the remainder can also be zero, so
564 * this may have to be done repeatedly.
565 *
566  IF( ilazro .OR. ilazr2 ) THEN
567  DO 20 jch = j, ilast - 1
568  ctemp = h( jch, jch )
569  CALL zlartg( ctemp, h( jch+1, jch ), c, s,
570  $ h( jch, jch ) )
571  h( jch+1, jch ) = czero
572  CALL zrot( ilastm-jch, h( jch, jch+1 ), ldh,
573  $ h( jch+1, jch+1 ), ldh, c, s )
574  CALL zrot( ilastm-jch, t( jch, jch+1 ), ldt,
575  $ t( jch+1, jch+1 ), ldt, c, s )
576  IF( ilq )
577  $ CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
578  $ c, dconjg( s ) )
579  IF( ilazr2 )
580  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
581  ilazr2 = .false.
582  IF( abs1( t( jch+1, jch+1 ) ).GE.btol ) THEN
583  IF( jch+1.GE.ilast ) THEN
584  go to 60
585  ELSE
586  ifirst = jch + 1
587  go to 70
588  END IF
589  END IF
590  t( jch+1, jch+1 ) = czero
591  20 CONTINUE
592  go to 50
593  ELSE
594 *
595 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
596 * Then process as in the case T(ILAST,ILAST)=0
597 *
598  DO 30 jch = j, ilast - 1
599  ctemp = t( jch, jch+1 )
600  CALL zlartg( ctemp, t( jch+1, jch+1 ), c, s,
601  $ t( jch, jch+1 ) )
602  t( jch+1, jch+1 ) = czero
603  IF( jch.LT.ilastm-1 )
604  $ CALL zrot( ilastm-jch-1, t( jch, jch+2 ), ldt,
605  $ t( jch+1, jch+2 ), ldt, c, s )
606  CALL zrot( ilastm-jch+2, h( jch, jch-1 ), ldh,
607  $ h( jch+1, jch-1 ), ldh, c, s )
608  IF( ilq )
609  $ CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
610  $ c, dconjg( s ) )
611  ctemp = h( jch+1, jch )
612  CALL zlartg( ctemp, h( jch+1, jch-1 ), c, s,
613  $ h( jch+1, jch ) )
614  h( jch+1, jch-1 ) = czero
615  CALL zrot( jch+1-ifrstm, h( ifrstm, jch ), 1,
616  $ h( ifrstm, jch-1 ), 1, c, s )
617  CALL zrot( jch-ifrstm, t( ifrstm, jch ), 1,
618  $ t( ifrstm, jch-1 ), 1, c, s )
619  IF( ilz )
620  $ CALL zrot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
621  $ c, s )
622  30 CONTINUE
623  go to 50
624  END IF
625  ELSE IF( ilazro ) THEN
626 *
627 * Only test 1 passed -- work on J:ILAST
628 *
629  ifirst = j
630  go to 70
631  END IF
632 *
633 * Neither test passed -- try next J
634 *
635  40 CONTINUE
636 *
637 * (Drop-through is "impossible")
638 *
639  info = 2*n + 1
640  go to 210
641 *
642 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
643 * 1x1 block.
644 *
645  50 CONTINUE
646  ctemp = h( ilast, ilast )
647  CALL zlartg( ctemp, h( ilast, ilast-1 ), c, s,
648  $ h( ilast, ilast ) )
649  h( ilast, ilast-1 ) = czero
650  CALL zrot( ilast-ifrstm, h( ifrstm, ilast ), 1,
651  $ h( ifrstm, ilast-1 ), 1, c, s )
652  CALL zrot( ilast-ifrstm, t( ifrstm, ilast ), 1,
653  $ t( ifrstm, ilast-1 ), 1, c, s )
654  IF( ilz )
655  $ CALL zrot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
656 *
657 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHA and BETA
658 *
659  60 CONTINUE
660  absb = abs( t( ilast, ilast ) )
661  IF( absb.GT.safmin ) THEN
662  signbc = dconjg( t( ilast, ilast ) / absb )
663  t( ilast, ilast ) = absb
664  IF( ilschr ) THEN
665  CALL zscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
666  CALL zscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
667  $ 1 )
668  ELSE
669  h( ilast, ilast ) = h( ilast, ilast )*signbc
670  END IF
671  IF( ilz )
672  $ CALL zscal( n, signbc, z( 1, ilast ), 1 )
673  ELSE
674  t( ilast, ilast ) = czero
675  END IF
676  alpha( ilast ) = h( ilast, ilast )
677  beta( ilast ) = t( ilast, ilast )
678 *
679 * Go to next block -- exit if finished.
680 *
681  ilast = ilast - 1
682  IF( ilast.LT.ilo )
683  $ go to 190
684 *
685 * Reset counters
686 *
687  iiter = 0
688  eshift = czero
689  IF( .NOT.ilschr ) THEN
690  ilastm = ilast
691  IF( ifrstm.GT.ilast )
692  $ ifrstm = ilo
693  END IF
694  go to 160
695 *
696 * QZ step
697 *
698 * This iteration only involves rows/columns IFIRST:ILAST. We
699 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
700 *
701  70 CONTINUE
702  iiter = iiter + 1
703  IF( .NOT.ilschr ) THEN
704  ifrstm = ifirst
705  END IF
706 *
707 * Compute the Shift.
708 *
709 * At this point, IFIRST < ILAST, and the diagonal elements of
710 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
711 * magnitude)
712 *
713  IF( ( iiter / 10 )*10.NE.iiter ) THEN
714 *
715 * The Wilkinson shift (AEP p.512), i.e., the eigenvalue of
716 * the bottom-right 2x2 block of A inv(B) which is nearest to
717 * the bottom-right element.
718 *
719 * We factor B as U*D, where U has unit diagonals, and
720 * compute (A*inv(D))*inv(U).
721 *
722  u12 = ( bscale*t( ilast-1, ilast ) ) /
723  $ ( bscale*t( ilast, ilast ) )
724  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
725  $ ( bscale*t( ilast-1, ilast-1 ) )
726  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
727  $ ( bscale*t( ilast-1, ilast-1 ) )
728  ad12 = ( ascale*h( ilast-1, ilast ) ) /
729  $ ( bscale*t( ilast, ilast ) )
730  ad22 = ( ascale*h( ilast, ilast ) ) /
731  $ ( bscale*t( ilast, ilast ) )
732  abi22 = ad22 - u12*ad21
733 *
734  t1 = half*( ad11+abi22 )
735  rtdisc = sqrt( t1**2+ad12*ad21-ad11*ad22 )
736  temp = dble( t1-abi22 )*dble( rtdisc ) +
737  $ dimag( t1-abi22 )*dimag( rtdisc )
738  IF( temp.LE.zero ) THEN
739  shift = t1 + rtdisc
740  ELSE
741  shift = t1 - rtdisc
742  END IF
743  ELSE
744 *
745 * Exceptional shift. Chosen for no particularly good reason.
746 *
747  eshift = eshift + (ascale*h(ilast,ilast-1))/
748  $ (bscale*t(ilast-1,ilast-1))
749  shift = eshift
750  END IF
751 *
752 * Now check for two consecutive small subdiagonals.
753 *
754  DO 80 j = ilast - 1, ifirst + 1, -1
755  istart = j
756  ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
757  temp = abs1( ctemp )
758  temp2 = ascale*abs1( h( j+1, j ) )
759  tempr = max( temp, temp2 )
760  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
761  temp = temp / tempr
762  temp2 = temp2 / tempr
763  END IF
764  IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
765  $ go to 90
766  80 CONTINUE
767 *
768  istart = ifirst
769  ctemp = ascale*h( ifirst, ifirst ) -
770  $ shift*( bscale*t( ifirst, ifirst ) )
771  90 CONTINUE
772 *
773 * Do an implicit-shift QZ sweep.
774 *
775 * Initial Q
776 *
777  ctemp2 = ascale*h( istart+1, istart )
778  CALL zlartg( ctemp, ctemp2, c, s, ctemp3 )
779 *
780 * Sweep
781 *
782  DO 150 j = istart, ilast - 1
783  IF( j.GT.istart ) THEN
784  ctemp = h( j, j-1 )
785  CALL zlartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
786  h( j+1, j-1 ) = czero
787  END IF
788 *
789  DO 100 jc = j, ilastm
790  ctemp = c*h( j, jc ) + s*h( j+1, jc )
791  h( j+1, jc ) = -dconjg( s )*h( j, jc ) + c*h( j+1, jc )
792  h( j, jc ) = ctemp
793  ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
794  t( j+1, jc ) = -dconjg( s )*t( j, jc ) + c*t( j+1, jc )
795  t( j, jc ) = ctemp2
796  100 CONTINUE
797  IF( ilq ) THEN
798  DO 110 jr = 1, n
799  ctemp = c*q( jr, j ) + dconjg( s )*q( jr, j+1 )
800  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
801  q( jr, j ) = ctemp
802  110 CONTINUE
803  END IF
804 *
805  ctemp = t( j+1, j+1 )
806  CALL zlartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
807  t( j+1, j ) = czero
808 *
809  DO 120 jr = ifrstm, min( j+2, ilast )
810  ctemp = c*h( jr, j+1 ) + s*h( jr, j )
811  h( jr, j ) = -dconjg( s )*h( jr, j+1 ) + c*h( jr, j )
812  h( jr, j+1 ) = ctemp
813  120 CONTINUE
814  DO 130 jr = ifrstm, j
815  ctemp = c*t( jr, j+1 ) + s*t( jr, j )
816  t( jr, j ) = -dconjg( s )*t( jr, j+1 ) + c*t( jr, j )
817  t( jr, j+1 ) = ctemp
818  130 CONTINUE
819  IF( ilz ) THEN
820  DO 140 jr = 1, n
821  ctemp = c*z( jr, j+1 ) + s*z( jr, j )
822  z( jr, j ) = -dconjg( s )*z( jr, j+1 ) + c*z( jr, j )
823  z( jr, j+1 ) = ctemp
824  140 CONTINUE
825  END IF
826  150 CONTINUE
827 *
828  160 CONTINUE
829 *
830  170 CONTINUE
831 *
832 * Drop-through = non-convergence
833 *
834  180 CONTINUE
835  info = ilast
836  go to 210
837 *
838 * Successful completion of all QZ steps
839 *
840  190 CONTINUE
841 *
842 * Set Eigenvalues 1:ILO-1
843 *
844  DO 200 j = 1, ilo - 1
845  absb = abs( t( j, j ) )
846  IF( absb.GT.safmin ) THEN
847  signbc = dconjg( t( j, j ) / absb )
848  t( j, j ) = absb
849  IF( ilschr ) THEN
850  CALL zscal( j-1, signbc, t( 1, j ), 1 )
851  CALL zscal( j, signbc, h( 1, j ), 1 )
852  ELSE
853  h( j, j ) = h( j, j )*signbc
854  END IF
855  IF( ilz )
856  $ CALL zscal( n, signbc, z( 1, j ), 1 )
857  ELSE
858  t( j, j ) = czero
859  END IF
860  alpha( j ) = h( j, j )
861  beta( j ) = t( j, j )
862  200 CONTINUE
863 *
864 * Normal Termination
865 *
866  info = 0
867 *
868 * Exit (other than argument error) -- return optimal workspace size
869 *
870  210 CONTINUE
871  work( 1 ) = dcmplx( n )
872  RETURN
873 *
874 * End of ZHGEQZ
875 *
876  END