LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dhgeqz.f
Go to the documentation of this file.
1 *> \brief \b DHGEQZ
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DHGEQZ + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhgeqz.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhgeqz.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhgeqz.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
22 * ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
23 * LWORK, 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 ALPHAI( * ), ALPHAR( * ), BETA( * ),
31 * $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
32 * $ WORK( * ), Z( LDZ, * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
42 *> where H is an upper Hessenberg matrix and T is upper triangular,
43 *> using the double-shift QZ method.
44 *> Matrix pairs of this type are produced by the reduction to
45 *> generalized upper Hessenberg form of a real matrix pair (A,B):
46 *>
47 *> A = Q1*H*Z1**T, B = Q1*T*Z1**T,
48 *>
49 *> as computed by DGGHRD.
50 *>
51 *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
52 *> also reduced to generalized Schur form,
53 *>
54 *> H = Q*S*Z**T, T = Q*P*Z**T,
55 *>
56 *> where Q and Z are orthogonal matrices, P is an upper triangular
57 *> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
58 *> diagonal blocks.
59 *>
60 *> The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
61 *> (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
62 *> eigenvalues.
63 *>
64 *> Additionally, the 2-by-2 upper triangular diagonal blocks of P
65 *> corresponding to 2-by-2 blocks of S are reduced to positive diagonal
66 *> form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
67 *> P(j,j) > 0, and P(j+1,j+1) > 0.
68 *>
69 *> Optionally, the orthogonal matrix Q from the generalized Schur
70 *> factorization may be postmultiplied into an input matrix Q1, and the
71 *> orthogonal matrix Z may be postmultiplied into an input matrix Z1.
72 *> If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
73 *> the matrix pair (A,B) to generalized upper Hessenberg form, then the
74 *> output matrices Q1*Q and Z1*Z are the orthogonal factors from the
75 *> generalized Schur factorization of (A,B):
76 *>
77 *> A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T.
78 *>
79 *> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
80 *> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
81 *> complex and beta real.
82 *> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
83 *> generalized nonsymmetric eigenvalue problem (GNEP)
84 *> A*x = lambda*B*x
85 *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
86 *> alternate form of the GNEP
87 *> mu*A*y = B*y.
88 *> Real eigenvalues can be read directly from the generalized Schur
89 *> form:
90 *> alpha = S(i,i), beta = P(i,i).
91 *>
92 *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
93 *> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
94 *> pp. 241--256.
95 *> \endverbatim
96 *
97 * Arguments:
98 * ==========
99 *
100 *> \param[in] JOB
101 *> \verbatim
102 *> JOB is CHARACTER*1
103 *> = 'E': Compute eigenvalues only;
104 *> = 'S': Compute eigenvalues and the Schur form.
105 *> \endverbatim
106 *>
107 *> \param[in] COMPQ
108 *> \verbatim
109 *> COMPQ is CHARACTER*1
110 *> = 'N': Left Schur vectors (Q) are not computed;
111 *> = 'I': Q is initialized to the unit matrix and the matrix Q
112 *> of left Schur vectors of (H,T) is returned;
113 *> = 'V': Q must contain an orthogonal matrix Q1 on entry and
114 *> the product Q1*Q is returned.
115 *> \endverbatim
116 *>
117 *> \param[in] COMPZ
118 *> \verbatim
119 *> COMPZ is CHARACTER*1
120 *> = 'N': Right Schur vectors (Z) are not computed;
121 *> = 'I': Z is initialized to the unit matrix and the matrix Z
122 *> of right Schur vectors of (H,T) is returned;
123 *> = 'V': Z must contain an orthogonal matrix Z1 on entry and
124 *> the product Z1*Z is returned.
125 *> \endverbatim
126 *>
127 *> \param[in] N
128 *> \verbatim
129 *> N is INTEGER
130 *> The order of the matrices H, T, Q, and Z. N >= 0.
131 *> \endverbatim
132 *>
133 *> \param[in] ILO
134 *> \verbatim
135 *> ILO is INTEGER
136 *> \endverbatim
137 *>
138 *> \param[in] IHI
139 *> \verbatim
140 *> IHI is INTEGER
141 *> ILO and IHI mark the rows and columns of H which are in
142 *> Hessenberg form. It is assumed that A is already upper
143 *> triangular in rows and columns 1:ILO-1 and IHI+1:N.
144 *> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
145 *> \endverbatim
146 *>
147 *> \param[in,out] H
148 *> \verbatim
149 *> H is DOUBLE PRECISION array, dimension (LDH, N)
150 *> On entry, the N-by-N upper Hessenberg matrix H.
151 *> On exit, if JOB = 'S', H contains the upper quasi-triangular
152 *> matrix S from the generalized Schur factorization.
153 *> If JOB = 'E', the diagonal blocks of H match those of S, but
154 *> the rest of H is unspecified.
155 *> \endverbatim
156 *>
157 *> \param[in] LDH
158 *> \verbatim
159 *> LDH is INTEGER
160 *> The leading dimension of the array H. LDH >= max( 1, N ).
161 *> \endverbatim
162 *>
163 *> \param[in,out] T
164 *> \verbatim
165 *> T is DOUBLE PRECISION array, dimension (LDT, N)
166 *> On entry, the N-by-N upper triangular matrix T.
167 *> On exit, if JOB = 'S', T contains the upper triangular
168 *> matrix P from the generalized Schur factorization;
169 *> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
170 *> are reduced to positive diagonal form, i.e., if H(j+1,j) is
171 *> non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
172 *> T(j+1,j+1) > 0.
173 *> If JOB = 'E', the diagonal blocks of T match those of P, but
174 *> the rest of T is unspecified.
175 *> \endverbatim
176 *>
177 *> \param[in] LDT
178 *> \verbatim
179 *> LDT is INTEGER
180 *> The leading dimension of the array T. LDT >= max( 1, N ).
181 *> \endverbatim
182 *>
183 *> \param[out] ALPHAR
184 *> \verbatim
185 *> ALPHAR is DOUBLE PRECISION array, dimension (N)
186 *> The real parts of each scalar alpha defining an eigenvalue
187 *> of GNEP.
188 *> \endverbatim
189 *>
190 *> \param[out] ALPHAI
191 *> \verbatim
192 *> ALPHAI is DOUBLE PRECISION array, dimension (N)
193 *> The imaginary parts of each scalar alpha defining an
194 *> eigenvalue of GNEP.
195 *> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
196 *> positive, then the j-th and (j+1)-st eigenvalues are a
197 *> complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
198 *> \endverbatim
199 *>
200 *> \param[out] BETA
201 *> \verbatim
202 *> BETA is DOUBLE PRECISION array, dimension (N)
203 *> The scalars beta that define the eigenvalues of GNEP.
204 *> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
205 *> beta = BETA(j) represent the j-th eigenvalue of the matrix
206 *> pair (A,B), in one of the forms lambda = alpha/beta or
207 *> mu = beta/alpha. Since either lambda or mu may overflow,
208 *> they should not, in general, be computed.
209 *> \endverbatim
210 *>
211 *> \param[in,out] Q
212 *> \verbatim
213 *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
214 *> On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in
215 *> the reduction of (A,B) to generalized Hessenberg form.
216 *> On exit, if COMPZ = 'I', the orthogonal matrix of left Schur
217 *> vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix
218 *> of left Schur vectors of (A,B).
219 *> Not referenced if COMPZ = 'N'.
220 *> \endverbatim
221 *>
222 *> \param[in] LDQ
223 *> \verbatim
224 *> LDQ is INTEGER
225 *> The leading dimension of the array Q. LDQ >= 1.
226 *> If COMPQ='V' or 'I', then LDQ >= N.
227 *> \endverbatim
228 *>
229 *> \param[in,out] Z
230 *> \verbatim
231 *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
232 *> On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
233 *> the reduction of (A,B) to generalized Hessenberg form.
234 *> On exit, if COMPZ = 'I', the orthogonal matrix of
235 *> right Schur vectors of (H,T), and if COMPZ = 'V', the
236 *> orthogonal matrix of right Schur vectors of (A,B).
237 *> Not referenced if COMPZ = 'N'.
238 *> \endverbatim
239 *>
240 *> \param[in] LDZ
241 *> \verbatim
242 *> LDZ is INTEGER
243 *> The leading dimension of the array Z. LDZ >= 1.
244 *> If COMPZ='V' or 'I', then LDZ >= N.
245 *> \endverbatim
246 *>
247 *> \param[out] WORK
248 *> \verbatim
249 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
250 *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
251 *> \endverbatim
252 *>
253 *> \param[in] LWORK
254 *> \verbatim
255 *> LWORK is INTEGER
256 *> The dimension of the array WORK. LWORK >= max(1,N).
257 *>
258 *> If LWORK = -1, then a workspace query is assumed; the routine
259 *> only calculates the optimal size of the WORK array, returns
260 *> this value as the first entry of the WORK array, and no error
261 *> message related to LWORK is issued by XERBLA.
262 *> \endverbatim
263 *>
264 *> \param[out] INFO
265 *> \verbatim
266 *> INFO is INTEGER
267 *> = 0: successful exit
268 *> < 0: if INFO = -i, the i-th argument had an illegal value
269 *> = 1,...,N: the QZ iteration did not converge. (H,T) is not
270 *> in Schur form, but ALPHAR(i), ALPHAI(i), and
271 *> BETA(i), i=INFO+1,...,N should be correct.
272 *> = N+1,...,2*N: the shift calculation failed. (H,T) is not
273 *> in Schur form, but ALPHAR(i), ALPHAI(i), and
274 *> BETA(i), i=INFO-N+1,...,N should be correct.
275 *> \endverbatim
276 *
277 * Authors:
278 * ========
279 *
280 *> \author Univ. of Tennessee
281 *> \author Univ. of California Berkeley
282 *> \author Univ. of Colorado Denver
283 *> \author NAG Ltd.
284 *
285 *> \date November 2013
286 *
287 *> \ingroup doubleGEcomputational
288 *
289 *> \par Further Details:
290 * =====================
291 *>
292 *> \verbatim
293 *>
294 *> Iteration counters:
295 *>
296 *> JITER -- counts iterations.
297 *> IITER -- counts iterations run since ILAST was last
298 *> changed. This is therefore reset only when a 1-by-1 or
299 *> 2-by-2 block deflates off the bottom.
300 *> \endverbatim
301 *>
302 * =====================================================================
303  SUBROUTINE dhgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
304  $ alphar, alphai, beta, q, ldq, z, ldz, work,
305  $ lwork, info )
306 *
307 * -- LAPACK computational routine (version 3.5.0) --
308 * -- LAPACK is a software package provided by Univ. of Tennessee, --
309 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310 * November 2013
311 *
312 * .. Scalar Arguments ..
313  CHARACTER compq, compz, job
314  INTEGER ihi, ilo, info, ldh, ldq, ldt, ldz, lwork, n
315 * ..
316 * .. Array Arguments ..
317  DOUBLE PRECISION alphai( * ), alphar( * ), beta( * ),
318  $ h( ldh, * ), q( ldq, * ), t( ldt, * ),
319  $ work( * ), z( ldz, * )
320 * ..
321 *
322 * =====================================================================
323 *
324 * .. Parameters ..
325 * $ SAFETY = 1.0E+0 )
326  DOUBLE PRECISION half, zero, one, safety
327  parameter( half = 0.5d+0, zero = 0.0d+0, one = 1.0d+0,
328  $ safety = 1.0d+2 )
329 * ..
330 * .. Local Scalars ..
331  LOGICAL ilazr2, ilazro, ilpivt, ilq, ilschr, ilz,
332  $ lquery
333  INTEGER icompq, icompz, ifirst, ifrstm, iiter, ilast,
334  $ ilastm, in, ischur, istart, j, jc, jch, jiter,
335  $ jr, maxit
336  DOUBLE PRECISION a11, a12, a1i, a1r, a21, a22, a2i, a2r, ad11,
337  $ ad11l, ad12, ad12l, ad21, ad21l, ad22, ad22l,
338  $ ad32l, an, anorm, ascale, atol, b11, b1a, b1i,
339  $ b1r, b22, b2a, b2i, b2r, bn, bnorm, bscale,
340  $ btol, c, c11i, c11r, c12, c21, c22i, c22r, cl,
341  $ cq, cr, cz, eshift, s, s1, s1inv, s2, safmax,
342  $ safmin, scale, sl, sqi, sqr, sr, szi, szr, t1,
343  $ tau, temp, temp2, tempi, tempr, u1, u12, u12l,
344  $ u2, ulp, vs, w11, w12, w21, w22, wabs, wi, wr,
345  $ wr2
346 * ..
347 * .. Local Arrays ..
348  DOUBLE PRECISION v( 3 )
349 * ..
350 * .. External Functions ..
351  LOGICAL lsame
352  DOUBLE PRECISION dlamch, dlanhs, dlapy2, dlapy3
353  EXTERNAL lsame, dlamch, dlanhs, dlapy2, dlapy3
354 * ..
355 * .. External Subroutines ..
356  EXTERNAL dlag2, dlarfg, dlartg, dlaset, dlasv2, drot,
357  $ xerbla
358 * ..
359 * .. Intrinsic Functions ..
360  INTRINSIC abs, dble, max, min, sqrt
361 * ..
362 * .. Executable Statements ..
363 *
364 * Decode JOB, COMPQ, COMPZ
365 *
366  IF( lsame( job, 'E' ) ) THEN
367  ilschr = .false.
368  ischur = 1
369  ELSE IF( lsame( job, 'S' ) ) THEN
370  ilschr = .true.
371  ischur = 2
372  ELSE
373  ischur = 0
374  END IF
375 *
376  IF( lsame( compq, 'N' ) ) THEN
377  ilq = .false.
378  icompq = 1
379  ELSE IF( lsame( compq, 'V' ) ) THEN
380  ilq = .true.
381  icompq = 2
382  ELSE IF( lsame( compq, 'I' ) ) THEN
383  ilq = .true.
384  icompq = 3
385  ELSE
386  icompq = 0
387  END IF
388 *
389  IF( lsame( compz, 'N' ) ) THEN
390  ilz = .false.
391  icompz = 1
392  ELSE IF( lsame( compz, 'V' ) ) THEN
393  ilz = .true.
394  icompz = 2
395  ELSE IF( lsame( compz, 'I' ) ) THEN
396  ilz = .true.
397  icompz = 3
398  ELSE
399  icompz = 0
400  END IF
401 *
402 * Check Argument Values
403 *
404  info = 0
405  work( 1 ) = max( 1, n )
406  lquery = ( lwork.EQ.-1 )
407  IF( ischur.EQ.0 ) THEN
408  info = -1
409  ELSE IF( icompq.EQ.0 ) THEN
410  info = -2
411  ELSE IF( icompz.EQ.0 ) THEN
412  info = -3
413  ELSE IF( n.LT.0 ) THEN
414  info = -4
415  ELSE IF( ilo.LT.1 ) THEN
416  info = -5
417  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
418  info = -6
419  ELSE IF( ldh.LT.n ) THEN
420  info = -8
421  ELSE IF( ldt.LT.n ) THEN
422  info = -10
423  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
424  info = -15
425  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
426  info = -17
427  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
428  info = -19
429  END IF
430  IF( info.NE.0 ) THEN
431  CALL xerbla( 'DHGEQZ', -info )
432  RETURN
433  ELSE IF( lquery ) THEN
434  RETURN
435  END IF
436 *
437 * Quick return if possible
438 *
439  IF( n.LE.0 ) THEN
440  work( 1 ) = dble( 1 )
441  RETURN
442  END IF
443 *
444 * Initialize Q and Z
445 *
446  IF( icompq.EQ.3 )
447  $ CALL dlaset( 'Full', n, n, zero, one, q, ldq )
448  IF( icompz.EQ.3 )
449  $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
450 *
451 * Machine Constants
452 *
453  in = ihi + 1 - ilo
454  safmin = dlamch( 'S' )
455  safmax = one / safmin
456  ulp = dlamch( 'E' )*dlamch( 'B' )
457  anorm = dlanhs( 'F', in, h( ilo, ilo ), ldh, work )
458  bnorm = dlanhs( 'F', in, t( ilo, ilo ), ldt, work )
459  atol = max( safmin, ulp*anorm )
460  btol = max( safmin, ulp*bnorm )
461  ascale = one / max( safmin, anorm )
462  bscale = one / max( safmin, bnorm )
463 *
464 * Set Eigenvalues IHI+1:N
465 *
466  DO 30 j = ihi + 1, n
467  IF( t( j, j ).LT.zero ) THEN
468  IF( ilschr ) THEN
469  DO 10 jr = 1, j
470  h( jr, j ) = -h( jr, j )
471  t( jr, j ) = -t( jr, j )
472  10 CONTINUE
473  ELSE
474  h( j, j ) = -h( j, j )
475  t( j, j ) = -t( j, j )
476  END IF
477  IF( ilz ) THEN
478  DO 20 jr = 1, n
479  z( jr, j ) = -z( jr, j )
480  20 CONTINUE
481  END IF
482  END IF
483  alphar( j ) = h( j, j )
484  alphai( j ) = zero
485  beta( j ) = t( j, j )
486  30 CONTINUE
487 *
488 * If IHI < ILO, skip QZ steps
489 *
490  IF( ihi.LT.ilo )
491  $ go to 380
492 *
493 * MAIN QZ ITERATION LOOP
494 *
495 * Initialize dynamic indices
496 *
497 * Eigenvalues ILAST+1:N have been found.
498 * Column operations modify rows IFRSTM:whatever.
499 * Row operations modify columns whatever:ILASTM.
500 *
501 * If only eigenvalues are being computed, then
502 * IFRSTM is the row of the last splitting row above row ILAST;
503 * this is always at least ILO.
504 * IITER counts iterations since the last eigenvalue was found,
505 * to tell when to use an extraordinary shift.
506 * MAXIT is the maximum number of QZ sweeps allowed.
507 *
508  ilast = ihi
509  IF( ilschr ) THEN
510  ifrstm = 1
511  ilastm = n
512  ELSE
513  ifrstm = ilo
514  ilastm = ihi
515  END IF
516  iiter = 0
517  eshift = zero
518  maxit = 30*( ihi-ilo+1 )
519 *
520  DO 360 jiter = 1, maxit
521 *
522 * Split the matrix if possible.
523 *
524 * Two tests:
525 * 1: H(j,j-1)=0 or j=ILO
526 * 2: T(j,j)=0
527 *
528  IF( ilast.EQ.ilo ) THEN
529 *
530 * Special case: j=ILAST
531 *
532  go to 80
533  ELSE
534  IF( abs( h( ilast, ilast-1 ) ).LE.atol ) THEN
535  h( ilast, ilast-1 ) = zero
536  go to 80
537  END IF
538  END IF
539 *
540  IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
541  t( ilast, ilast ) = zero
542  go to 70
543  END IF
544 *
545 * General case: j<ILAST
546 *
547  DO 60 j = ilast - 1, ilo, -1
548 *
549 * Test 1: for H(j,j-1)=0 or j=ILO
550 *
551  IF( j.EQ.ilo ) THEN
552  ilazro = .true.
553  ELSE
554  IF( abs( h( j, j-1 ) ).LE.atol ) THEN
555  h( j, j-1 ) = zero
556  ilazro = .true.
557  ELSE
558  ilazro = .false.
559  END IF
560  END IF
561 *
562 * Test 2: for T(j,j)=0
563 *
564  IF( abs( t( j, j ) ).LT.btol ) THEN
565  t( j, j ) = zero
566 *
567 * Test 1a: Check for 2 consecutive small subdiagonals in A
568 *
569  ilazr2 = .false.
570  IF( .NOT.ilazro ) THEN
571  temp = abs( h( j, j-1 ) )
572  temp2 = abs( h( j, j ) )
573  tempr = max( temp, temp2 )
574  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
575  temp = temp / tempr
576  temp2 = temp2 / tempr
577  END IF
578  IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
579  $ ( ascale*atol ) )ilazr2 = .true.
580  END IF
581 *
582 * If both tests pass (1 & 2), i.e., the leading diagonal
583 * element of B in the block is zero, split a 1x1 block off
584 * at the top. (I.e., at the J-th row/column) The leading
585 * diagonal element of the remainder can also be zero, so
586 * this may have to be done repeatedly.
587 *
588  IF( ilazro .OR. ilazr2 ) THEN
589  DO 40 jch = j, ilast - 1
590  temp = h( jch, jch )
591  CALL dlartg( temp, h( jch+1, jch ), c, s,
592  $ h( jch, jch ) )
593  h( jch+1, jch ) = zero
594  CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
595  $ h( jch+1, jch+1 ), ldh, c, s )
596  CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
597  $ t( jch+1, jch+1 ), ldt, c, s )
598  IF( ilq )
599  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
600  $ c, s )
601  IF( ilazr2 )
602  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
603  ilazr2 = .false.
604  IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
605  IF( jch+1.GE.ilast ) THEN
606  go to 80
607  ELSE
608  ifirst = jch + 1
609  go to 110
610  END IF
611  END IF
612  t( jch+1, jch+1 ) = zero
613  40 CONTINUE
614  go to 70
615  ELSE
616 *
617 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
618 * Then process as in the case T(ILAST,ILAST)=0
619 *
620  DO 50 jch = j, ilast - 1
621  temp = t( jch, jch+1 )
622  CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
623  $ t( jch, jch+1 ) )
624  t( jch+1, jch+1 ) = zero
625  IF( jch.LT.ilastm-1 )
626  $ CALL drot( ilastm-jch-1, t( jch, jch+2 ), ldt,
627  $ t( jch+1, jch+2 ), ldt, c, s )
628  CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
629  $ h( jch+1, jch-1 ), ldh, c, s )
630  IF( ilq )
631  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
632  $ c, s )
633  temp = h( jch+1, jch )
634  CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
635  $ h( jch+1, jch ) )
636  h( jch+1, jch-1 ) = zero
637  CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
638  $ h( ifrstm, jch-1 ), 1, c, s )
639  CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
640  $ t( ifrstm, jch-1 ), 1, c, s )
641  IF( ilz )
642  $ CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
643  $ c, s )
644  50 CONTINUE
645  go to 70
646  END IF
647  ELSE IF( ilazro ) THEN
648 *
649 * Only test 1 passed -- work on J:ILAST
650 *
651  ifirst = j
652  go to 110
653  END IF
654 *
655 * Neither test passed -- try next J
656 *
657  60 CONTINUE
658 *
659 * (Drop-through is "impossible")
660 *
661  info = n + 1
662  go to 420
663 *
664 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
665 * 1x1 block.
666 *
667  70 CONTINUE
668  temp = h( ilast, ilast )
669  CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
670  $ h( ilast, ilast ) )
671  h( ilast, ilast-1 ) = zero
672  CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
673  $ h( ifrstm, ilast-1 ), 1, c, s )
674  CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
675  $ t( ifrstm, ilast-1 ), 1, c, s )
676  IF( ilz )
677  $ CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
678 *
679 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
680 * and BETA
681 *
682  80 CONTINUE
683  IF( t( ilast, ilast ).LT.zero ) THEN
684  IF( ilschr ) THEN
685  DO 90 j = ifrstm, ilast
686  h( j, ilast ) = -h( j, ilast )
687  t( j, ilast ) = -t( j, ilast )
688  90 CONTINUE
689  ELSE
690  h( ilast, ilast ) = -h( ilast, ilast )
691  t( ilast, ilast ) = -t( ilast, ilast )
692  END IF
693  IF( ilz ) THEN
694  DO 100 j = 1, n
695  z( j, ilast ) = -z( j, ilast )
696  100 CONTINUE
697  END IF
698  END IF
699  alphar( ilast ) = h( ilast, ilast )
700  alphai( ilast ) = zero
701  beta( ilast ) = t( ilast, ilast )
702 *
703 * Go to next block -- exit if finished.
704 *
705  ilast = ilast - 1
706  IF( ilast.LT.ilo )
707  $ go to 380
708 *
709 * Reset counters
710 *
711  iiter = 0
712  eshift = zero
713  IF( .NOT.ilschr ) THEN
714  ilastm = ilast
715  IF( ifrstm.GT.ilast )
716  $ ifrstm = ilo
717  END IF
718  go to 350
719 *
720 * QZ step
721 *
722 * This iteration only involves rows/columns IFIRST:ILAST. We
723 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
724 *
725  110 CONTINUE
726  iiter = iiter + 1
727  IF( .NOT.ilschr ) THEN
728  ifrstm = ifirst
729  END IF
730 *
731 * Compute single shifts.
732 *
733 * At this point, IFIRST < ILAST, and the diagonal elements of
734 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
735 * magnitude)
736 *
737  IF( ( iiter / 10 )*10.EQ.iiter ) THEN
738 *
739 * Exceptional shift. Chosen for no particularly good reason.
740 * (Single shift only.)
741 *
742  IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
743  $ abs( t( ilast-1, ilast-1 ) ) ) THEN
744  eshift = h( ilast, ilast-1 ) /
745  $ t( ilast-1, ilast-1 )
746  ELSE
747  eshift = eshift + one / ( safmin*dble( maxit ) )
748  END IF
749  s1 = one
750  wr = eshift
751 *
752  ELSE
753 *
754 * Shifts based on the generalized eigenvalues of the
755 * bottom-right 2x2 block of A and B. The first eigenvalue
756 * returned by DLAG2 is the Wilkinson shift (AEP p.512),
757 *
758  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
759  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
760  $ s2, wr, wr2, wi )
761 *
762  IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
763  $ .GT. abs( (wr2/s2)*t( ilast, ilast )
764  $ - h( ilast, ilast ) ) ) THEN
765  temp = wr
766  wr = wr2
767  wr2 = temp
768  temp = s1
769  s1 = s2
770  s2 = temp
771  END IF
772  temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
773  IF( wi.NE.zero )
774  $ go to 200
775  END IF
776 *
777 * Fiddle with shift to avoid overflow
778 *
779  temp = min( ascale, one )*( half*safmax )
780  IF( s1.GT.temp ) THEN
781  scale = temp / s1
782  ELSE
783  scale = one
784  END IF
785 *
786  temp = min( bscale, one )*( half*safmax )
787  IF( abs( wr ).GT.temp )
788  $ scale = min( scale, temp / abs( wr ) )
789  s1 = scale*s1
790  wr = scale*wr
791 *
792 * Now check for two consecutive small subdiagonals.
793 *
794  DO 120 j = ilast - 1, ifirst + 1, -1
795  istart = j
796  temp = abs( s1*h( j, j-1 ) )
797  temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
798  tempr = max( temp, temp2 )
799  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
800  temp = temp / tempr
801  temp2 = temp2 / tempr
802  END IF
803  IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
804  $ temp2 )go to 130
805  120 CONTINUE
806 *
807  istart = ifirst
808  130 CONTINUE
809 *
810 * Do an implicit single-shift QZ sweep.
811 *
812 * Initial Q
813 *
814  temp = s1*h( istart, istart ) - wr*t( istart, istart )
815  temp2 = s1*h( istart+1, istart )
816  CALL dlartg( temp, temp2, c, s, tempr )
817 *
818 * Sweep
819 *
820  DO 190 j = istart, ilast - 1
821  IF( j.GT.istart ) THEN
822  temp = h( j, j-1 )
823  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
824  h( j+1, j-1 ) = zero
825  END IF
826 *
827  DO 140 jc = j, ilastm
828  temp = c*h( j, jc ) + s*h( j+1, jc )
829  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
830  h( j, jc ) = temp
831  temp2 = c*t( j, jc ) + s*t( j+1, jc )
832  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
833  t( j, jc ) = temp2
834  140 CONTINUE
835  IF( ilq ) THEN
836  DO 150 jr = 1, n
837  temp = c*q( jr, j ) + s*q( jr, j+1 )
838  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
839  q( jr, j ) = temp
840  150 CONTINUE
841  END IF
842 *
843  temp = t( j+1, j+1 )
844  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
845  t( j+1, j ) = zero
846 *
847  DO 160 jr = ifrstm, min( j+2, ilast )
848  temp = c*h( jr, j+1 ) + s*h( jr, j )
849  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
850  h( jr, j+1 ) = temp
851  160 CONTINUE
852  DO 170 jr = ifrstm, j
853  temp = c*t( jr, j+1 ) + s*t( jr, j )
854  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
855  t( jr, j+1 ) = temp
856  170 CONTINUE
857  IF( ilz ) THEN
858  DO 180 jr = 1, n
859  temp = c*z( jr, j+1 ) + s*z( jr, j )
860  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
861  z( jr, j+1 ) = temp
862  180 CONTINUE
863  END IF
864  190 CONTINUE
865 *
866  go to 350
867 *
868 * Use Francis double-shift
869 *
870 * Note: the Francis double-shift should work with real shifts,
871 * but only if the block is at least 3x3.
872 * This code may break if this point is reached with
873 * a 2x2 block with real eigenvalues.
874 *
875  200 CONTINUE
876  IF( ifirst+1.EQ.ilast ) THEN
877 *
878 * Special case -- 2x2 block with complex eigenvectors
879 *
880 * Step 1: Standardize, that is, rotate so that
881 *
882 * ( B11 0 )
883 * B = ( ) with B11 non-negative.
884 * ( 0 B22 )
885 *
886  CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
887  $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
888 *
889  IF( b11.LT.zero ) THEN
890  cr = -cr
891  sr = -sr
892  b11 = -b11
893  b22 = -b22
894  END IF
895 *
896  CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
897  $ h( ilast, ilast-1 ), ldh, cl, sl )
898  CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
899  $ h( ifrstm, ilast ), 1, cr, sr )
900 *
901  IF( ilast.LT.ilastm )
902  $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
903  $ t( ilast, ilast+1 ), ldt, cl, sl )
904  IF( ifrstm.LT.ilast-1 )
905  $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
906  $ t( ifrstm, ilast ), 1, cr, sr )
907 *
908  IF( ilq )
909  $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
910  $ sl )
911  IF( ilz )
912  $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
913  $ sr )
914 *
915  t( ilast-1, ilast-1 ) = b11
916  t( ilast-1, ilast ) = zero
917  t( ilast, ilast-1 ) = zero
918  t( ilast, ilast ) = b22
919 *
920 * If B22 is negative, negate column ILAST
921 *
922  IF( b22.LT.zero ) THEN
923  DO 210 j = ifrstm, ilast
924  h( j, ilast ) = -h( j, ilast )
925  t( j, ilast ) = -t( j, ilast )
926  210 CONTINUE
927 *
928  IF( ilz ) THEN
929  DO 220 j = 1, n
930  z( j, ilast ) = -z( j, ilast )
931  220 CONTINUE
932  END IF
933  b22 = -b22
934  END IF
935 *
936 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
937 *
938 * Recompute shift
939 *
940  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
941  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
942  $ temp, wr, temp2, wi )
943 *
944 * If standardization has perturbed the shift onto real line,
945 * do another (real single-shift) QR step.
946 *
947  IF( wi.EQ.zero )
948  $ go to 350
949  s1inv = one / s1
950 *
951 * Do EISPACK (QZVAL) computation of alpha and beta
952 *
953  a11 = h( ilast-1, ilast-1 )
954  a21 = h( ilast, ilast-1 )
955  a12 = h( ilast-1, ilast )
956  a22 = h( ilast, ilast )
957 *
958 * Compute complex Givens rotation on right
959 * (Assume some element of C = (sA - wB) > unfl )
960 * __
961 * (sA - wB) ( CZ -SZ )
962 * ( SZ CZ )
963 *
964  c11r = s1*a11 - wr*b11
965  c11i = -wi*b11
966  c12 = s1*a12
967  c21 = s1*a21
968  c22r = s1*a22 - wr*b22
969  c22i = -wi*b22
970 *
971  IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
972  $ abs( c22r )+abs( c22i ) ) THEN
973  t1 = dlapy3( c12, c11r, c11i )
974  cz = c12 / t1
975  szr = -c11r / t1
976  szi = -c11i / t1
977  ELSE
978  cz = dlapy2( c22r, c22i )
979  IF( cz.LE.safmin ) THEN
980  cz = zero
981  szr = one
982  szi = zero
983  ELSE
984  tempr = c22r / cz
985  tempi = c22i / cz
986  t1 = dlapy2( cz, c21 )
987  cz = cz / t1
988  szr = -c21*tempr / t1
989  szi = c21*tempi / t1
990  END IF
991  END IF
992 *
993 * Compute Givens rotation on left
994 *
995 * ( CQ SQ )
996 * ( __ ) A or B
997 * ( -SQ CQ )
998 *
999  an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1000  bn = abs( b11 ) + abs( b22 )
1001  wabs = abs( wr ) + abs( wi )
1002  IF( s1*an.GT.wabs*bn ) THEN
1003  cq = cz*b11
1004  sqr = szr*b22
1005  sqi = -szi*b22
1006  ELSE
1007  a1r = cz*a11 + szr*a12
1008  a1i = szi*a12
1009  a2r = cz*a21 + szr*a22
1010  a2i = szi*a22
1011  cq = dlapy2( a1r, a1i )
1012  IF( cq.LE.safmin ) THEN
1013  cq = zero
1014  sqr = one
1015  sqi = zero
1016  ELSE
1017  tempr = a1r / cq
1018  tempi = a1i / cq
1019  sqr = tempr*a2r + tempi*a2i
1020  sqi = tempi*a2r - tempr*a2i
1021  END IF
1022  END IF
1023  t1 = dlapy3( cq, sqr, sqi )
1024  cq = cq / t1
1025  sqr = sqr / t1
1026  sqi = sqi / t1
1027 *
1028 * Compute diagonal elements of QBZ
1029 *
1030  tempr = sqr*szr - sqi*szi
1031  tempi = sqr*szi + sqi*szr
1032  b1r = cq*cz*b11 + tempr*b22
1033  b1i = tempi*b22
1034  b1a = dlapy2( b1r, b1i )
1035  b2r = cq*cz*b22 + tempr*b11
1036  b2i = -tempi*b11
1037  b2a = dlapy2( b2r, b2i )
1038 *
1039 * Normalize so beta > 0, and Im( alpha1 ) > 0
1040 *
1041  beta( ilast-1 ) = b1a
1042  beta( ilast ) = b2a
1043  alphar( ilast-1 ) = ( wr*b1a )*s1inv
1044  alphai( ilast-1 ) = ( wi*b1a )*s1inv
1045  alphar( ilast ) = ( wr*b2a )*s1inv
1046  alphai( ilast ) = -( wi*b2a )*s1inv
1047 *
1048 * Step 3: Go to next block -- exit if finished.
1049 *
1050  ilast = ifirst - 1
1051  IF( ilast.LT.ilo )
1052  $ go to 380
1053 *
1054 * Reset counters
1055 *
1056  iiter = 0
1057  eshift = zero
1058  IF( .NOT.ilschr ) THEN
1059  ilastm = ilast
1060  IF( ifrstm.GT.ilast )
1061  $ ifrstm = ilo
1062  END IF
1063  go to 350
1064  ELSE
1065 *
1066 * Usual case: 3x3 or larger block, using Francis implicit
1067 * double-shift
1068 *
1069 * 2
1070 * Eigenvalue equation is w - c w + d = 0,
1071 *
1072 * -1 2 -1
1073 * so compute 1st column of (A B ) - c A B + d
1074 * using the formula in QZIT (from EISPACK)
1075 *
1076 * We assume that the block is at least 3x3
1077 *
1078  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1079  $ ( bscale*t( ilast-1, ilast-1 ) )
1080  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1081  $ ( bscale*t( ilast-1, ilast-1 ) )
1082  ad12 = ( ascale*h( ilast-1, ilast ) ) /
1083  $ ( bscale*t( ilast, ilast ) )
1084  ad22 = ( ascale*h( ilast, ilast ) ) /
1085  $ ( bscale*t( ilast, ilast ) )
1086  u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1087  ad11l = ( ascale*h( ifirst, ifirst ) ) /
1088  $ ( bscale*t( ifirst, ifirst ) )
1089  ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1090  $ ( bscale*t( ifirst, ifirst ) )
1091  ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1092  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1093  ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1094  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1095  ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1096  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1097  u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1098 *
1099  v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1100  $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1101  v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1102  $ ( ad22-ad11l )+ad21*u12 )*ad21l
1103  v( 3 ) = ad32l*ad21l
1104 *
1105  istart = ifirst
1106 *
1107  CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1108  v( 1 ) = one
1109 *
1110 * Sweep
1111 *
1112  DO 290 j = istart, ilast - 2
1113 *
1114 * All but last elements: use 3x3 Householder transforms.
1115 *
1116 * Zero (j-1)st column of A
1117 *
1118  IF( j.GT.istart ) THEN
1119  v( 1 ) = h( j, j-1 )
1120  v( 2 ) = h( j+1, j-1 )
1121  v( 3 ) = h( j+2, j-1 )
1122 *
1123  CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1124  v( 1 ) = one
1125  h( j+1, j-1 ) = zero
1126  h( j+2, j-1 ) = zero
1127  END IF
1128 *
1129  DO 230 jc = j, ilastm
1130  temp = tau*( h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1131  $ h( j+2, jc ) )
1132  h( j, jc ) = h( j, jc ) - temp
1133  h( j+1, jc ) = h( j+1, jc ) - temp*v( 2 )
1134  h( j+2, jc ) = h( j+2, jc ) - temp*v( 3 )
1135  temp2 = tau*( t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1136  $ t( j+2, jc ) )
1137  t( j, jc ) = t( j, jc ) - temp2
1138  t( j+1, jc ) = t( j+1, jc ) - temp2*v( 2 )
1139  t( j+2, jc ) = t( j+2, jc ) - temp2*v( 3 )
1140  230 CONTINUE
1141  IF( ilq ) THEN
1142  DO 240 jr = 1, n
1143  temp = tau*( q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1144  $ q( jr, j+2 ) )
1145  q( jr, j ) = q( jr, j ) - temp
1146  q( jr, j+1 ) = q( jr, j+1 ) - temp*v( 2 )
1147  q( jr, j+2 ) = q( jr, j+2 ) - temp*v( 3 )
1148  240 CONTINUE
1149  END IF
1150 *
1151 * Zero j-th column of B (see DLAGBC for details)
1152 *
1153 * Swap rows to pivot
1154 *
1155  ilpivt = .false.
1156  temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1157  temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1158  IF( max( temp, temp2 ).LT.safmin ) THEN
1159  scale = zero
1160  u1 = one
1161  u2 = zero
1162  go to 250
1163  ELSE IF( temp.GE.temp2 ) THEN
1164  w11 = t( j+1, j+1 )
1165  w21 = t( j+2, j+1 )
1166  w12 = t( j+1, j+2 )
1167  w22 = t( j+2, j+2 )
1168  u1 = t( j+1, j )
1169  u2 = t( j+2, j )
1170  ELSE
1171  w21 = t( j+1, j+1 )
1172  w11 = t( j+2, j+1 )
1173  w22 = t( j+1, j+2 )
1174  w12 = t( j+2, j+2 )
1175  u2 = t( j+1, j )
1176  u1 = t( j+2, j )
1177  END IF
1178 *
1179 * Swap columns if nec.
1180 *
1181  IF( abs( w12 ).GT.abs( w11 ) ) THEN
1182  ilpivt = .true.
1183  temp = w12
1184  temp2 = w22
1185  w12 = w11
1186  w22 = w21
1187  w11 = temp
1188  w21 = temp2
1189  END IF
1190 *
1191 * LU-factor
1192 *
1193  temp = w21 / w11
1194  u2 = u2 - temp*u1
1195  w22 = w22 - temp*w12
1196  w21 = zero
1197 *
1198 * Compute SCALE
1199 *
1200  scale = one
1201  IF( abs( w22 ).LT.safmin ) THEN
1202  scale = zero
1203  u2 = one
1204  u1 = -w12 / w11
1205  go to 250
1206  END IF
1207  IF( abs( w22 ).LT.abs( u2 ) )
1208  $ scale = abs( w22 / u2 )
1209  IF( abs( w11 ).LT.abs( u1 ) )
1210  $ scale = min( scale, abs( w11 / u1 ) )
1211 *
1212 * Solve
1213 *
1214  u2 = ( scale*u2 ) / w22
1215  u1 = ( scale*u1-w12*u2 ) / w11
1216 *
1217  250 CONTINUE
1218  IF( ilpivt ) THEN
1219  temp = u2
1220  u2 = u1
1221  u1 = temp
1222  END IF
1223 *
1224 * Compute Householder Vector
1225 *
1226  t1 = sqrt( scale**2+u1**2+u2**2 )
1227  tau = one + scale / t1
1228  vs = -one / ( scale+t1 )
1229  v( 1 ) = one
1230  v( 2 ) = vs*u1
1231  v( 3 ) = vs*u2
1232 *
1233 * Apply transformations from the right.
1234 *
1235  DO 260 jr = ifrstm, min( j+3, ilast )
1236  temp = tau*( h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1237  $ h( jr, j+2 ) )
1238  h( jr, j ) = h( jr, j ) - temp
1239  h( jr, j+1 ) = h( jr, j+1 ) - temp*v( 2 )
1240  h( jr, j+2 ) = h( jr, j+2 ) - temp*v( 3 )
1241  260 CONTINUE
1242  DO 270 jr = ifrstm, j + 2
1243  temp = tau*( t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1244  $ t( jr, j+2 ) )
1245  t( jr, j ) = t( jr, j ) - temp
1246  t( jr, j+1 ) = t( jr, j+1 ) - temp*v( 2 )
1247  t( jr, j+2 ) = t( jr, j+2 ) - temp*v( 3 )
1248  270 CONTINUE
1249  IF( ilz ) THEN
1250  DO 280 jr = 1, n
1251  temp = tau*( z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1252  $ z( jr, j+2 ) )
1253  z( jr, j ) = z( jr, j ) - temp
1254  z( jr, j+1 ) = z( jr, j+1 ) - temp*v( 2 )
1255  z( jr, j+2 ) = z( jr, j+2 ) - temp*v( 3 )
1256  280 CONTINUE
1257  END IF
1258  t( j+1, j ) = zero
1259  t( j+2, j ) = zero
1260  290 CONTINUE
1261 *
1262 * Last elements: Use Givens rotations
1263 *
1264 * Rotations from the left
1265 *
1266  j = ilast - 1
1267  temp = h( j, j-1 )
1268  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1269  h( j+1, j-1 ) = zero
1270 *
1271  DO 300 jc = j, ilastm
1272  temp = c*h( j, jc ) + s*h( j+1, jc )
1273  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1274  h( j, jc ) = temp
1275  temp2 = c*t( j, jc ) + s*t( j+1, jc )
1276  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1277  t( j, jc ) = temp2
1278  300 CONTINUE
1279  IF( ilq ) THEN
1280  DO 310 jr = 1, n
1281  temp = c*q( jr, j ) + s*q( jr, j+1 )
1282  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1283  q( jr, j ) = temp
1284  310 CONTINUE
1285  END IF
1286 *
1287 * Rotations from the right.
1288 *
1289  temp = t( j+1, j+1 )
1290  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1291  t( j+1, j ) = zero
1292 *
1293  DO 320 jr = ifrstm, ilast
1294  temp = c*h( jr, j+1 ) + s*h( jr, j )
1295  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1296  h( jr, j+1 ) = temp
1297  320 CONTINUE
1298  DO 330 jr = ifrstm, ilast - 1
1299  temp = c*t( jr, j+1 ) + s*t( jr, j )
1300  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1301  t( jr, j+1 ) = temp
1302  330 CONTINUE
1303  IF( ilz ) THEN
1304  DO 340 jr = 1, n
1305  temp = c*z( jr, j+1 ) + s*z( jr, j )
1306  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1307  z( jr, j+1 ) = temp
1308  340 CONTINUE
1309  END IF
1310 *
1311 * End of Double-Shift code
1312 *
1313  END IF
1314 *
1315  go to 350
1316 *
1317 * End of iteration loop
1318 *
1319  350 CONTINUE
1320  360 CONTINUE
1321 *
1322 * Drop-through = non-convergence
1323 *
1324  info = ilast
1325  go to 420
1326 *
1327 * Successful completion of all QZ steps
1328 *
1329  380 CONTINUE
1330 *
1331 * Set Eigenvalues 1:ILO-1
1332 *
1333  DO 410 j = 1, ilo - 1
1334  IF( t( j, j ).LT.zero ) THEN
1335  IF( ilschr ) THEN
1336  DO 390 jr = 1, j
1337  h( jr, j ) = -h( jr, j )
1338  t( jr, j ) = -t( jr, j )
1339  390 CONTINUE
1340  ELSE
1341  h( j, j ) = -h( j, j )
1342  t( j, j ) = -t( j, j )
1343  END IF
1344  IF( ilz ) THEN
1345  DO 400 jr = 1, n
1346  z( jr, j ) = -z( jr, j )
1347  400 CONTINUE
1348  END IF
1349  END IF
1350  alphar( j ) = h( j, j )
1351  alphai( j ) = zero
1352  beta( j ) = t( j, j )
1353  410 CONTINUE
1354 *
1355 * Normal Termination
1356 *
1357  info = 0
1358 *
1359 * Exit (other than argument error) -- return optimal workspace size
1360 *
1361  420 CONTINUE
1362  work( 1 ) = dble( n )
1363  RETURN
1364 *
1365 * End of DHGEQZ
1366 *
1367  END