LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 COMPQ = 'V', the orthogonal matrix Q1 used in
215 *> the reduction of (A,B) to generalized Hessenberg form.
216 *> On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
217 *> vectors of (H,T), and if COMPQ = 'V', the orthogonal matrix
218 *> of left Schur vectors of (A,B).
219 *> Not referenced if COMPQ = '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 *> \ingroup doubleGEcomputational
286 *
287 *> \par Further Details:
288 * =====================
289 *>
290 *> \verbatim
291 *>
292 *> Iteration counters:
293 *>
294 *> JITER -- counts iterations.
295 *> IITER -- counts iterations run since ILAST was last
296 *> changed. This is therefore reset only when a 1-by-1 or
297 *> 2-by-2 block deflates off the bottom.
298 *> \endverbatim
299 *>
300 * =====================================================================
301  SUBROUTINE dhgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
302  $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
303  $ LWORK, INFO )
304 *
305 * -- LAPACK computational routine --
306 * -- LAPACK is a software package provided by Univ. of Tennessee, --
307 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
308 *
309 * .. Scalar Arguments ..
310  CHARACTER COMPQ, COMPZ, JOB
311  INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
312 * ..
313 * .. Array Arguments ..
314  DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ),
315  $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
316  $ work( * ), z( ldz, * )
317 * ..
318 *
319 * =====================================================================
320 *
321 * .. Parameters ..
322 * $ SAFETY = 1.0E+0 )
323  DOUBLE PRECISION HALF, ZERO, ONE, SAFETY
324  PARAMETER ( HALF = 0.5d+0, zero = 0.0d+0, one = 1.0d+0,
325  $ safety = 1.0d+2 )
326 * ..
327 * .. Local Scalars ..
328  LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
329  $ LQUERY
330  INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
331  $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
332  $ jr, maxit
333  DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
334  $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
335  $ ad32l, an, anorm, ascale, atol, b11, b1a, b1i,
336  $ b1r, b22, b2a, b2i, b2r, bn, bnorm, bscale,
337  $ btol, c, c11i, c11r, c12, c21, c22i, c22r, cl,
338  $ cq, cr, cz, eshift, s, s1, s1inv, s2, safmax,
339  $ safmin, scale, sl, sqi, sqr, sr, szi, szr, t1,
340  $ tau, temp, temp2, tempi, tempr, u1, u12, u12l,
341  $ u2, ulp, vs, w11, w12, w21, w22, wabs, wi, wr,
342  $ wr2
343 * ..
344 * .. Local Arrays ..
345  DOUBLE PRECISION V( 3 )
346 * ..
347 * .. External Functions ..
348  LOGICAL LSAME
349  DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3
350  EXTERNAL lsame, dlamch, dlanhs, dlapy2, dlapy3
351 * ..
352 * .. External Subroutines ..
353  EXTERNAL dlag2, dlarfg, dlartg, dlaset, dlasv2, drot,
354  $ xerbla
355 * ..
356 * .. Intrinsic Functions ..
357  INTRINSIC abs, dble, max, min, sqrt
358 * ..
359 * .. Executable Statements ..
360 *
361 * Decode JOB, COMPQ, COMPZ
362 *
363  IF( lsame( job, 'E' ) ) THEN
364  ilschr = .false.
365  ischur = 1
366  ELSE IF( lsame( job, 'S' ) ) THEN
367  ilschr = .true.
368  ischur = 2
369  ELSE
370  ischur = 0
371  END IF
372 *
373  IF( lsame( compq, 'N' ) ) THEN
374  ilq = .false.
375  icompq = 1
376  ELSE IF( lsame( compq, 'V' ) ) THEN
377  ilq = .true.
378  icompq = 2
379  ELSE IF( lsame( compq, 'I' ) ) THEN
380  ilq = .true.
381  icompq = 3
382  ELSE
383  icompq = 0
384  END IF
385 *
386  IF( lsame( compz, 'N' ) ) THEN
387  ilz = .false.
388  icompz = 1
389  ELSE IF( lsame( compz, 'V' ) ) THEN
390  ilz = .true.
391  icompz = 2
392  ELSE IF( lsame( compz, 'I' ) ) THEN
393  ilz = .true.
394  icompz = 3
395  ELSE
396  icompz = 0
397  END IF
398 *
399 * Check Argument Values
400 *
401  info = 0
402  work( 1 ) = max( 1, n )
403  lquery = ( lwork.EQ.-1 )
404  IF( ischur.EQ.0 ) THEN
405  info = -1
406  ELSE IF( icompq.EQ.0 ) THEN
407  info = -2
408  ELSE IF( icompz.EQ.0 ) THEN
409  info = -3
410  ELSE IF( n.LT.0 ) THEN
411  info = -4
412  ELSE IF( ilo.LT.1 ) THEN
413  info = -5
414  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
415  info = -6
416  ELSE IF( ldh.LT.n ) THEN
417  info = -8
418  ELSE IF( ldt.LT.n ) THEN
419  info = -10
420  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
421  info = -15
422  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
423  info = -17
424  ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
425  info = -19
426  END IF
427  IF( info.NE.0 ) THEN
428  CALL xerbla( 'DHGEQZ', -info )
429  RETURN
430  ELSE IF( lquery ) THEN
431  RETURN
432  END IF
433 *
434 * Quick return if possible
435 *
436  IF( n.LE.0 ) THEN
437  work( 1 ) = dble( 1 )
438  RETURN
439  END IF
440 *
441 * Initialize Q and Z
442 *
443  IF( icompq.EQ.3 )
444  $ CALL dlaset( 'Full', n, n, zero, one, q, ldq )
445  IF( icompz.EQ.3 )
446  $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
447 *
448 * Machine Constants
449 *
450  in = ihi + 1 - ilo
451  safmin = dlamch( 'S' )
452  safmax = one / safmin
453  ulp = dlamch( 'E' )*dlamch( 'B' )
454  anorm = dlanhs( 'F', in, h( ilo, ilo ), ldh, work )
455  bnorm = dlanhs( 'F', in, t( ilo, ilo ), ldt, work )
456  atol = max( safmin, ulp*anorm )
457  btol = max( safmin, ulp*bnorm )
458  ascale = one / max( safmin, anorm )
459  bscale = one / max( safmin, bnorm )
460 *
461 * Set Eigenvalues IHI+1:N
462 *
463  DO 30 j = ihi + 1, n
464  IF( t( j, j ).LT.zero ) THEN
465  IF( ilschr ) THEN
466  DO 10 jr = 1, j
467  h( jr, j ) = -h( jr, j )
468  t( jr, j ) = -t( jr, j )
469  10 CONTINUE
470  ELSE
471  h( j, j ) = -h( j, j )
472  t( j, j ) = -t( j, j )
473  END IF
474  IF( ilz ) THEN
475  DO 20 jr = 1, n
476  z( jr, j ) = -z( jr, j )
477  20 CONTINUE
478  END IF
479  END IF
480  alphar( j ) = h( j, j )
481  alphai( j ) = zero
482  beta( j ) = t( j, j )
483  30 CONTINUE
484 *
485 * If IHI < ILO, skip QZ steps
486 *
487  IF( ihi.LT.ilo )
488  $ GO TO 380
489 *
490 * MAIN QZ ITERATION LOOP
491 *
492 * Initialize dynamic indices
493 *
494 * Eigenvalues ILAST+1:N have been found.
495 * Column operations modify rows IFRSTM:whatever.
496 * Row operations modify columns whatever:ILASTM.
497 *
498 * If only eigenvalues are being computed, then
499 * IFRSTM is the row of the last splitting row above row ILAST;
500 * this is always at least ILO.
501 * IITER counts iterations since the last eigenvalue was found,
502 * to tell when to use an extraordinary shift.
503 * MAXIT is the maximum number of QZ sweeps allowed.
504 *
505  ilast = ihi
506  IF( ilschr ) THEN
507  ifrstm = 1
508  ilastm = n
509  ELSE
510  ifrstm = ilo
511  ilastm = ihi
512  END IF
513  iiter = 0
514  eshift = zero
515  maxit = 30*( ihi-ilo+1 )
516 *
517  DO 360 jiter = 1, maxit
518 *
519 * Split the matrix if possible.
520 *
521 * Two tests:
522 * 1: H(j,j-1)=0 or j=ILO
523 * 2: T(j,j)=0
524 *
525  IF( ilast.EQ.ilo ) THEN
526 *
527 * Special case: j=ILAST
528 *
529  GO TO 80
530  ELSE
531  IF( abs( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*(
532  $ abs( h( ilast, ilast ) ) + abs( h( ilast-1, ilast-1 ) )
533  $ ) ) ) THEN
534  h( ilast, ilast-1 ) = zero
535  GO TO 80
536  END IF
537  END IF
538 *
539  IF( abs( t( ilast, ilast ) ).LE.max( safmin, ulp*(
540  $ abs( t( ilast - 1, ilast ) ) + abs( t( ilast-1, ilast-1 )
541  $ ) ) ) ) THEN
542  t( ilast, ilast ) = zero
543  GO TO 70
544  END IF
545 *
546 * General case: j<ILAST
547 *
548  DO 60 j = ilast - 1, ilo, -1
549 *
550 * Test 1: for H(j,j-1)=0 or j=ILO
551 *
552  IF( j.EQ.ilo ) THEN
553  ilazro = .true.
554  ELSE
555  IF( abs( h( j, j-1 ) ).LE.max( safmin, ulp*(
556  $ abs( h( j, j ) ) + abs( h( j-1, j-1 ) )
557  $ ) ) ) THEN
558  h( j, j-1 ) = zero
559  ilazro = .true.
560  ELSE
561  ilazro = .false.
562  END IF
563  END IF
564 *
565 * Test 2: for T(j,j)=0
566 *
567  temp = abs( t( j, j + 1 ) )
568  IF ( j .GT. ilo )
569  $ temp = temp + abs( t( j - 1, j ) )
570  IF( abs( t( j, j ) ).LT.max( safmin,ulp*temp ) ) THEN
571  t( j, j ) = zero
572 *
573 * Test 1a: Check for 2 consecutive small subdiagonals in A
574 *
575  ilazr2 = .false.
576  IF( .NOT.ilazro ) THEN
577  temp = abs( h( j, j-1 ) )
578  temp2 = abs( h( j, j ) )
579  tempr = max( temp, temp2 )
580  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
581  temp = temp / tempr
582  temp2 = temp2 / tempr
583  END IF
584  IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
585  $ ( ascale*atol ) )ilazr2 = .true.
586  END IF
587 *
588 * If both tests pass (1 & 2), i.e., the leading diagonal
589 * element of B in the block is zero, split a 1x1 block off
590 * at the top. (I.e., at the J-th row/column) The leading
591 * diagonal element of the remainder can also be zero, so
592 * this may have to be done repeatedly.
593 *
594  IF( ilazro .OR. ilazr2 ) THEN
595  DO 40 jch = j, ilast - 1
596  temp = h( jch, jch )
597  CALL dlartg( temp, h( jch+1, jch ), c, s,
598  $ h( jch, jch ) )
599  h( jch+1, jch ) = zero
600  CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
601  $ h( jch+1, jch+1 ), ldh, c, s )
602  CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
603  $ t( jch+1, jch+1 ), ldt, c, s )
604  IF( ilq )
605  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
606  $ c, s )
607  IF( ilazr2 )
608  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
609  ilazr2 = .false.
610  IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
611  IF( jch+1.GE.ilast ) THEN
612  GO TO 80
613  ELSE
614  ifirst = jch + 1
615  GO TO 110
616  END IF
617  END IF
618  t( jch+1, jch+1 ) = zero
619  40 CONTINUE
620  GO TO 70
621  ELSE
622 *
623 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
624 * Then process as in the case T(ILAST,ILAST)=0
625 *
626  DO 50 jch = j, ilast - 1
627  temp = t( jch, jch+1 )
628  CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
629  $ t( jch, jch+1 ) )
630  t( jch+1, jch+1 ) = zero
631  IF( jch.LT.ilastm-1 )
632  $ CALL drot( ilastm-jch-1, t( jch, jch+2 ), ldt,
633  $ t( jch+1, jch+2 ), ldt, c, s )
634  CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
635  $ h( jch+1, jch-1 ), ldh, c, s )
636  IF( ilq )
637  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
638  $ c, s )
639  temp = h( jch+1, jch )
640  CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
641  $ h( jch+1, jch ) )
642  h( jch+1, jch-1 ) = zero
643  CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
644  $ h( ifrstm, jch-1 ), 1, c, s )
645  CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
646  $ t( ifrstm, jch-1 ), 1, c, s )
647  IF( ilz )
648  $ CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
649  $ c, s )
650  50 CONTINUE
651  GO TO 70
652  END IF
653  ELSE IF( ilazro ) THEN
654 *
655 * Only test 1 passed -- work on J:ILAST
656 *
657  ifirst = j
658  GO TO 110
659  END IF
660 *
661 * Neither test passed -- try next J
662 *
663  60 CONTINUE
664 *
665 * (Drop-through is "impossible")
666 *
667  info = n + 1
668  GO TO 420
669 *
670 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
671 * 1x1 block.
672 *
673  70 CONTINUE
674  temp = h( ilast, ilast )
675  CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
676  $ h( ilast, ilast ) )
677  h( ilast, ilast-1 ) = zero
678  CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
679  $ h( ifrstm, ilast-1 ), 1, c, s )
680  CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
681  $ t( ifrstm, ilast-1 ), 1, c, s )
682  IF( ilz )
683  $ CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
684 *
685 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
686 * and BETA
687 *
688  80 CONTINUE
689  IF( t( ilast, ilast ).LT.zero ) THEN
690  IF( ilschr ) THEN
691  DO 90 j = ifrstm, ilast
692  h( j, ilast ) = -h( j, ilast )
693  t( j, ilast ) = -t( j, ilast )
694  90 CONTINUE
695  ELSE
696  h( ilast, ilast ) = -h( ilast, ilast )
697  t( ilast, ilast ) = -t( ilast, ilast )
698  END IF
699  IF( ilz ) THEN
700  DO 100 j = 1, n
701  z( j, ilast ) = -z( j, ilast )
702  100 CONTINUE
703  END IF
704  END IF
705  alphar( ilast ) = h( ilast, ilast )
706  alphai( ilast ) = zero
707  beta( ilast ) = t( ilast, ilast )
708 *
709 * Go to next block -- exit if finished.
710 *
711  ilast = ilast - 1
712  IF( ilast.LT.ilo )
713  $ GO TO 380
714 *
715 * Reset counters
716 *
717  iiter = 0
718  eshift = zero
719  IF( .NOT.ilschr ) THEN
720  ilastm = ilast
721  IF( ifrstm.GT.ilast )
722  $ ifrstm = ilo
723  END IF
724  GO TO 350
725 *
726 * QZ step
727 *
728 * This iteration only involves rows/columns IFIRST:ILAST. We
729 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
730 *
731  110 CONTINUE
732  iiter = iiter + 1
733  IF( .NOT.ilschr ) THEN
734  ifrstm = ifirst
735  END IF
736 *
737 * Compute single shifts.
738 *
739 * At this point, IFIRST < ILAST, and the diagonal elements of
740 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
741 * magnitude)
742 *
743  IF( ( iiter / 10 )*10.EQ.iiter ) THEN
744 *
745 * Exceptional shift. Chosen for no particularly good reason.
746 * (Single shift only.)
747 *
748  IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
749  $ abs( t( ilast-1, ilast-1 ) ) ) THEN
750  eshift = h( ilast, ilast-1 ) /
751  $ t( ilast-1, ilast-1 )
752  ELSE
753  eshift = eshift + one / ( safmin*dble( maxit ) )
754  END IF
755  s1 = one
756  wr = eshift
757 *
758  ELSE
759 *
760 * Shifts based on the generalized eigenvalues of the
761 * bottom-right 2x2 block of A and B. The first eigenvalue
762 * returned by DLAG2 is the Wilkinson shift (AEP p.512),
763 *
764  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
765  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
766  $ s2, wr, wr2, wi )
767 *
768  IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
769  $ .GT. abs( (wr2/s2)*t( ilast, ilast )
770  $ - h( ilast, ilast ) ) ) THEN
771  temp = wr
772  wr = wr2
773  wr2 = temp
774  temp = s1
775  s1 = s2
776  s2 = temp
777  END IF
778  temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
779  IF( wi.NE.zero )
780  $ GO TO 200
781  END IF
782 *
783 * Fiddle with shift to avoid overflow
784 *
785  temp = min( ascale, one )*( half*safmax )
786  IF( s1.GT.temp ) THEN
787  scale = temp / s1
788  ELSE
789  scale = one
790  END IF
791 *
792  temp = min( bscale, one )*( half*safmax )
793  IF( abs( wr ).GT.temp )
794  $ scale = min( scale, temp / abs( wr ) )
795  s1 = scale*s1
796  wr = scale*wr
797 *
798 * Now check for two consecutive small subdiagonals.
799 *
800  DO 120 j = ilast - 1, ifirst + 1, -1
801  istart = j
802  temp = abs( s1*h( j, j-1 ) )
803  temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
804  tempr = max( temp, temp2 )
805  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
806  temp = temp / tempr
807  temp2 = temp2 / tempr
808  END IF
809  IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
810  $ temp2 )GO TO 130
811  120 CONTINUE
812 *
813  istart = ifirst
814  130 CONTINUE
815 *
816 * Do an implicit single-shift QZ sweep.
817 *
818 * Initial Q
819 *
820  temp = s1*h( istart, istart ) - wr*t( istart, istart )
821  temp2 = s1*h( istart+1, istart )
822  CALL dlartg( temp, temp2, c, s, tempr )
823 *
824 * Sweep
825 *
826  DO 190 j = istart, ilast - 1
827  IF( j.GT.istart ) THEN
828  temp = h( j, j-1 )
829  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
830  h( j+1, j-1 ) = zero
831  END IF
832 *
833  DO 140 jc = j, ilastm
834  temp = c*h( j, jc ) + s*h( j+1, jc )
835  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
836  h( j, jc ) = temp
837  temp2 = c*t( j, jc ) + s*t( j+1, jc )
838  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
839  t( j, jc ) = temp2
840  140 CONTINUE
841  IF( ilq ) THEN
842  DO 150 jr = 1, n
843  temp = c*q( jr, j ) + s*q( jr, j+1 )
844  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
845  q( jr, j ) = temp
846  150 CONTINUE
847  END IF
848 *
849  temp = t( j+1, j+1 )
850  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
851  t( j+1, j ) = zero
852 *
853  DO 160 jr = ifrstm, min( j+2, ilast )
854  temp = c*h( jr, j+1 ) + s*h( jr, j )
855  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
856  h( jr, j+1 ) = temp
857  160 CONTINUE
858  DO 170 jr = ifrstm, j
859  temp = c*t( jr, j+1 ) + s*t( jr, j )
860  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
861  t( jr, j+1 ) = temp
862  170 CONTINUE
863  IF( ilz ) THEN
864  DO 180 jr = 1, n
865  temp = c*z( jr, j+1 ) + s*z( jr, j )
866  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
867  z( jr, j+1 ) = temp
868  180 CONTINUE
869  END IF
870  190 CONTINUE
871 *
872  GO TO 350
873 *
874 * Use Francis double-shift
875 *
876 * Note: the Francis double-shift should work with real shifts,
877 * but only if the block is at least 3x3.
878 * This code may break if this point is reached with
879 * a 2x2 block with real eigenvalues.
880 *
881  200 CONTINUE
882  IF( ifirst+1.EQ.ilast ) THEN
883 *
884 * Special case -- 2x2 block with complex eigenvectors
885 *
886 * Step 1: Standardize, that is, rotate so that
887 *
888 * ( B11 0 )
889 * B = ( ) with B11 non-negative.
890 * ( 0 B22 )
891 *
892  CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
893  $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
894 *
895  IF( b11.LT.zero ) THEN
896  cr = -cr
897  sr = -sr
898  b11 = -b11
899  b22 = -b22
900  END IF
901 *
902  CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
903  $ h( ilast, ilast-1 ), ldh, cl, sl )
904  CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
905  $ h( ifrstm, ilast ), 1, cr, sr )
906 *
907  IF( ilast.LT.ilastm )
908  $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
909  $ t( ilast, ilast+1 ), ldt, cl, sl )
910  IF( ifrstm.LT.ilast-1 )
911  $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
912  $ t( ifrstm, ilast ), 1, cr, sr )
913 *
914  IF( ilq )
915  $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
916  $ sl )
917  IF( ilz )
918  $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
919  $ sr )
920 *
921  t( ilast-1, ilast-1 ) = b11
922  t( ilast-1, ilast ) = zero
923  t( ilast, ilast-1 ) = zero
924  t( ilast, ilast ) = b22
925 *
926 * If B22 is negative, negate column ILAST
927 *
928  IF( b22.LT.zero ) THEN
929  DO 210 j = ifrstm, ilast
930  h( j, ilast ) = -h( j, ilast )
931  t( j, ilast ) = -t( j, ilast )
932  210 CONTINUE
933 *
934  IF( ilz ) THEN
935  DO 220 j = 1, n
936  z( j, ilast ) = -z( j, ilast )
937  220 CONTINUE
938  END IF
939  b22 = -b22
940  END IF
941 *
942 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
943 *
944 * Recompute shift
945 *
946  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
947  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
948  $ temp, wr, temp2, wi )
949 *
950 * If standardization has perturbed the shift onto real line,
951 * do another (real single-shift) QR step.
952 *
953  IF( wi.EQ.zero )
954  $ GO TO 350
955  s1inv = one / s1
956 *
957 * Do EISPACK (QZVAL) computation of alpha and beta
958 *
959  a11 = h( ilast-1, ilast-1 )
960  a21 = h( ilast, ilast-1 )
961  a12 = h( ilast-1, ilast )
962  a22 = h( ilast, ilast )
963 *
964 * Compute complex Givens rotation on right
965 * (Assume some element of C = (sA - wB) > unfl )
966 * __
967 * (sA - wB) ( CZ -SZ )
968 * ( SZ CZ )
969 *
970  c11r = s1*a11 - wr*b11
971  c11i = -wi*b11
972  c12 = s1*a12
973  c21 = s1*a21
974  c22r = s1*a22 - wr*b22
975  c22i = -wi*b22
976 *
977  IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
978  $ abs( c22r )+abs( c22i ) ) THEN
979  t1 = dlapy3( c12, c11r, c11i )
980  cz = c12 / t1
981  szr = -c11r / t1
982  szi = -c11i / t1
983  ELSE
984  cz = dlapy2( c22r, c22i )
985  IF( cz.LE.safmin ) THEN
986  cz = zero
987  szr = one
988  szi = zero
989  ELSE
990  tempr = c22r / cz
991  tempi = c22i / cz
992  t1 = dlapy2( cz, c21 )
993  cz = cz / t1
994  szr = -c21*tempr / t1
995  szi = c21*tempi / t1
996  END IF
997  END IF
998 *
999 * Compute Givens rotation on left
1000 *
1001 * ( CQ SQ )
1002 * ( __ ) A or B
1003 * ( -SQ CQ )
1004 *
1005  an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1006  bn = abs( b11 ) + abs( b22 )
1007  wabs = abs( wr ) + abs( wi )
1008  IF( s1*an.GT.wabs*bn ) THEN
1009  cq = cz*b11
1010  sqr = szr*b22
1011  sqi = -szi*b22
1012  ELSE
1013  a1r = cz*a11 + szr*a12
1014  a1i = szi*a12
1015  a2r = cz*a21 + szr*a22
1016  a2i = szi*a22
1017  cq = dlapy2( a1r, a1i )
1018  IF( cq.LE.safmin ) THEN
1019  cq = zero
1020  sqr = one
1021  sqi = zero
1022  ELSE
1023  tempr = a1r / cq
1024  tempi = a1i / cq
1025  sqr = tempr*a2r + tempi*a2i
1026  sqi = tempi*a2r - tempr*a2i
1027  END IF
1028  END IF
1029  t1 = dlapy3( cq, sqr, sqi )
1030  cq = cq / t1
1031  sqr = sqr / t1
1032  sqi = sqi / t1
1033 *
1034 * Compute diagonal elements of QBZ
1035 *
1036  tempr = sqr*szr - sqi*szi
1037  tempi = sqr*szi + sqi*szr
1038  b1r = cq*cz*b11 + tempr*b22
1039  b1i = tempi*b22
1040  b1a = dlapy2( b1r, b1i )
1041  b2r = cq*cz*b22 + tempr*b11
1042  b2i = -tempi*b11
1043  b2a = dlapy2( b2r, b2i )
1044 *
1045 * Normalize so beta > 0, and Im( alpha1 ) > 0
1046 *
1047  beta( ilast-1 ) = b1a
1048  beta( ilast ) = b2a
1049  alphar( ilast-1 ) = ( wr*b1a )*s1inv
1050  alphai( ilast-1 ) = ( wi*b1a )*s1inv
1051  alphar( ilast ) = ( wr*b2a )*s1inv
1052  alphai( ilast ) = -( wi*b2a )*s1inv
1053 *
1054 * Step 3: Go to next block -- exit if finished.
1055 *
1056  ilast = ifirst - 1
1057  IF( ilast.LT.ilo )
1058  $ GO TO 380
1059 *
1060 * Reset counters
1061 *
1062  iiter = 0
1063  eshift = zero
1064  IF( .NOT.ilschr ) THEN
1065  ilastm = ilast
1066  IF( ifrstm.GT.ilast )
1067  $ ifrstm = ilo
1068  END IF
1069  GO TO 350
1070  ELSE
1071 *
1072 * Usual case: 3x3 or larger block, using Francis implicit
1073 * double-shift
1074 *
1075 * 2
1076 * Eigenvalue equation is w - c w + d = 0,
1077 *
1078 * -1 2 -1
1079 * so compute 1st column of (A B ) - c A B + d
1080 * using the formula in QZIT (from EISPACK)
1081 *
1082 * We assume that the block is at least 3x3
1083 *
1084  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1085  $ ( bscale*t( ilast-1, ilast-1 ) )
1086  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1087  $ ( bscale*t( ilast-1, ilast-1 ) )
1088  ad12 = ( ascale*h( ilast-1, ilast ) ) /
1089  $ ( bscale*t( ilast, ilast ) )
1090  ad22 = ( ascale*h( ilast, ilast ) ) /
1091  $ ( bscale*t( ilast, ilast ) )
1092  u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1093  ad11l = ( ascale*h( ifirst, ifirst ) ) /
1094  $ ( bscale*t( ifirst, ifirst ) )
1095  ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1096  $ ( bscale*t( ifirst, ifirst ) )
1097  ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1098  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1099  ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1100  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1101  ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1102  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1103  u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1104 *
1105  v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1106  $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1107  v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1108  $ ( ad22-ad11l )+ad21*u12 )*ad21l
1109  v( 3 ) = ad32l*ad21l
1110 *
1111  istart = ifirst
1112 *
1113  CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1114  v( 1 ) = one
1115 *
1116 * Sweep
1117 *
1118  DO 290 j = istart, ilast - 2
1119 *
1120 * All but last elements: use 3x3 Householder transforms.
1121 *
1122 * Zero (j-1)st column of A
1123 *
1124  IF( j.GT.istart ) THEN
1125  v( 1 ) = h( j, j-1 )
1126  v( 2 ) = h( j+1, j-1 )
1127  v( 3 ) = h( j+2, j-1 )
1128 *
1129  CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1130  v( 1 ) = one
1131  h( j+1, j-1 ) = zero
1132  h( j+2, j-1 ) = zero
1133  END IF
1134 *
1135  DO 230 jc = j, ilastm
1136  temp = tau*( h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1137  $ h( j+2, jc ) )
1138  h( j, jc ) = h( j, jc ) - temp
1139  h( j+1, jc ) = h( j+1, jc ) - temp*v( 2 )
1140  h( j+2, jc ) = h( j+2, jc ) - temp*v( 3 )
1141  temp2 = tau*( t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1142  $ t( j+2, jc ) )
1143  t( j, jc ) = t( j, jc ) - temp2
1144  t( j+1, jc ) = t( j+1, jc ) - temp2*v( 2 )
1145  t( j+2, jc ) = t( j+2, jc ) - temp2*v( 3 )
1146  230 CONTINUE
1147  IF( ilq ) THEN
1148  DO 240 jr = 1, n
1149  temp = tau*( q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1150  $ q( jr, j+2 ) )
1151  q( jr, j ) = q( jr, j ) - temp
1152  q( jr, j+1 ) = q( jr, j+1 ) - temp*v( 2 )
1153  q( jr, j+2 ) = q( jr, j+2 ) - temp*v( 3 )
1154  240 CONTINUE
1155  END IF
1156 *
1157 * Zero j-th column of B (see DLAGBC for details)
1158 *
1159 * Swap rows to pivot
1160 *
1161  ilpivt = .false.
1162  temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1163  temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1164  IF( max( temp, temp2 ).LT.safmin ) THEN
1165  scale = zero
1166  u1 = one
1167  u2 = zero
1168  GO TO 250
1169  ELSE IF( temp.GE.temp2 ) THEN
1170  w11 = t( j+1, j+1 )
1171  w21 = t( j+2, j+1 )
1172  w12 = t( j+1, j+2 )
1173  w22 = t( j+2, j+2 )
1174  u1 = t( j+1, j )
1175  u2 = t( j+2, j )
1176  ELSE
1177  w21 = t( j+1, j+1 )
1178  w11 = t( j+2, j+1 )
1179  w22 = t( j+1, j+2 )
1180  w12 = t( j+2, j+2 )
1181  u2 = t( j+1, j )
1182  u1 = t( j+2, j )
1183  END IF
1184 *
1185 * Swap columns if nec.
1186 *
1187  IF( abs( w12 ).GT.abs( w11 ) ) THEN
1188  ilpivt = .true.
1189  temp = w12
1190  temp2 = w22
1191  w12 = w11
1192  w22 = w21
1193  w11 = temp
1194  w21 = temp2
1195  END IF
1196 *
1197 * LU-factor
1198 *
1199  temp = w21 / w11
1200  u2 = u2 - temp*u1
1201  w22 = w22 - temp*w12
1202  w21 = zero
1203 *
1204 * Compute SCALE
1205 *
1206  scale = one
1207  IF( abs( w22 ).LT.safmin ) THEN
1208  scale = zero
1209  u2 = one
1210  u1 = -w12 / w11
1211  GO TO 250
1212  END IF
1213  IF( abs( w22 ).LT.abs( u2 ) )
1214  $ scale = abs( w22 / u2 )
1215  IF( abs( w11 ).LT.abs( u1 ) )
1216  $ scale = min( scale, abs( w11 / u1 ) )
1217 *
1218 * Solve
1219 *
1220  u2 = ( scale*u2 ) / w22
1221  u1 = ( scale*u1-w12*u2 ) / w11
1222 *
1223  250 CONTINUE
1224  IF( ilpivt ) THEN
1225  temp = u2
1226  u2 = u1
1227  u1 = temp
1228  END IF
1229 *
1230 * Compute Householder Vector
1231 *
1232  t1 = sqrt( scale**2+u1**2+u2**2 )
1233  tau = one + scale / t1
1234  vs = -one / ( scale+t1 )
1235  v( 1 ) = one
1236  v( 2 ) = vs*u1
1237  v( 3 ) = vs*u2
1238 *
1239 * Apply transformations from the right.
1240 *
1241  DO 260 jr = ifrstm, min( j+3, ilast )
1242  temp = tau*( h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1243  $ h( jr, j+2 ) )
1244  h( jr, j ) = h( jr, j ) - temp
1245  h( jr, j+1 ) = h( jr, j+1 ) - temp*v( 2 )
1246  h( jr, j+2 ) = h( jr, j+2 ) - temp*v( 3 )
1247  260 CONTINUE
1248  DO 270 jr = ifrstm, j + 2
1249  temp = tau*( t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1250  $ t( jr, j+2 ) )
1251  t( jr, j ) = t( jr, j ) - temp
1252  t( jr, j+1 ) = t( jr, j+1 ) - temp*v( 2 )
1253  t( jr, j+2 ) = t( jr, j+2 ) - temp*v( 3 )
1254  270 CONTINUE
1255  IF( ilz ) THEN
1256  DO 280 jr = 1, n
1257  temp = tau*( z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1258  $ z( jr, j+2 ) )
1259  z( jr, j ) = z( jr, j ) - temp
1260  z( jr, j+1 ) = z( jr, j+1 ) - temp*v( 2 )
1261  z( jr, j+2 ) = z( jr, j+2 ) - temp*v( 3 )
1262  280 CONTINUE
1263  END IF
1264  t( j+1, j ) = zero
1265  t( j+2, j ) = zero
1266  290 CONTINUE
1267 *
1268 * Last elements: Use Givens rotations
1269 *
1270 * Rotations from the left
1271 *
1272  j = ilast - 1
1273  temp = h( j, j-1 )
1274  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1275  h( j+1, j-1 ) = zero
1276 *
1277  DO 300 jc = j, ilastm
1278  temp = c*h( j, jc ) + s*h( j+1, jc )
1279  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1280  h( j, jc ) = temp
1281  temp2 = c*t( j, jc ) + s*t( j+1, jc )
1282  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1283  t( j, jc ) = temp2
1284  300 CONTINUE
1285  IF( ilq ) THEN
1286  DO 310 jr = 1, n
1287  temp = c*q( jr, j ) + s*q( jr, j+1 )
1288  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1289  q( jr, j ) = temp
1290  310 CONTINUE
1291  END IF
1292 *
1293 * Rotations from the right.
1294 *
1295  temp = t( j+1, j+1 )
1296  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1297  t( j+1, j ) = zero
1298 *
1299  DO 320 jr = ifrstm, ilast
1300  temp = c*h( jr, j+1 ) + s*h( jr, j )
1301  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1302  h( jr, j+1 ) = temp
1303  320 CONTINUE
1304  DO 330 jr = ifrstm, ilast - 1
1305  temp = c*t( jr, j+1 ) + s*t( jr, j )
1306  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1307  t( jr, j+1 ) = temp
1308  330 CONTINUE
1309  IF( ilz ) THEN
1310  DO 340 jr = 1, n
1311  temp = c*z( jr, j+1 ) + s*z( jr, j )
1312  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1313  z( jr, j+1 ) = temp
1314  340 CONTINUE
1315  END IF
1316 *
1317 * End of Double-Shift code
1318 *
1319  END IF
1320 *
1321  GO TO 350
1322 *
1323 * End of iteration loop
1324 *
1325  350 CONTINUE
1326  360 CONTINUE
1327 *
1328 * Drop-through = non-convergence
1329 *
1330  info = ilast
1331  GO TO 420
1332 *
1333 * Successful completion of all QZ steps
1334 *
1335  380 CONTINUE
1336 *
1337 * Set Eigenvalues 1:ILO-1
1338 *
1339  DO 410 j = 1, ilo - 1
1340  IF( t( j, j ).LT.zero ) THEN
1341  IF( ilschr ) THEN
1342  DO 390 jr = 1, j
1343  h( jr, j ) = -h( jr, j )
1344  t( jr, j ) = -t( jr, j )
1345  390 CONTINUE
1346  ELSE
1347  h( j, j ) = -h( j, j )
1348  t( j, j ) = -t( j, j )
1349  END IF
1350  IF( ilz ) THEN
1351  DO 400 jr = 1, n
1352  z( jr, j ) = -z( jr, j )
1353  400 CONTINUE
1354  END IF
1355  END IF
1356  alphar( j ) = h( j, j )
1357  alphai( j ) = zero
1358  beta( j ) = t( j, j )
1359  410 CONTINUE
1360 *
1361 * Normal Termination
1362 *
1363  info = 0
1364 *
1365 * Exit (other than argument error) -- return optimal workspace size
1366 *
1367  420 CONTINUE
1368  work( 1 ) = dble( n )
1369  RETURN
1370 *
1371 * End of DHGEQZ
1372 *
1373  END
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine dlasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: dlasv2.f:138
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
Definition: dhgeqz.f:304
subroutine dlag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: dlag2.f:156
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106