LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 April 2012
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.4.1) --
308 * -- LAPACK is a software package provided by Univ. of Tennessee, --
309 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310 * April 2012
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-1, ilast ) ).LT.
743  $ abs( t( ilast-1, ilast-1 ) ) ) THEN
744  eshift = 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  temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
763  IF( wi.NE.zero )
764  $ go to 200
765  END IF
766 *
767 * Fiddle with shift to avoid overflow
768 *
769  temp = min( ascale, one )*( half*safmax )
770  IF( s1.GT.temp ) THEN
771  scale = temp / s1
772  ELSE
773  scale = one
774  END IF
775 *
776  temp = min( bscale, one )*( half*safmax )
777  IF( abs( wr ).GT.temp )
778  $ scale = min( scale, temp / abs( wr ) )
779  s1 = scale*s1
780  wr = scale*wr
781 *
782 * Now check for two consecutive small subdiagonals.
783 *
784  DO 120 j = ilast - 1, ifirst + 1, -1
785  istart = j
786  temp = abs( s1*h( j, j-1 ) )
787  temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
788  tempr = max( temp, temp2 )
789  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
790  temp = temp / tempr
791  temp2 = temp2 / tempr
792  END IF
793  IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
794  $ temp2 )go to 130
795  120 continue
796 *
797  istart = ifirst
798  130 continue
799 *
800 * Do an implicit single-shift QZ sweep.
801 *
802 * Initial Q
803 *
804  temp = s1*h( istart, istart ) - wr*t( istart, istart )
805  temp2 = s1*h( istart+1, istart )
806  CALL dlartg( temp, temp2, c, s, tempr )
807 *
808 * Sweep
809 *
810  DO 190 j = istart, ilast - 1
811  IF( j.GT.istart ) THEN
812  temp = h( j, j-1 )
813  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
814  h( j+1, j-1 ) = zero
815  END IF
816 *
817  DO 140 jc = j, ilastm
818  temp = c*h( j, jc ) + s*h( j+1, jc )
819  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
820  h( j, jc ) = temp
821  temp2 = c*t( j, jc ) + s*t( j+1, jc )
822  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
823  t( j, jc ) = temp2
824  140 continue
825  IF( ilq ) THEN
826  DO 150 jr = 1, n
827  temp = c*q( jr, j ) + s*q( jr, j+1 )
828  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
829  q( jr, j ) = temp
830  150 continue
831  END IF
832 *
833  temp = t( j+1, j+1 )
834  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
835  t( j+1, j ) = zero
836 *
837  DO 160 jr = ifrstm, min( j+2, ilast )
838  temp = c*h( jr, j+1 ) + s*h( jr, j )
839  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
840  h( jr, j+1 ) = temp
841  160 continue
842  DO 170 jr = ifrstm, j
843  temp = c*t( jr, j+1 ) + s*t( jr, j )
844  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
845  t( jr, j+1 ) = temp
846  170 continue
847  IF( ilz ) THEN
848  DO 180 jr = 1, n
849  temp = c*z( jr, j+1 ) + s*z( jr, j )
850  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
851  z( jr, j+1 ) = temp
852  180 continue
853  END IF
854  190 continue
855 *
856  go to 350
857 *
858 * Use Francis double-shift
859 *
860 * Note: the Francis double-shift should work with real shifts,
861 * but only if the block is at least 3x3.
862 * This code may break if this point is reached with
863 * a 2x2 block with real eigenvalues.
864 *
865  200 continue
866  IF( ifirst+1.EQ.ilast ) THEN
867 *
868 * Special case -- 2x2 block with complex eigenvectors
869 *
870 * Step 1: Standardize, that is, rotate so that
871 *
872 * ( B11 0 )
873 * B = ( ) with B11 non-negative.
874 * ( 0 B22 )
875 *
876  CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
877  $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
878 *
879  IF( b11.LT.zero ) THEN
880  cr = -cr
881  sr = -sr
882  b11 = -b11
883  b22 = -b22
884  END IF
885 *
886  CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
887  $ h( ilast, ilast-1 ), ldh, cl, sl )
888  CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
889  $ h( ifrstm, ilast ), 1, cr, sr )
890 *
891  IF( ilast.LT.ilastm )
892  $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
893  $ t( ilast, ilast+1 ), ldt, cl, sl )
894  IF( ifrstm.LT.ilast-1 )
895  $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
896  $ t( ifrstm, ilast ), 1, cr, sr )
897 *
898  IF( ilq )
899  $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
900  $ sl )
901  IF( ilz )
902  $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
903  $ sr )
904 *
905  t( ilast-1, ilast-1 ) = b11
906  t( ilast-1, ilast ) = zero
907  t( ilast, ilast-1 ) = zero
908  t( ilast, ilast ) = b22
909 *
910 * If B22 is negative, negate column ILAST
911 *
912  IF( b22.LT.zero ) THEN
913  DO 210 j = ifrstm, ilast
914  h( j, ilast ) = -h( j, ilast )
915  t( j, ilast ) = -t( j, ilast )
916  210 continue
917 *
918  IF( ilz ) THEN
919  DO 220 j = 1, n
920  z( j, ilast ) = -z( j, ilast )
921  220 continue
922  END IF
923  b22 = -b22
924  END IF
925 *
926 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
927 *
928 * Recompute shift
929 *
930  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
931  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
932  $ temp, wr, temp2, wi )
933 *
934 * If standardization has perturbed the shift onto real line,
935 * do another (real single-shift) QR step.
936 *
937  IF( wi.EQ.zero )
938  $ go to 350
939  s1inv = one / s1
940 *
941 * Do EISPACK (QZVAL) computation of alpha and beta
942 *
943  a11 = h( ilast-1, ilast-1 )
944  a21 = h( ilast, ilast-1 )
945  a12 = h( ilast-1, ilast )
946  a22 = h( ilast, ilast )
947 *
948 * Compute complex Givens rotation on right
949 * (Assume some element of C = (sA - wB) > unfl )
950 * __
951 * (sA - wB) ( CZ -SZ )
952 * ( SZ CZ )
953 *
954  c11r = s1*a11 - wr*b11
955  c11i = -wi*b11
956  c12 = s1*a12
957  c21 = s1*a21
958  c22r = s1*a22 - wr*b22
959  c22i = -wi*b22
960 *
961  IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
962  $ abs( c22r )+abs( c22i ) ) THEN
963  t1 = dlapy3( c12, c11r, c11i )
964  cz = c12 / t1
965  szr = -c11r / t1
966  szi = -c11i / t1
967  ELSE
968  cz = dlapy2( c22r, c22i )
969  IF( cz.LE.safmin ) THEN
970  cz = zero
971  szr = one
972  szi = zero
973  ELSE
974  tempr = c22r / cz
975  tempi = c22i / cz
976  t1 = dlapy2( cz, c21 )
977  cz = cz / t1
978  szr = -c21*tempr / t1
979  szi = c21*tempi / t1
980  END IF
981  END IF
982 *
983 * Compute Givens rotation on left
984 *
985 * ( CQ SQ )
986 * ( __ ) A or B
987 * ( -SQ CQ )
988 *
989  an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
990  bn = abs( b11 ) + abs( b22 )
991  wabs = abs( wr ) + abs( wi )
992  IF( s1*an.GT.wabs*bn ) THEN
993  cq = cz*b11
994  sqr = szr*b22
995  sqi = -szi*b22
996  ELSE
997  a1r = cz*a11 + szr*a12
998  a1i = szi*a12
999  a2r = cz*a21 + szr*a22
1000  a2i = szi*a22
1001  cq = dlapy2( a1r, a1i )
1002  IF( cq.LE.safmin ) THEN
1003  cq = zero
1004  sqr = one
1005  sqi = zero
1006  ELSE
1007  tempr = a1r / cq
1008  tempi = a1i / cq
1009  sqr = tempr*a2r + tempi*a2i
1010  sqi = tempi*a2r - tempr*a2i
1011  END IF
1012  END IF
1013  t1 = dlapy3( cq, sqr, sqi )
1014  cq = cq / t1
1015  sqr = sqr / t1
1016  sqi = sqi / t1
1017 *
1018 * Compute diagonal elements of QBZ
1019 *
1020  tempr = sqr*szr - sqi*szi
1021  tempi = sqr*szi + sqi*szr
1022  b1r = cq*cz*b11 + tempr*b22
1023  b1i = tempi*b22
1024  b1a = dlapy2( b1r, b1i )
1025  b2r = cq*cz*b22 + tempr*b11
1026  b2i = -tempi*b11
1027  b2a = dlapy2( b2r, b2i )
1028 *
1029 * Normalize so beta > 0, and Im( alpha1 ) > 0
1030 *
1031  beta( ilast-1 ) = b1a
1032  beta( ilast ) = b2a
1033  alphar( ilast-1 ) = ( wr*b1a )*s1inv
1034  alphai( ilast-1 ) = ( wi*b1a )*s1inv
1035  alphar( ilast ) = ( wr*b2a )*s1inv
1036  alphai( ilast ) = -( wi*b2a )*s1inv
1037 *
1038 * Step 3: Go to next block -- exit if finished.
1039 *
1040  ilast = ifirst - 1
1041  IF( ilast.LT.ilo )
1042  $ go to 380
1043 *
1044 * Reset counters
1045 *
1046  iiter = 0
1047  eshift = zero
1048  IF( .NOT.ilschr ) THEN
1049  ilastm = ilast
1050  IF( ifrstm.GT.ilast )
1051  $ ifrstm = ilo
1052  END IF
1053  go to 350
1054  ELSE
1055 *
1056 * Usual case: 3x3 or larger block, using Francis implicit
1057 * double-shift
1058 *
1059 * 2
1060 * Eigenvalue equation is w - c w + d = 0,
1061 *
1062 * -1 2 -1
1063 * so compute 1st column of (A B ) - c A B + d
1064 * using the formula in QZIT (from EISPACK)
1065 *
1066 * We assume that the block is at least 3x3
1067 *
1068  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1069  $ ( bscale*t( ilast-1, ilast-1 ) )
1070  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1071  $ ( bscale*t( ilast-1, ilast-1 ) )
1072  ad12 = ( ascale*h( ilast-1, ilast ) ) /
1073  $ ( bscale*t( ilast, ilast ) )
1074  ad22 = ( ascale*h( ilast, ilast ) ) /
1075  $ ( bscale*t( ilast, ilast ) )
1076  u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1077  ad11l = ( ascale*h( ifirst, ifirst ) ) /
1078  $ ( bscale*t( ifirst, ifirst ) )
1079  ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1080  $ ( bscale*t( ifirst, ifirst ) )
1081  ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1082  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1083  ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1084  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1085  ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1086  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1087  u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1088 *
1089  v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1090  $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1091  v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1092  $ ( ad22-ad11l )+ad21*u12 )*ad21l
1093  v( 3 ) = ad32l*ad21l
1094 *
1095  istart = ifirst
1096 *
1097  CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1098  v( 1 ) = one
1099 *
1100 * Sweep
1101 *
1102  DO 290 j = istart, ilast - 2
1103 *
1104 * All but last elements: use 3x3 Householder transforms.
1105 *
1106 * Zero (j-1)st column of A
1107 *
1108  IF( j.GT.istart ) THEN
1109  v( 1 ) = h( j, j-1 )
1110  v( 2 ) = h( j+1, j-1 )
1111  v( 3 ) = h( j+2, j-1 )
1112 *
1113  CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1114  v( 1 ) = one
1115  h( j+1, j-1 ) = zero
1116  h( j+2, j-1 ) = zero
1117  END IF
1118 *
1119  DO 230 jc = j, ilastm
1120  temp = tau*( h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1121  $ h( j+2, jc ) )
1122  h( j, jc ) = h( j, jc ) - temp
1123  h( j+1, jc ) = h( j+1, jc ) - temp*v( 2 )
1124  h( j+2, jc ) = h( j+2, jc ) - temp*v( 3 )
1125  temp2 = tau*( t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1126  $ t( j+2, jc ) )
1127  t( j, jc ) = t( j, jc ) - temp2
1128  t( j+1, jc ) = t( j+1, jc ) - temp2*v( 2 )
1129  t( j+2, jc ) = t( j+2, jc ) - temp2*v( 3 )
1130  230 continue
1131  IF( ilq ) THEN
1132  DO 240 jr = 1, n
1133  temp = tau*( q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1134  $ q( jr, j+2 ) )
1135  q( jr, j ) = q( jr, j ) - temp
1136  q( jr, j+1 ) = q( jr, j+1 ) - temp*v( 2 )
1137  q( jr, j+2 ) = q( jr, j+2 ) - temp*v( 3 )
1138  240 continue
1139  END IF
1140 *
1141 * Zero j-th column of B (see DLAGBC for details)
1142 *
1143 * Swap rows to pivot
1144 *
1145  ilpivt = .false.
1146  temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1147  temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1148  IF( max( temp, temp2 ).LT.safmin ) THEN
1149  scale = zero
1150  u1 = one
1151  u2 = zero
1152  go to 250
1153  ELSE IF( temp.GE.temp2 ) THEN
1154  w11 = t( j+1, j+1 )
1155  w21 = t( j+2, j+1 )
1156  w12 = t( j+1, j+2 )
1157  w22 = t( j+2, j+2 )
1158  u1 = t( j+1, j )
1159  u2 = t( j+2, j )
1160  ELSE
1161  w21 = t( j+1, j+1 )
1162  w11 = t( j+2, j+1 )
1163  w22 = t( j+1, j+2 )
1164  w12 = t( j+2, j+2 )
1165  u2 = t( j+1, j )
1166  u1 = t( j+2, j )
1167  END IF
1168 *
1169 * Swap columns if nec.
1170 *
1171  IF( abs( w12 ).GT.abs( w11 ) ) THEN
1172  ilpivt = .true.
1173  temp = w12
1174  temp2 = w22
1175  w12 = w11
1176  w22 = w21
1177  w11 = temp
1178  w21 = temp2
1179  END IF
1180 *
1181 * LU-factor
1182 *
1183  temp = w21 / w11
1184  u2 = u2 - temp*u1
1185  w22 = w22 - temp*w12
1186  w21 = zero
1187 *
1188 * Compute SCALE
1189 *
1190  scale = one
1191  IF( abs( w22 ).LT.safmin ) THEN
1192  scale = zero
1193  u2 = one
1194  u1 = -w12 / w11
1195  go to 250
1196  END IF
1197  IF( abs( w22 ).LT.abs( u2 ) )
1198  $ scale = abs( w22 / u2 )
1199  IF( abs( w11 ).LT.abs( u1 ) )
1200  $ scale = min( scale, abs( w11 / u1 ) )
1201 *
1202 * Solve
1203 *
1204  u2 = ( scale*u2 ) / w22
1205  u1 = ( scale*u1-w12*u2 ) / w11
1206 *
1207  250 continue
1208  IF( ilpivt ) THEN
1209  temp = u2
1210  u2 = u1
1211  u1 = temp
1212  END IF
1213 *
1214 * Compute Householder Vector
1215 *
1216  t1 = sqrt( scale**2+u1**2+u2**2 )
1217  tau = one + scale / t1
1218  vs = -one / ( scale+t1 )
1219  v( 1 ) = one
1220  v( 2 ) = vs*u1
1221  v( 3 ) = vs*u2
1222 *
1223 * Apply transformations from the right.
1224 *
1225  DO 260 jr = ifrstm, min( j+3, ilast )
1226  temp = tau*( h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1227  $ h( jr, j+2 ) )
1228  h( jr, j ) = h( jr, j ) - temp
1229  h( jr, j+1 ) = h( jr, j+1 ) - temp*v( 2 )
1230  h( jr, j+2 ) = h( jr, j+2 ) - temp*v( 3 )
1231  260 continue
1232  DO 270 jr = ifrstm, j + 2
1233  temp = tau*( t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1234  $ t( jr, j+2 ) )
1235  t( jr, j ) = t( jr, j ) - temp
1236  t( jr, j+1 ) = t( jr, j+1 ) - temp*v( 2 )
1237  t( jr, j+2 ) = t( jr, j+2 ) - temp*v( 3 )
1238  270 continue
1239  IF( ilz ) THEN
1240  DO 280 jr = 1, n
1241  temp = tau*( z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1242  $ z( jr, j+2 ) )
1243  z( jr, j ) = z( jr, j ) - temp
1244  z( jr, j+1 ) = z( jr, j+1 ) - temp*v( 2 )
1245  z( jr, j+2 ) = z( jr, j+2 ) - temp*v( 3 )
1246  280 continue
1247  END IF
1248  t( j+1, j ) = zero
1249  t( j+2, j ) = zero
1250  290 continue
1251 *
1252 * Last elements: Use Givens rotations
1253 *
1254 * Rotations from the left
1255 *
1256  j = ilast - 1
1257  temp = h( j, j-1 )
1258  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1259  h( j+1, j-1 ) = zero
1260 *
1261  DO 300 jc = j, ilastm
1262  temp = c*h( j, jc ) + s*h( j+1, jc )
1263  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1264  h( j, jc ) = temp
1265  temp2 = c*t( j, jc ) + s*t( j+1, jc )
1266  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1267  t( j, jc ) = temp2
1268  300 continue
1269  IF( ilq ) THEN
1270  DO 310 jr = 1, n
1271  temp = c*q( jr, j ) + s*q( jr, j+1 )
1272  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1273  q( jr, j ) = temp
1274  310 continue
1275  END IF
1276 *
1277 * Rotations from the right.
1278 *
1279  temp = t( j+1, j+1 )
1280  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1281  t( j+1, j ) = zero
1282 *
1283  DO 320 jr = ifrstm, ilast
1284  temp = c*h( jr, j+1 ) + s*h( jr, j )
1285  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1286  h( jr, j+1 ) = temp
1287  320 continue
1288  DO 330 jr = ifrstm, ilast - 1
1289  temp = c*t( jr, j+1 ) + s*t( jr, j )
1290  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1291  t( jr, j+1 ) = temp
1292  330 continue
1293  IF( ilz ) THEN
1294  DO 340 jr = 1, n
1295  temp = c*z( jr, j+1 ) + s*z( jr, j )
1296  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1297  z( jr, j+1 ) = temp
1298  340 continue
1299  END IF
1300 *
1301 * End of Double-Shift code
1302 *
1303  END IF
1304 *
1305  go to 350
1306 *
1307 * End of iteration loop
1308 *
1309  350 continue
1310  360 continue
1311 *
1312 * Drop-through = non-convergence
1313 *
1314  info = ilast
1315  go to 420
1316 *
1317 * Successful completion of all QZ steps
1318 *
1319  380 continue
1320 *
1321 * Set Eigenvalues 1:ILO-1
1322 *
1323  DO 410 j = 1, ilo - 1
1324  IF( t( j, j ).LT.zero ) THEN
1325  IF( ilschr ) THEN
1326  DO 390 jr = 1, j
1327  h( jr, j ) = -h( jr, j )
1328  t( jr, j ) = -t( jr, j )
1329  390 continue
1330  ELSE
1331  h( j, j ) = -h( j, j )
1332  t( j, j ) = -t( j, j )
1333  END IF
1334  IF( ilz ) THEN
1335  DO 400 jr = 1, n
1336  z( jr, j ) = -z( jr, j )
1337  400 continue
1338  END IF
1339  END IF
1340  alphar( j ) = h( j, j )
1341  alphai( j ) = zero
1342  beta( j ) = t( j, j )
1343  410 continue
1344 *
1345 * Normal Termination
1346 *
1347  info = 0
1348 *
1349 * Exit (other than argument error) -- return optimal workspace size
1350 *
1351  420 continue
1352  work( 1 ) = dble( n )
1353  return
1354 *
1355 * End of DHGEQZ
1356 *
1357  END