LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dhgeqz.f
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 hgeqz
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.btol ) THEN
540  t( ilast, ilast ) = zero
541  GO TO 70
542  END IF
543 *
544 * General case: j<ILAST
545 *
546  DO 60 j = ilast - 1, ilo, -1
547 *
548 * Test 1: for H(j,j-1)=0 or j=ILO
549 *
550  IF( j.EQ.ilo ) THEN
551  ilazro = .true.
552  ELSE
553  IF( abs( h( j, j-1 ) ).LE.max( safmin, ulp*(
554  $ abs( h( j, j ) ) + abs( h( j-1, j-1 ) )
555  $ ) ) ) THEN
556  h( j, j-1 ) = zero
557  ilazro = .true.
558  ELSE
559  ilazro = .false.
560  END IF
561  END IF
562 *
563 * Test 2: for T(j,j)=0
564 *
565  IF( abs( t( j, j ) ).LT.btol ) THEN
566  t( j, j ) = zero
567 *
568 * Test 1a: Check for 2 consecutive small subdiagonals in A
569 *
570  ilazr2 = .false.
571  IF( .NOT.ilazro ) THEN
572  temp = abs( h( j, j-1 ) )
573  temp2 = abs( h( j, j ) )
574  tempr = max( temp, temp2 )
575  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
576  temp = temp / tempr
577  temp2 = temp2 / tempr
578  END IF
579  IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
580  $ ( ascale*atol ) )ilazr2 = .true.
581  END IF
582 *
583 * If both tests pass (1 & 2), i.e., the leading diagonal
584 * element of B in the block is zero, split a 1x1 block off
585 * at the top. (I.e., at the J-th row/column) The leading
586 * diagonal element of the remainder can also be zero, so
587 * this may have to be done repeatedly.
588 *
589  IF( ilazro .OR. ilazr2 ) THEN
590  DO 40 jch = j, ilast - 1
591  temp = h( jch, jch )
592  CALL dlartg( temp, h( jch+1, jch ), c, s,
593  $ h( jch, jch ) )
594  h( jch+1, jch ) = zero
595  CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
596  $ h( jch+1, jch+1 ), ldh, c, s )
597  CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
598  $ t( jch+1, jch+1 ), ldt, c, s )
599  IF( ilq )
600  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
601  $ c, s )
602  IF( ilazr2 )
603  $ h( jch, jch-1 ) = h( jch, jch-1 )*c
604  ilazr2 = .false.
605  IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
606  IF( jch+1.GE.ilast ) THEN
607  GO TO 80
608  ELSE
609  ifirst = jch + 1
610  GO TO 110
611  END IF
612  END IF
613  t( jch+1, jch+1 ) = zero
614  40 CONTINUE
615  GO TO 70
616  ELSE
617 *
618 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
619 * Then process as in the case T(ILAST,ILAST)=0
620 *
621  DO 50 jch = j, ilast - 1
622  temp = t( jch, jch+1 )
623  CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
624  $ t( jch, jch+1 ) )
625  t( jch+1, jch+1 ) = zero
626  IF( jch.LT.ilastm-1 )
627  $ CALL drot( ilastm-jch-1, t( jch, jch+2 ), ldt,
628  $ t( jch+1, jch+2 ), ldt, c, s )
629  CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
630  $ h( jch+1, jch-1 ), ldh, c, s )
631  IF( ilq )
632  $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
633  $ c, s )
634  temp = h( jch+1, jch )
635  CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
636  $ h( jch+1, jch ) )
637  h( jch+1, jch-1 ) = zero
638  CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
639  $ h( ifrstm, jch-1 ), 1, c, s )
640  CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
641  $ t( ifrstm, jch-1 ), 1, c, s )
642  IF( ilz )
643  $ CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
644  $ c, s )
645  50 CONTINUE
646  GO TO 70
647  END IF
648  ELSE IF( ilazro ) THEN
649 *
650 * Only test 1 passed -- work on J:ILAST
651 *
652  ifirst = j
653  GO TO 110
654  END IF
655 *
656 * Neither test passed -- try next J
657 *
658  60 CONTINUE
659 *
660 * (Drop-through is "impossible")
661 *
662  info = n + 1
663  GO TO 420
664 *
665 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
666 * 1x1 block.
667 *
668  70 CONTINUE
669  temp = h( ilast, ilast )
670  CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
671  $ h( ilast, ilast ) )
672  h( ilast, ilast-1 ) = zero
673  CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
674  $ h( ifrstm, ilast-1 ), 1, c, s )
675  CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
676  $ t( ifrstm, ilast-1 ), 1, c, s )
677  IF( ilz )
678  $ CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
679 *
680 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
681 * and BETA
682 *
683  80 CONTINUE
684  IF( t( ilast, ilast ).LT.zero ) THEN
685  IF( ilschr ) THEN
686  DO 90 j = ifrstm, ilast
687  h( j, ilast ) = -h( j, ilast )
688  t( j, ilast ) = -t( j, ilast )
689  90 CONTINUE
690  ELSE
691  h( ilast, ilast ) = -h( ilast, ilast )
692  t( ilast, ilast ) = -t( ilast, ilast )
693  END IF
694  IF( ilz ) THEN
695  DO 100 j = 1, n
696  z( j, ilast ) = -z( j, ilast )
697  100 CONTINUE
698  END IF
699  END IF
700  alphar( ilast ) = h( ilast, ilast )
701  alphai( ilast ) = zero
702  beta( ilast ) = t( ilast, ilast )
703 *
704 * Go to next block -- exit if finished.
705 *
706  ilast = ilast - 1
707  IF( ilast.LT.ilo )
708  $ GO TO 380
709 *
710 * Reset counters
711 *
712  iiter = 0
713  eshift = zero
714  IF( .NOT.ilschr ) THEN
715  ilastm = ilast
716  IF( ifrstm.GT.ilast )
717  $ ifrstm = ilo
718  END IF
719  GO TO 350
720 *
721 * QZ step
722 *
723 * This iteration only involves rows/columns IFIRST:ILAST. We
724 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
725 *
726  110 CONTINUE
727  iiter = iiter + 1
728  IF( .NOT.ilschr ) THEN
729  ifrstm = ifirst
730  END IF
731 *
732 * Compute single shifts.
733 *
734 * At this point, IFIRST < ILAST, and the diagonal elements of
735 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
736 * magnitude)
737 *
738  IF( ( iiter / 10 )*10.EQ.iiter ) THEN
739 *
740 * Exceptional shift. Chosen for no particularly good reason.
741 * (Single shift only.)
742 *
743  IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
744  $ abs( t( ilast-1, ilast-1 ) ) ) THEN
745  eshift = h( ilast, ilast-1 ) /
746  $ t( ilast-1, ilast-1 )
747  ELSE
748  eshift = eshift + one / ( safmin*dble( maxit ) )
749  END IF
750  s1 = one
751  wr = eshift
752 *
753  ELSE
754 *
755 * Shifts based on the generalized eigenvalues of the
756 * bottom-right 2x2 block of A and B. The first eigenvalue
757 * returned by DLAG2 is the Wilkinson shift (AEP p.512),
758 *
759  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
760  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
761  $ s2, wr, wr2, wi )
762 *
763  IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
764  $ .GT. abs( (wr2/s2)*t( ilast, ilast )
765  $ - h( ilast, ilast ) ) ) THEN
766  temp = wr
767  wr = wr2
768  wr2 = temp
769  temp = s1
770  s1 = s2
771  s2 = temp
772  END IF
773  temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
774  IF( wi.NE.zero )
775  $ GO TO 200
776  END IF
777 *
778 * Fiddle with shift to avoid overflow
779 *
780  temp = min( ascale, one )*( half*safmax )
781  IF( s1.GT.temp ) THEN
782  scale = temp / s1
783  ELSE
784  scale = one
785  END IF
786 *
787  temp = min( bscale, one )*( half*safmax )
788  IF( abs( wr ).GT.temp )
789  $ scale = min( scale, temp / abs( wr ) )
790  s1 = scale*s1
791  wr = scale*wr
792 *
793 * Now check for two consecutive small subdiagonals.
794 *
795  DO 120 j = ilast - 1, ifirst + 1, -1
796  istart = j
797  temp = abs( s1*h( j, j-1 ) )
798  temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
799  tempr = max( temp, temp2 )
800  IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
801  temp = temp / tempr
802  temp2 = temp2 / tempr
803  END IF
804  IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
805  $ temp2 )GO TO 130
806  120 CONTINUE
807 *
808  istart = ifirst
809  130 CONTINUE
810 *
811 * Do an implicit single-shift QZ sweep.
812 *
813 * Initial Q
814 *
815  temp = s1*h( istart, istart ) - wr*t( istart, istart )
816  temp2 = s1*h( istart+1, istart )
817  CALL dlartg( temp, temp2, c, s, tempr )
818 *
819 * Sweep
820 *
821  DO 190 j = istart, ilast - 1
822  IF( j.GT.istart ) THEN
823  temp = h( j, j-1 )
824  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
825  h( j+1, j-1 ) = zero
826  END IF
827 *
828  DO 140 jc = j, ilastm
829  temp = c*h( j, jc ) + s*h( j+1, jc )
830  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
831  h( j, jc ) = temp
832  temp2 = c*t( j, jc ) + s*t( j+1, jc )
833  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
834  t( j, jc ) = temp2
835  140 CONTINUE
836  IF( ilq ) THEN
837  DO 150 jr = 1, n
838  temp = c*q( jr, j ) + s*q( jr, j+1 )
839  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
840  q( jr, j ) = temp
841  150 CONTINUE
842  END IF
843 *
844  temp = t( j+1, j+1 )
845  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
846  t( j+1, j ) = zero
847 *
848  DO 160 jr = ifrstm, min( j+2, ilast )
849  temp = c*h( jr, j+1 ) + s*h( jr, j )
850  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
851  h( jr, j+1 ) = temp
852  160 CONTINUE
853  DO 170 jr = ifrstm, j
854  temp = c*t( jr, j+1 ) + s*t( jr, j )
855  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
856  t( jr, j+1 ) = temp
857  170 CONTINUE
858  IF( ilz ) THEN
859  DO 180 jr = 1, n
860  temp = c*z( jr, j+1 ) + s*z( jr, j )
861  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
862  z( jr, j+1 ) = temp
863  180 CONTINUE
864  END IF
865  190 CONTINUE
866 *
867  GO TO 350
868 *
869 * Use Francis double-shift
870 *
871 * Note: the Francis double-shift should work with real shifts,
872 * but only if the block is at least 3x3.
873 * This code may break if this point is reached with
874 * a 2x2 block with real eigenvalues.
875 *
876  200 CONTINUE
877  IF( ifirst+1.EQ.ilast ) THEN
878 *
879 * Special case -- 2x2 block with complex eigenvectors
880 *
881 * Step 1: Standardize, that is, rotate so that
882 *
883 * ( B11 0 )
884 * B = ( ) with B11 non-negative.
885 * ( 0 B22 )
886 *
887  CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
888  $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
889 *
890  IF( b11.LT.zero ) THEN
891  cr = -cr
892  sr = -sr
893  b11 = -b11
894  b22 = -b22
895  END IF
896 *
897  CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
898  $ h( ilast, ilast-1 ), ldh, cl, sl )
899  CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
900  $ h( ifrstm, ilast ), 1, cr, sr )
901 *
902  IF( ilast.LT.ilastm )
903  $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
904  $ t( ilast, ilast+1 ), ldt, cl, sl )
905  IF( ifrstm.LT.ilast-1 )
906  $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
907  $ t( ifrstm, ilast ), 1, cr, sr )
908 *
909  IF( ilq )
910  $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
911  $ sl )
912  IF( ilz )
913  $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
914  $ sr )
915 *
916  t( ilast-1, ilast-1 ) = b11
917  t( ilast-1, ilast ) = zero
918  t( ilast, ilast-1 ) = zero
919  t( ilast, ilast ) = b22
920 *
921 * If B22 is negative, negate column ILAST
922 *
923  IF( b22.LT.zero ) THEN
924  DO 210 j = ifrstm, ilast
925  h( j, ilast ) = -h( j, ilast )
926  t( j, ilast ) = -t( j, ilast )
927  210 CONTINUE
928 *
929  IF( ilz ) THEN
930  DO 220 j = 1, n
931  z( j, ilast ) = -z( j, ilast )
932  220 CONTINUE
933  END IF
934  b22 = -b22
935  END IF
936 *
937 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
938 *
939 * Recompute shift
940 *
941  CALL dlag2( h( ilast-1, ilast-1 ), ldh,
942  $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
943  $ temp, wr, temp2, wi )
944 *
945 * If standardization has perturbed the shift onto real line,
946 * do another (real single-shift) QR step.
947 *
948  IF( wi.EQ.zero )
949  $ GO TO 350
950  s1inv = one / s1
951 *
952 * Do EISPACK (QZVAL) computation of alpha and beta
953 *
954  a11 = h( ilast-1, ilast-1 )
955  a21 = h( ilast, ilast-1 )
956  a12 = h( ilast-1, ilast )
957  a22 = h( ilast, ilast )
958 *
959 * Compute complex Givens rotation on right
960 * (Assume some element of C = (sA - wB) > unfl )
961 * __
962 * (sA - wB) ( CZ -SZ )
963 * ( SZ CZ )
964 *
965  c11r = s1*a11 - wr*b11
966  c11i = -wi*b11
967  c12 = s1*a12
968  c21 = s1*a21
969  c22r = s1*a22 - wr*b22
970  c22i = -wi*b22
971 *
972  IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
973  $ abs( c22r )+abs( c22i ) ) THEN
974  t1 = dlapy3( c12, c11r, c11i )
975  cz = c12 / t1
976  szr = -c11r / t1
977  szi = -c11i / t1
978  ELSE
979  cz = dlapy2( c22r, c22i )
980  IF( cz.LE.safmin ) THEN
981  cz = zero
982  szr = one
983  szi = zero
984  ELSE
985  tempr = c22r / cz
986  tempi = c22i / cz
987  t1 = dlapy2( cz, c21 )
988  cz = cz / t1
989  szr = -c21*tempr / t1
990  szi = c21*tempi / t1
991  END IF
992  END IF
993 *
994 * Compute Givens rotation on left
995 *
996 * ( CQ SQ )
997 * ( __ ) A or B
998 * ( -SQ CQ )
999 *
1000  an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1001  bn = abs( b11 ) + abs( b22 )
1002  wabs = abs( wr ) + abs( wi )
1003  IF( s1*an.GT.wabs*bn ) THEN
1004  cq = cz*b11
1005  sqr = szr*b22
1006  sqi = -szi*b22
1007  ELSE
1008  a1r = cz*a11 + szr*a12
1009  a1i = szi*a12
1010  a2r = cz*a21 + szr*a22
1011  a2i = szi*a22
1012  cq = dlapy2( a1r, a1i )
1013  IF( cq.LE.safmin ) THEN
1014  cq = zero
1015  sqr = one
1016  sqi = zero
1017  ELSE
1018  tempr = a1r / cq
1019  tempi = a1i / cq
1020  sqr = tempr*a2r + tempi*a2i
1021  sqi = tempi*a2r - tempr*a2i
1022  END IF
1023  END IF
1024  t1 = dlapy3( cq, sqr, sqi )
1025  cq = cq / t1
1026  sqr = sqr / t1
1027  sqi = sqi / t1
1028 *
1029 * Compute diagonal elements of QBZ
1030 *
1031  tempr = sqr*szr - sqi*szi
1032  tempi = sqr*szi + sqi*szr
1033  b1r = cq*cz*b11 + tempr*b22
1034  b1i = tempi*b22
1035  b1a = dlapy2( b1r, b1i )
1036  b2r = cq*cz*b22 + tempr*b11
1037  b2i = -tempi*b11
1038  b2a = dlapy2( b2r, b2i )
1039 *
1040 * Normalize so beta > 0, and Im( alpha1 ) > 0
1041 *
1042  beta( ilast-1 ) = b1a
1043  beta( ilast ) = b2a
1044  alphar( ilast-1 ) = ( wr*b1a )*s1inv
1045  alphai( ilast-1 ) = ( wi*b1a )*s1inv
1046  alphar( ilast ) = ( wr*b2a )*s1inv
1047  alphai( ilast ) = -( wi*b2a )*s1inv
1048 *
1049 * Step 3: Go to next block -- exit if finished.
1050 *
1051  ilast = ifirst - 1
1052  IF( ilast.LT.ilo )
1053  $ GO TO 380
1054 *
1055 * Reset counters
1056 *
1057  iiter = 0
1058  eshift = zero
1059  IF( .NOT.ilschr ) THEN
1060  ilastm = ilast
1061  IF( ifrstm.GT.ilast )
1062  $ ifrstm = ilo
1063  END IF
1064  GO TO 350
1065  ELSE
1066 *
1067 * Usual case: 3x3 or larger block, using Francis implicit
1068 * double-shift
1069 *
1070 * 2
1071 * Eigenvalue equation is w - c w + d = 0,
1072 *
1073 * -1 2 -1
1074 * so compute 1st column of (A B ) - c A B + d
1075 * using the formula in QZIT (from EISPACK)
1076 *
1077 * We assume that the block is at least 3x3
1078 *
1079  ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1080  $ ( bscale*t( ilast-1, ilast-1 ) )
1081  ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1082  $ ( bscale*t( ilast-1, ilast-1 ) )
1083  ad12 = ( ascale*h( ilast-1, ilast ) ) /
1084  $ ( bscale*t( ilast, ilast ) )
1085  ad22 = ( ascale*h( ilast, ilast ) ) /
1086  $ ( bscale*t( ilast, ilast ) )
1087  u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1088  ad11l = ( ascale*h( ifirst, ifirst ) ) /
1089  $ ( bscale*t( ifirst, ifirst ) )
1090  ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1091  $ ( bscale*t( ifirst, ifirst ) )
1092  ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1093  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1094  ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1095  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1096  ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1097  $ ( bscale*t( ifirst+1, ifirst+1 ) )
1098  u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1099 *
1100  v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1101  $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1102  v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1103  $ ( ad22-ad11l )+ad21*u12 )*ad21l
1104  v( 3 ) = ad32l*ad21l
1105 *
1106  istart = ifirst
1107 *
1108  CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1109  v( 1 ) = one
1110 *
1111 * Sweep
1112 *
1113  DO 290 j = istart, ilast - 2
1114 *
1115 * All but last elements: use 3x3 Householder transforms.
1116 *
1117 * Zero (j-1)st column of A
1118 *
1119  IF( j.GT.istart ) THEN
1120  v( 1 ) = h( j, j-1 )
1121  v( 2 ) = h( j+1, j-1 )
1122  v( 3 ) = h( j+2, j-1 )
1123 *
1124  CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1125  v( 1 ) = one
1126  h( j+1, j-1 ) = zero
1127  h( j+2, j-1 ) = zero
1128  END IF
1129 *
1130  DO 230 jc = j, ilastm
1131  temp = tau*( h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1132  $ h( j+2, jc ) )
1133  h( j, jc ) = h( j, jc ) - temp
1134  h( j+1, jc ) = h( j+1, jc ) - temp*v( 2 )
1135  h( j+2, jc ) = h( j+2, jc ) - temp*v( 3 )
1136  temp2 = tau*( t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1137  $ t( j+2, jc ) )
1138  t( j, jc ) = t( j, jc ) - temp2
1139  t( j+1, jc ) = t( j+1, jc ) - temp2*v( 2 )
1140  t( j+2, jc ) = t( j+2, jc ) - temp2*v( 3 )
1141  230 CONTINUE
1142  IF( ilq ) THEN
1143  DO 240 jr = 1, n
1144  temp = tau*( q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1145  $ q( jr, j+2 ) )
1146  q( jr, j ) = q( jr, j ) - temp
1147  q( jr, j+1 ) = q( jr, j+1 ) - temp*v( 2 )
1148  q( jr, j+2 ) = q( jr, j+2 ) - temp*v( 3 )
1149  240 CONTINUE
1150  END IF
1151 *
1152 * Zero j-th column of B (see DLAGBC for details)
1153 *
1154 * Swap rows to pivot
1155 *
1156  ilpivt = .false.
1157  temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1158  temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1159  IF( max( temp, temp2 ).LT.safmin ) THEN
1160  scale = zero
1161  u1 = one
1162  u2 = zero
1163  GO TO 250
1164  ELSE IF( temp.GE.temp2 ) THEN
1165  w11 = t( j+1, j+1 )
1166  w21 = t( j+2, j+1 )
1167  w12 = t( j+1, j+2 )
1168  w22 = t( j+2, j+2 )
1169  u1 = t( j+1, j )
1170  u2 = t( j+2, j )
1171  ELSE
1172  w21 = t( j+1, j+1 )
1173  w11 = t( j+2, j+1 )
1174  w22 = t( j+1, j+2 )
1175  w12 = t( j+2, j+2 )
1176  u2 = t( j+1, j )
1177  u1 = t( j+2, j )
1178  END IF
1179 *
1180 * Swap columns if nec.
1181 *
1182  IF( abs( w12 ).GT.abs( w11 ) ) THEN
1183  ilpivt = .true.
1184  temp = w12
1185  temp2 = w22
1186  w12 = w11
1187  w22 = w21
1188  w11 = temp
1189  w21 = temp2
1190  END IF
1191 *
1192 * LU-factor
1193 *
1194  temp = w21 / w11
1195  u2 = u2 - temp*u1
1196  w22 = w22 - temp*w12
1197  w21 = zero
1198 *
1199 * Compute SCALE
1200 *
1201  scale = one
1202  IF( abs( w22 ).LT.safmin ) THEN
1203  scale = zero
1204  u2 = one
1205  u1 = -w12 / w11
1206  GO TO 250
1207  END IF
1208  IF( abs( w22 ).LT.abs( u2 ) )
1209  $ scale = abs( w22 / u2 )
1210  IF( abs( w11 ).LT.abs( u1 ) )
1211  $ scale = min( scale, abs( w11 / u1 ) )
1212 *
1213 * Solve
1214 *
1215  u2 = ( scale*u2 ) / w22
1216  u1 = ( scale*u1-w12*u2 ) / w11
1217 *
1218  250 CONTINUE
1219  IF( ilpivt ) THEN
1220  temp = u2
1221  u2 = u1
1222  u1 = temp
1223  END IF
1224 *
1225 * Compute Householder Vector
1226 *
1227  t1 = sqrt( scale**2+u1**2+u2**2 )
1228  tau = one + scale / t1
1229  vs = -one / ( scale+t1 )
1230  v( 1 ) = one
1231  v( 2 ) = vs*u1
1232  v( 3 ) = vs*u2
1233 *
1234 * Apply transformations from the right.
1235 *
1236  DO 260 jr = ifrstm, min( j+3, ilast )
1237  temp = tau*( h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1238  $ h( jr, j+2 ) )
1239  h( jr, j ) = h( jr, j ) - temp
1240  h( jr, j+1 ) = h( jr, j+1 ) - temp*v( 2 )
1241  h( jr, j+2 ) = h( jr, j+2 ) - temp*v( 3 )
1242  260 CONTINUE
1243  DO 270 jr = ifrstm, j + 2
1244  temp = tau*( t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1245  $ t( jr, j+2 ) )
1246  t( jr, j ) = t( jr, j ) - temp
1247  t( jr, j+1 ) = t( jr, j+1 ) - temp*v( 2 )
1248  t( jr, j+2 ) = t( jr, j+2 ) - temp*v( 3 )
1249  270 CONTINUE
1250  IF( ilz ) THEN
1251  DO 280 jr = 1, n
1252  temp = tau*( z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1253  $ z( jr, j+2 ) )
1254  z( jr, j ) = z( jr, j ) - temp
1255  z( jr, j+1 ) = z( jr, j+1 ) - temp*v( 2 )
1256  z( jr, j+2 ) = z( jr, j+2 ) - temp*v( 3 )
1257  280 CONTINUE
1258  END IF
1259  t( j+1, j ) = zero
1260  t( j+2, j ) = zero
1261  290 CONTINUE
1262 *
1263 * Last elements: Use Givens rotations
1264 *
1265 * Rotations from the left
1266 *
1267  j = ilast - 1
1268  temp = h( j, j-1 )
1269  CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1270  h( j+1, j-1 ) = zero
1271 *
1272  DO 300 jc = j, ilastm
1273  temp = c*h( j, jc ) + s*h( j+1, jc )
1274  h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1275  h( j, jc ) = temp
1276  temp2 = c*t( j, jc ) + s*t( j+1, jc )
1277  t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1278  t( j, jc ) = temp2
1279  300 CONTINUE
1280  IF( ilq ) THEN
1281  DO 310 jr = 1, n
1282  temp = c*q( jr, j ) + s*q( jr, j+1 )
1283  q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1284  q( jr, j ) = temp
1285  310 CONTINUE
1286  END IF
1287 *
1288 * Rotations from the right.
1289 *
1290  temp = t( j+1, j+1 )
1291  CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1292  t( j+1, j ) = zero
1293 *
1294  DO 320 jr = ifrstm, ilast
1295  temp = c*h( jr, j+1 ) + s*h( jr, j )
1296  h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1297  h( jr, j+1 ) = temp
1298  320 CONTINUE
1299  DO 330 jr = ifrstm, ilast - 1
1300  temp = c*t( jr, j+1 ) + s*t( jr, j )
1301  t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1302  t( jr, j+1 ) = temp
1303  330 CONTINUE
1304  IF( ilz ) THEN
1305  DO 340 jr = 1, n
1306  temp = c*z( jr, j+1 ) + s*z( jr, j )
1307  z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1308  z( jr, j+1 ) = temp
1309  340 CONTINUE
1310  END IF
1311 *
1312 * End of Double-Shift code
1313 *
1314  END IF
1315 *
1316  GO TO 350
1317 *
1318 * End of iteration loop
1319 *
1320  350 CONTINUE
1321  360 CONTINUE
1322 *
1323 * Drop-through = non-convergence
1324 *
1325  info = ilast
1326  GO TO 420
1327 *
1328 * Successful completion of all QZ steps
1329 *
1330  380 CONTINUE
1331 *
1332 * Set Eigenvalues 1:ILO-1
1333 *
1334  DO 410 j = 1, ilo - 1
1335  IF( t( j, j ).LT.zero ) THEN
1336  IF( ilschr ) THEN
1337  DO 390 jr = 1, j
1338  h( jr, j ) = -h( jr, j )
1339  t( jr, j ) = -t( jr, j )
1340  390 CONTINUE
1341  ELSE
1342  h( j, j ) = -h( j, j )
1343  t( j, j ) = -t( j, j )
1344  END IF
1345  IF( ilz ) THEN
1346  DO 400 jr = 1, n
1347  z( jr, j ) = -z( jr, j )
1348  400 CONTINUE
1349  END IF
1350  END IF
1351  alphar( j ) = h( j, j )
1352  alphai( j ) = zero
1353  beta( j ) = t( j, j )
1354  410 CONTINUE
1355 *
1356 * Normal Termination
1357 *
1358  info = 0
1359 *
1360 * Exit (other than argument error) -- return optimal workspace size
1361 *
1362  420 CONTINUE
1363  work( 1 ) = dble( n )
1364  RETURN
1365 *
1366 * End of DHGEQZ
1367 *
1368  END
dlarfg
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
dlag2
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
dlasv2
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
dhgeqz
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
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
dlaset
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