1 SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
2 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
3 $ LWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2009 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER COMPQ, COMPZ, JOB
12 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ),
16 $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
17 $ WORK( * ), Z( LDZ, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
24 * where H is an upper Hessenberg matrix and T is upper triangular,
25 * using the double-shift QZ method.
26 * Matrix pairs of this type are produced by the reduction to
27 * generalized upper Hessenberg form of a real matrix pair (A,B):
28 *
29 * A = Q1*H*Z1**T, B = Q1*T*Z1**T,
30 *
31 * as computed by DGGHRD.
32 *
33 * If JOB='S', then the Hessenberg-triangular pair (H,T) is
34 * also reduced to generalized Schur form,
35 *
36 * H = Q*S*Z**T, T = Q*P*Z**T,
37 *
38 * where Q and Z are orthogonal matrices, P is an upper triangular
39 * matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
40 * diagonal blocks.
41 *
42 * The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
43 * (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
44 * eigenvalues.
45 *
46 * Additionally, the 2-by-2 upper triangular diagonal blocks of P
47 * corresponding to 2-by-2 blocks of S are reduced to positive diagonal
48 * form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
49 * P(j,j) > 0, and P(j+1,j+1) > 0.
50 *
51 * Optionally, the orthogonal matrix Q from the generalized Schur
52 * factorization may be postmultiplied into an input matrix Q1, and the
53 * orthogonal matrix Z may be postmultiplied into an input matrix Z1.
54 * If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
55 * the matrix pair (A,B) to generalized upper Hessenberg form, then the
56 * output matrices Q1*Q and Z1*Z are the orthogonal factors from the
57 * generalized Schur factorization of (A,B):
58 *
59 * A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T.
60 *
61 * To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
62 * of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
63 * complex and beta real.
64 * If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
65 * generalized nonsymmetric eigenvalue problem (GNEP)
66 * A*x = lambda*B*x
67 * and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
68 * alternate form of the GNEP
69 * mu*A*y = B*y.
70 * Real eigenvalues can be read directly from the generalized Schur
71 * form:
72 * alpha = S(i,i), beta = P(i,i).
73 *
74 * Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
75 * Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
76 * pp. 241--256.
77 *
78 * Arguments
79 * =========
80 *
81 * JOB (input) CHARACTER*1
82 * = 'E': Compute eigenvalues only;
83 * = 'S': Compute eigenvalues and the Schur form.
84 *
85 * COMPQ (input) CHARACTER*1
86 * = 'N': Left Schur vectors (Q) are not computed;
87 * = 'I': Q is initialized to the unit matrix and the matrix Q
88 * of left Schur vectors of (H,T) is returned;
89 * = 'V': Q must contain an orthogonal matrix Q1 on entry and
90 * the product Q1*Q is returned.
91 *
92 * COMPZ (input) CHARACTER*1
93 * = 'N': Right Schur vectors (Z) are not computed;
94 * = 'I': Z is initialized to the unit matrix and the matrix Z
95 * of right Schur vectors of (H,T) is returned;
96 * = 'V': Z must contain an orthogonal matrix Z1 on entry and
97 * the product Z1*Z is returned.
98 *
99 * N (input) INTEGER
100 * The order of the matrices H, T, Q, and Z. N >= 0.
101 *
102 * ILO (input) INTEGER
103 * IHI (input) INTEGER
104 * ILO and IHI mark the rows and columns of H which are in
105 * Hessenberg form. It is assumed that A is already upper
106 * triangular in rows and columns 1:ILO-1 and IHI+1:N.
107 * If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
108 *
109 * H (input/output) DOUBLE PRECISION array, dimension (LDH, N)
110 * On entry, the N-by-N upper Hessenberg matrix H.
111 * On exit, if JOB = 'S', H contains the upper quasi-triangular
112 * matrix S from the generalized Schur factorization.
113 * If JOB = 'E', the diagonal blocks of H match those of S, but
114 * the rest of H is unspecified.
115 *
116 * LDH (input) INTEGER
117 * The leading dimension of the array H. LDH >= max( 1, N ).
118 *
119 * T (input/output) DOUBLE PRECISION array, dimension (LDT, N)
120 * On entry, the N-by-N upper triangular matrix T.
121 * On exit, if JOB = 'S', T contains the upper triangular
122 * matrix P from the generalized Schur factorization;
123 * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
124 * are reduced to positive diagonal form, i.e., if H(j+1,j) is
125 * non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
126 * T(j+1,j+1) > 0.
127 * If JOB = 'E', the diagonal blocks of T match those of P, but
128 * the rest of T is unspecified.
129 *
130 * LDT (input) INTEGER
131 * The leading dimension of the array T. LDT >= max( 1, N ).
132 *
133 * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
134 * The real parts of each scalar alpha defining an eigenvalue
135 * of GNEP.
136 *
137 * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
138 * The imaginary parts of each scalar alpha defining an
139 * eigenvalue of GNEP.
140 * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
141 * positive, then the j-th and (j+1)-st eigenvalues are a
142 * complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
143 *
144 * BETA (output) DOUBLE PRECISION array, dimension (N)
145 * The scalars beta that define the eigenvalues of GNEP.
146 * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
147 * beta = BETA(j) represent the j-th eigenvalue of the matrix
148 * pair (A,B), in one of the forms lambda = alpha/beta or
149 * mu = beta/alpha. Since either lambda or mu may overflow,
150 * they should not, in general, be computed.
151 *
152 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
153 * On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in
154 * the reduction of (A,B) to generalized Hessenberg form.
155 * On exit, if COMPZ = 'I', the orthogonal matrix of left Schur
156 * vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix
157 * of left Schur vectors of (A,B).
158 * Not referenced if COMPZ = 'N'.
159 *
160 * LDQ (input) INTEGER
161 * The leading dimension of the array Q. LDQ >= 1.
162 * If COMPQ='V' or 'I', then LDQ >= N.
163 *
164 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
165 * On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
166 * the reduction of (A,B) to generalized Hessenberg form.
167 * On exit, if COMPZ = 'I', the orthogonal matrix of
168 * right Schur vectors of (H,T), and if COMPZ = 'V', the
169 * orthogonal matrix of right Schur vectors of (A,B).
170 * Not referenced if COMPZ = 'N'.
171 *
172 * LDZ (input) INTEGER
173 * The leading dimension of the array Z. LDZ >= 1.
174 * If COMPZ='V' or 'I', then LDZ >= N.
175 *
176 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
177 * On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
178 *
179 * LWORK (input) INTEGER
180 * The dimension of the array WORK. LWORK >= max(1,N).
181 *
182 * If LWORK = -1, then a workspace query is assumed; the routine
183 * only calculates the optimal size of the WORK array, returns
184 * this value as the first entry of the WORK array, and no error
185 * message related to LWORK is issued by XERBLA.
186 *
187 * INFO (output) INTEGER
188 * = 0: successful exit
189 * < 0: if INFO = -i, the i-th argument had an illegal value
190 * = 1,...,N: the QZ iteration did not converge. (H,T) is not
191 * in Schur form, but ALPHAR(i), ALPHAI(i), and
192 * BETA(i), i=INFO+1,...,N should be correct.
193 * = N+1,...,2*N: the shift calculation failed. (H,T) is not
194 * in Schur form, but ALPHAR(i), ALPHAI(i), and
195 * BETA(i), i=INFO-N+1,...,N should be correct.
196 *
197 * Further Details
198 * ===============
199 *
200 * Iteration counters:
201 *
202 * JITER -- counts iterations.
203 * IITER -- counts iterations run since ILAST was last
204 * changed. This is therefore reset only when a 1-by-1 or
205 * 2-by-2 block deflates off the bottom.
206 *
207 * =====================================================================
208 *
209 * .. Parameters ..
210 * $ SAFETY = 1.0E+0 )
211 DOUBLE PRECISION HALF, ZERO, ONE, SAFETY
212 PARAMETER ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0,
213 $ SAFETY = 1.0D+2 )
214 * ..
215 * .. Local Scalars ..
216 LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
217 $ LQUERY
218 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
219 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
220 $ JR, MAXIT
221 DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
222 $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
223 $ AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
224 $ B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
225 $ BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
226 $ CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
227 $ SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
228 $ TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
229 $ U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
230 $ WR2
231 * ..
232 * .. Local Arrays ..
233 DOUBLE PRECISION V( 3 )
234 * ..
235 * .. External Functions ..
236 LOGICAL LSAME
237 DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3
238 EXTERNAL LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3
239 * ..
240 * .. External Subroutines ..
241 EXTERNAL DLAG2, DLARFG, DLARTG, DLASET, DLASV2, DROT,
242 $ XERBLA
243 * ..
244 * .. Intrinsic Functions ..
245 INTRINSIC ABS, DBLE, MAX, MIN, SQRT
246 * ..
247 * .. Executable Statements ..
248 *
249 * Decode JOB, COMPQ, COMPZ
250 *
251 IF( LSAME( JOB, 'E' ) ) THEN
252 ILSCHR = .FALSE.
253 ISCHUR = 1
254 ELSE IF( LSAME( JOB, 'S' ) ) THEN
255 ILSCHR = .TRUE.
256 ISCHUR = 2
257 ELSE
258 ISCHUR = 0
259 END IF
260 *
261 IF( LSAME( COMPQ, 'N' ) ) THEN
262 ILQ = .FALSE.
263 ICOMPQ = 1
264 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
265 ILQ = .TRUE.
266 ICOMPQ = 2
267 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
268 ILQ = .TRUE.
269 ICOMPQ = 3
270 ELSE
271 ICOMPQ = 0
272 END IF
273 *
274 IF( LSAME( COMPZ, 'N' ) ) THEN
275 ILZ = .FALSE.
276 ICOMPZ = 1
277 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
278 ILZ = .TRUE.
279 ICOMPZ = 2
280 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
281 ILZ = .TRUE.
282 ICOMPZ = 3
283 ELSE
284 ICOMPZ = 0
285 END IF
286 *
287 * Check Argument Values
288 *
289 INFO = 0
290 WORK( 1 ) = MAX( 1, N )
291 LQUERY = ( LWORK.EQ.-1 )
292 IF( ISCHUR.EQ.0 ) THEN
293 INFO = -1
294 ELSE IF( ICOMPQ.EQ.0 ) THEN
295 INFO = -2
296 ELSE IF( ICOMPZ.EQ.0 ) THEN
297 INFO = -3
298 ELSE IF( N.LT.0 ) THEN
299 INFO = -4
300 ELSE IF( ILO.LT.1 ) THEN
301 INFO = -5
302 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
303 INFO = -6
304 ELSE IF( LDH.LT.N ) THEN
305 INFO = -8
306 ELSE IF( LDT.LT.N ) THEN
307 INFO = -10
308 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
309 INFO = -15
310 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
311 INFO = -17
312 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
313 INFO = -19
314 END IF
315 IF( INFO.NE.0 ) THEN
316 CALL XERBLA( 'DHGEQZ', -INFO )
317 RETURN
318 ELSE IF( LQUERY ) THEN
319 RETURN
320 END IF
321 *
322 * Quick return if possible
323 *
324 IF( N.LE.0 ) THEN
325 WORK( 1 ) = DBLE( 1 )
326 RETURN
327 END IF
328 *
329 * Initialize Q and Z
330 *
331 IF( ICOMPQ.EQ.3 )
332 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
333 IF( ICOMPZ.EQ.3 )
334 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
335 *
336 * Machine Constants
337 *
338 IN = IHI + 1 - ILO
339 SAFMIN = DLAMCH( 'S' )
340 SAFMAX = ONE / SAFMIN
341 ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
342 ANORM = DLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK )
343 BNORM = DLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK )
344 ATOL = MAX( SAFMIN, ULP*ANORM )
345 BTOL = MAX( SAFMIN, ULP*BNORM )
346 ASCALE = ONE / MAX( SAFMIN, ANORM )
347 BSCALE = ONE / MAX( SAFMIN, BNORM )
348 *
349 * Set Eigenvalues IHI+1:N
350 *
351 DO 30 J = IHI + 1, N
352 IF( T( J, J ).LT.ZERO ) THEN
353 IF( ILSCHR ) THEN
354 DO 10 JR = 1, J
355 H( JR, J ) = -H( JR, J )
356 T( JR, J ) = -T( JR, J )
357 10 CONTINUE
358 ELSE
359 H( J, J ) = -H( J, J )
360 T( J, J ) = -T( J, J )
361 END IF
362 IF( ILZ ) THEN
363 DO 20 JR = 1, N
364 Z( JR, J ) = -Z( JR, J )
365 20 CONTINUE
366 END IF
367 END IF
368 ALPHAR( J ) = H( J, J )
369 ALPHAI( J ) = ZERO
370 BETA( J ) = T( J, J )
371 30 CONTINUE
372 *
373 * If IHI < ILO, skip QZ steps
374 *
375 IF( IHI.LT.ILO )
376 $ GO TO 380
377 *
378 * MAIN QZ ITERATION LOOP
379 *
380 * Initialize dynamic indices
381 *
382 * Eigenvalues ILAST+1:N have been found.
383 * Column operations modify rows IFRSTM:whatever.
384 * Row operations modify columns whatever:ILASTM.
385 *
386 * If only eigenvalues are being computed, then
387 * IFRSTM is the row of the last splitting row above row ILAST;
388 * this is always at least ILO.
389 * IITER counts iterations since the last eigenvalue was found,
390 * to tell when to use an extraordinary shift.
391 * MAXIT is the maximum number of QZ sweeps allowed.
392 *
393 ILAST = IHI
394 IF( ILSCHR ) THEN
395 IFRSTM = 1
396 ILASTM = N
397 ELSE
398 IFRSTM = ILO
399 ILASTM = IHI
400 END IF
401 IITER = 0
402 ESHIFT = ZERO
403 MAXIT = 30*( IHI-ILO+1 )
404 *
405 DO 360 JITER = 1, MAXIT
406 *
407 * Split the matrix if possible.
408 *
409 * Two tests:
410 * 1: H(j,j-1)=0 or j=ILO
411 * 2: T(j,j)=0
412 *
413 IF( ILAST.EQ.ILO ) THEN
414 *
415 * Special case: j=ILAST
416 *
417 GO TO 80
418 ELSE
419 IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
420 H( ILAST, ILAST-1 ) = ZERO
421 GO TO 80
422 END IF
423 END IF
424 *
425 IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
426 T( ILAST, ILAST ) = ZERO
427 GO TO 70
428 END IF
429 *
430 * General case: j<ILAST
431 *
432 DO 60 J = ILAST - 1, ILO, -1
433 *
434 * Test 1: for H(j,j-1)=0 or j=ILO
435 *
436 IF( J.EQ.ILO ) THEN
437 ILAZRO = .TRUE.
438 ELSE
439 IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN
440 H( J, J-1 ) = ZERO
441 ILAZRO = .TRUE.
442 ELSE
443 ILAZRO = .FALSE.
444 END IF
445 END IF
446 *
447 * Test 2: for T(j,j)=0
448 *
449 IF( ABS( T( J, J ) ).LT.BTOL ) THEN
450 T( J, J ) = ZERO
451 *
452 * Test 1a: Check for 2 consecutive small subdiagonals in A
453 *
454 ILAZR2 = .FALSE.
455 IF( .NOT.ILAZRO ) THEN
456 TEMP = ABS( H( J, J-1 ) )
457 TEMP2 = ABS( H( J, J ) )
458 TEMPR = MAX( TEMP, TEMP2 )
459 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
460 TEMP = TEMP / TEMPR
461 TEMP2 = TEMP2 / TEMPR
462 END IF
463 IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2*
464 $ ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
465 END IF
466 *
467 * If both tests pass (1 & 2), i.e., the leading diagonal
468 * element of B in the block is zero, split a 1x1 block off
469 * at the top. (I.e., at the J-th row/column) The leading
470 * diagonal element of the remainder can also be zero, so
471 * this may have to be done repeatedly.
472 *
473 IF( ILAZRO .OR. ILAZR2 ) THEN
474 DO 40 JCH = J, ILAST - 1
475 TEMP = H( JCH, JCH )
476 CALL DLARTG( TEMP, H( JCH+1, JCH ), C, S,
477 $ H( JCH, JCH ) )
478 H( JCH+1, JCH ) = ZERO
479 CALL DROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
480 $ H( JCH+1, JCH+1 ), LDH, C, S )
481 CALL DROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
482 $ T( JCH+1, JCH+1 ), LDT, C, S )
483 IF( ILQ )
484 $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
485 $ C, S )
486 IF( ILAZR2 )
487 $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
488 ILAZR2 = .FALSE.
489 IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
490 IF( JCH+1.GE.ILAST ) THEN
491 GO TO 80
492 ELSE
493 IFIRST = JCH + 1
494 GO TO 110
495 END IF
496 END IF
497 T( JCH+1, JCH+1 ) = ZERO
498 40 CONTINUE
499 GO TO 70
500 ELSE
501 *
502 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
503 * Then process as in the case T(ILAST,ILAST)=0
504 *
505 DO 50 JCH = J, ILAST - 1
506 TEMP = T( JCH, JCH+1 )
507 CALL DLARTG( TEMP, T( JCH+1, JCH+1 ), C, S,
508 $ T( JCH, JCH+1 ) )
509 T( JCH+1, JCH+1 ) = ZERO
510 IF( JCH.LT.ILASTM-1 )
511 $ CALL DROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
512 $ T( JCH+1, JCH+2 ), LDT, C, S )
513 CALL DROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
514 $ H( JCH+1, JCH-1 ), LDH, C, S )
515 IF( ILQ )
516 $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
517 $ C, S )
518 TEMP = H( JCH+1, JCH )
519 CALL DLARTG( TEMP, H( JCH+1, JCH-1 ), C, S,
520 $ H( JCH+1, JCH ) )
521 H( JCH+1, JCH-1 ) = ZERO
522 CALL DROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
523 $ H( IFRSTM, JCH-1 ), 1, C, S )
524 CALL DROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
525 $ T( IFRSTM, JCH-1 ), 1, C, S )
526 IF( ILZ )
527 $ CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
528 $ C, S )
529 50 CONTINUE
530 GO TO 70
531 END IF
532 ELSE IF( ILAZRO ) THEN
533 *
534 * Only test 1 passed -- work on J:ILAST
535 *
536 IFIRST = J
537 GO TO 110
538 END IF
539 *
540 * Neither test passed -- try next J
541 *
542 60 CONTINUE
543 *
544 * (Drop-through is "impossible")
545 *
546 INFO = N + 1
547 GO TO 420
548 *
549 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
550 * 1x1 block.
551 *
552 70 CONTINUE
553 TEMP = H( ILAST, ILAST )
554 CALL DLARTG( TEMP, H( ILAST, ILAST-1 ), C, S,
555 $ H( ILAST, ILAST ) )
556 H( ILAST, ILAST-1 ) = ZERO
557 CALL DROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
558 $ H( IFRSTM, ILAST-1 ), 1, C, S )
559 CALL DROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
560 $ T( IFRSTM, ILAST-1 ), 1, C, S )
561 IF( ILZ )
562 $ CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
563 *
564 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
565 * and BETA
566 *
567 80 CONTINUE
568 IF( T( ILAST, ILAST ).LT.ZERO ) THEN
569 IF( ILSCHR ) THEN
570 DO 90 J = IFRSTM, ILAST
571 H( J, ILAST ) = -H( J, ILAST )
572 T( J, ILAST ) = -T( J, ILAST )
573 90 CONTINUE
574 ELSE
575 H( ILAST, ILAST ) = -H( ILAST, ILAST )
576 T( ILAST, ILAST ) = -T( ILAST, ILAST )
577 END IF
578 IF( ILZ ) THEN
579 DO 100 J = 1, N
580 Z( J, ILAST ) = -Z( J, ILAST )
581 100 CONTINUE
582 END IF
583 END IF
584 ALPHAR( ILAST ) = H( ILAST, ILAST )
585 ALPHAI( ILAST ) = ZERO
586 BETA( ILAST ) = T( ILAST, ILAST )
587 *
588 * Go to next block -- exit if finished.
589 *
590 ILAST = ILAST - 1
591 IF( ILAST.LT.ILO )
592 $ GO TO 380
593 *
594 * Reset counters
595 *
596 IITER = 0
597 ESHIFT = ZERO
598 IF( .NOT.ILSCHR ) THEN
599 ILASTM = ILAST
600 IF( IFRSTM.GT.ILAST )
601 $ IFRSTM = ILO
602 END IF
603 GO TO 350
604 *
605 * QZ step
606 *
607 * This iteration only involves rows/columns IFIRST:ILAST. We
608 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
609 *
610 110 CONTINUE
611 IITER = IITER + 1
612 IF( .NOT.ILSCHR ) THEN
613 IFRSTM = IFIRST
614 END IF
615 *
616 * Compute single shifts.
617 *
618 * At this point, IFIRST < ILAST, and the diagonal elements of
619 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
620 * magnitude)
621 *
622 IF( ( IITER / 10 )*10.EQ.IITER ) THEN
623 *
624 * Exceptional shift. Chosen for no particularly good reason.
625 * (Single shift only.)
626 *
627 IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
628 $ ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
629 ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /
630 $ T( ILAST-1, ILAST-1 )
631 ELSE
632 ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
633 END IF
634 S1 = ONE
635 WR = ESHIFT
636 *
637 ELSE
638 *
639 * Shifts based on the generalized eigenvalues of the
640 * bottom-right 2x2 block of A and B. The first eigenvalue
641 * returned by DLAG2 is the Wilkinson shift (AEP p.512),
642 *
643 CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
644 $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
645 $ S2, WR, WR2, WI )
646 *
647 TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
648 IF( WI.NE.ZERO )
649 $ GO TO 200
650 END IF
651 *
652 * Fiddle with shift to avoid overflow
653 *
654 TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
655 IF( S1.GT.TEMP ) THEN
656 SCALE = TEMP / S1
657 ELSE
658 SCALE = ONE
659 END IF
660 *
661 TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
662 IF( ABS( WR ).GT.TEMP )
663 $ SCALE = MIN( SCALE, TEMP / ABS( WR ) )
664 S1 = SCALE*S1
665 WR = SCALE*WR
666 *
667 * Now check for two consecutive small subdiagonals.
668 *
669 DO 120 J = ILAST - 1, IFIRST + 1, -1
670 ISTART = J
671 TEMP = ABS( S1*H( J, J-1 ) )
672 TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) )
673 TEMPR = MAX( TEMP, TEMP2 )
674 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
675 TEMP = TEMP / TEMPR
676 TEMP2 = TEMP2 / TEMPR
677 END IF
678 IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
679 $ TEMP2 )GO TO 130
680 120 CONTINUE
681 *
682 ISTART = IFIRST
683 130 CONTINUE
684 *
685 * Do an implicit single-shift QZ sweep.
686 *
687 * Initial Q
688 *
689 TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART )
690 TEMP2 = S1*H( ISTART+1, ISTART )
691 CALL DLARTG( TEMP, TEMP2, C, S, TEMPR )
692 *
693 * Sweep
694 *
695 DO 190 J = ISTART, ILAST - 1
696 IF( J.GT.ISTART ) THEN
697 TEMP = H( J, J-1 )
698 CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
699 H( J+1, J-1 ) = ZERO
700 END IF
701 *
702 DO 140 JC = J, ILASTM
703 TEMP = C*H( J, JC ) + S*H( J+1, JC )
704 H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
705 H( J, JC ) = TEMP
706 TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
707 T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
708 T( J, JC ) = TEMP2
709 140 CONTINUE
710 IF( ILQ ) THEN
711 DO 150 JR = 1, N
712 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
713 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
714 Q( JR, J ) = TEMP
715 150 CONTINUE
716 END IF
717 *
718 TEMP = T( J+1, J+1 )
719 CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
720 T( J+1, J ) = ZERO
721 *
722 DO 160 JR = IFRSTM, MIN( J+2, ILAST )
723 TEMP = C*H( JR, J+1 ) + S*H( JR, J )
724 H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
725 H( JR, J+1 ) = TEMP
726 160 CONTINUE
727 DO 170 JR = IFRSTM, J
728 TEMP = C*T( JR, J+1 ) + S*T( JR, J )
729 T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
730 T( JR, J+1 ) = TEMP
731 170 CONTINUE
732 IF( ILZ ) THEN
733 DO 180 JR = 1, N
734 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
735 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
736 Z( JR, J+1 ) = TEMP
737 180 CONTINUE
738 END IF
739 190 CONTINUE
740 *
741 GO TO 350
742 *
743 * Use Francis double-shift
744 *
745 * Note: the Francis double-shift should work with real shifts,
746 * but only if the block is at least 3x3.
747 * This code may break if this point is reached with
748 * a 2x2 block with real eigenvalues.
749 *
750 200 CONTINUE
751 IF( IFIRST+1.EQ.ILAST ) THEN
752 *
753 * Special case -- 2x2 block with complex eigenvectors
754 *
755 * Step 1: Standardize, that is, rotate so that
756 *
757 * ( B11 0 )
758 * B = ( ) with B11 non-negative.
759 * ( 0 B22 )
760 *
761 CALL DLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ),
762 $ T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
763 *
764 IF( B11.LT.ZERO ) THEN
765 CR = -CR
766 SR = -SR
767 B11 = -B11
768 B22 = -B22
769 END IF
770 *
771 CALL DROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,
772 $ H( ILAST, ILAST-1 ), LDH, CL, SL )
773 CALL DROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
774 $ H( IFRSTM, ILAST ), 1, CR, SR )
775 *
776 IF( ILAST.LT.ILASTM )
777 $ CALL DROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
778 $ T( ILAST, ILAST+1 ), LDT, CL, SL )
779 IF( IFRSTM.LT.ILAST-1 )
780 $ CALL DROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
781 $ T( IFRSTM, ILAST ), 1, CR, SR )
782 *
783 IF( ILQ )
784 $ CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
785 $ SL )
786 IF( ILZ )
787 $ CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
788 $ SR )
789 *
790 T( ILAST-1, ILAST-1 ) = B11
791 T( ILAST-1, ILAST ) = ZERO
792 T( ILAST, ILAST-1 ) = ZERO
793 T( ILAST, ILAST ) = B22
794 *
795 * If B22 is negative, negate column ILAST
796 *
797 IF( B22.LT.ZERO ) THEN
798 DO 210 J = IFRSTM, ILAST
799 H( J, ILAST ) = -H( J, ILAST )
800 T( J, ILAST ) = -T( J, ILAST )
801 210 CONTINUE
802 *
803 IF( ILZ ) THEN
804 DO 220 J = 1, N
805 Z( J, ILAST ) = -Z( J, ILAST )
806 220 CONTINUE
807 END IF
808 END IF
809 *
810 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
811 *
812 * Recompute shift
813 *
814 CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
815 $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
816 $ TEMP, WR, TEMP2, WI )
817 *
818 * If standardization has perturbed the shift onto real line,
819 * do another (real single-shift) QR step.
820 *
821 IF( WI.EQ.ZERO )
822 $ GO TO 350
823 S1INV = ONE / S1
824 *
825 * Do EISPACK (QZVAL) computation of alpha and beta
826 *
827 A11 = H( ILAST-1, ILAST-1 )
828 A21 = H( ILAST, ILAST-1 )
829 A12 = H( ILAST-1, ILAST )
830 A22 = H( ILAST, ILAST )
831 *
832 * Compute complex Givens rotation on right
833 * (Assume some element of C = (sA - wB) > unfl )
834 * __
835 * (sA - wB) ( CZ -SZ )
836 * ( SZ CZ )
837 *
838 C11R = S1*A11 - WR*B11
839 C11I = -WI*B11
840 C12 = S1*A12
841 C21 = S1*A21
842 C22R = S1*A22 - WR*B22
843 C22I = -WI*B22
844 *
845 IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
846 $ ABS( C22R )+ABS( C22I ) ) THEN
847 T1 = DLAPY3( C12, C11R, C11I )
848 CZ = C12 / T1
849 SZR = -C11R / T1
850 SZI = -C11I / T1
851 ELSE
852 CZ = DLAPY2( C22R, C22I )
853 IF( CZ.LE.SAFMIN ) THEN
854 CZ = ZERO
855 SZR = ONE
856 SZI = ZERO
857 ELSE
858 TEMPR = C22R / CZ
859 TEMPI = C22I / CZ
860 T1 = DLAPY2( CZ, C21 )
861 CZ = CZ / T1
862 SZR = -C21*TEMPR / T1
863 SZI = C21*TEMPI / T1
864 END IF
865 END IF
866 *
867 * Compute Givens rotation on left
868 *
869 * ( CQ SQ )
870 * ( __ ) A or B
871 * ( -SQ CQ )
872 *
873 AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
874 BN = ABS( B11 ) + ABS( B22 )
875 WABS = ABS( WR ) + ABS( WI )
876 IF( S1*AN.GT.WABS*BN ) THEN
877 CQ = CZ*B11
878 SQR = SZR*B22
879 SQI = -SZI*B22
880 ELSE
881 A1R = CZ*A11 + SZR*A12
882 A1I = SZI*A12
883 A2R = CZ*A21 + SZR*A22
884 A2I = SZI*A22
885 CQ = DLAPY2( A1R, A1I )
886 IF( CQ.LE.SAFMIN ) THEN
887 CQ = ZERO
888 SQR = ONE
889 SQI = ZERO
890 ELSE
891 TEMPR = A1R / CQ
892 TEMPI = A1I / CQ
893 SQR = TEMPR*A2R + TEMPI*A2I
894 SQI = TEMPI*A2R - TEMPR*A2I
895 END IF
896 END IF
897 T1 = DLAPY3( CQ, SQR, SQI )
898 CQ = CQ / T1
899 SQR = SQR / T1
900 SQI = SQI / T1
901 *
902 * Compute diagonal elements of QBZ
903 *
904 TEMPR = SQR*SZR - SQI*SZI
905 TEMPI = SQR*SZI + SQI*SZR
906 B1R = CQ*CZ*B11 + TEMPR*B22
907 B1I = TEMPI*B22
908 B1A = DLAPY2( B1R, B1I )
909 B2R = CQ*CZ*B22 + TEMPR*B11
910 B2I = -TEMPI*B11
911 B2A = DLAPY2( B2R, B2I )
912 *
913 * Normalize so beta > 0, and Im( alpha1 ) > 0
914 *
915 BETA( ILAST-1 ) = B1A
916 BETA( ILAST ) = B2A
917 ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
918 ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
919 ALPHAR( ILAST ) = ( WR*B2A )*S1INV
920 ALPHAI( ILAST ) = -( WI*B2A )*S1INV
921 *
922 * Step 3: Go to next block -- exit if finished.
923 *
924 ILAST = IFIRST - 1
925 IF( ILAST.LT.ILO )
926 $ GO TO 380
927 *
928 * Reset counters
929 *
930 IITER = 0
931 ESHIFT = ZERO
932 IF( .NOT.ILSCHR ) THEN
933 ILASTM = ILAST
934 IF( IFRSTM.GT.ILAST )
935 $ IFRSTM = ILO
936 END IF
937 GO TO 350
938 ELSE
939 *
940 * Usual case: 3x3 or larger block, using Francis implicit
941 * double-shift
942 *
943 * 2
944 * Eigenvalue equation is w - c w + d = 0,
945 *
946 * -1 2 -1
947 * so compute 1st column of (A B ) - c A B + d
948 * using the formula in QZIT (from EISPACK)
949 *
950 * We assume that the block is at least 3x3
951 *
952 AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
953 $ ( BSCALE*T( ILAST-1, ILAST-1 ) )
954 AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
955 $ ( BSCALE*T( ILAST-1, ILAST-1 ) )
956 AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
957 $ ( BSCALE*T( ILAST, ILAST ) )
958 AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
959 $ ( BSCALE*T( ILAST, ILAST ) )
960 U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST )
961 AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) /
962 $ ( BSCALE*T( IFIRST, IFIRST ) )
963 AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) /
964 $ ( BSCALE*T( IFIRST, IFIRST ) )
965 AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) /
966 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
967 AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) /
968 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
969 AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) /
970 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
971 U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 )
972 *
973 V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 +
974 $ AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L
975 V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )-
976 $ ( AD22-AD11L )+AD21*U12 )*AD21L
977 V( 3 ) = AD32L*AD21L
978 *
979 ISTART = IFIRST
980 *
981 CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU )
982 V( 1 ) = ONE
983 *
984 * Sweep
985 *
986 DO 290 J = ISTART, ILAST - 2
987 *
988 * All but last elements: use 3x3 Householder transforms.
989 *
990 * Zero (j-1)st column of A
991 *
992 IF( J.GT.ISTART ) THEN
993 V( 1 ) = H( J, J-1 )
994 V( 2 ) = H( J+1, J-1 )
995 V( 3 ) = H( J+2, J-1 )
996 *
997 CALL DLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU )
998 V( 1 ) = ONE
999 H( J+1, J-1 ) = ZERO
1000 H( J+2, J-1 ) = ZERO
1001 END IF
1002 *
1003 DO 230 JC = J, ILASTM
1004 TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )*
1005 $ H( J+2, JC ) )
1006 H( J, JC ) = H( J, JC ) - TEMP
1007 H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 )
1008 H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 )
1009 TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )*
1010 $ T( J+2, JC ) )
1011 T( J, JC ) = T( J, JC ) - TEMP2
1012 T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 )
1013 T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 )
1014 230 CONTINUE
1015 IF( ILQ ) THEN
1016 DO 240 JR = 1, N
1017 TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )*
1018 $ Q( JR, J+2 ) )
1019 Q( JR, J ) = Q( JR, J ) - TEMP
1020 Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 )
1021 Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 )
1022 240 CONTINUE
1023 END IF
1024 *
1025 * Zero j-th column of B (see DLAGBC for details)
1026 *
1027 * Swap rows to pivot
1028 *
1029 ILPIVT = .FALSE.
1030 TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) )
1031 TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) )
1032 IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN
1033 SCALE = ZERO
1034 U1 = ONE
1035 U2 = ZERO
1036 GO TO 250
1037 ELSE IF( TEMP.GE.TEMP2 ) THEN
1038 W11 = T( J+1, J+1 )
1039 W21 = T( J+2, J+1 )
1040 W12 = T( J+1, J+2 )
1041 W22 = T( J+2, J+2 )
1042 U1 = T( J+1, J )
1043 U2 = T( J+2, J )
1044 ELSE
1045 W21 = T( J+1, J+1 )
1046 W11 = T( J+2, J+1 )
1047 W22 = T( J+1, J+2 )
1048 W12 = T( J+2, J+2 )
1049 U2 = T( J+1, J )
1050 U1 = T( J+2, J )
1051 END IF
1052 *
1053 * Swap columns if nec.
1054 *
1055 IF( ABS( W12 ).GT.ABS( W11 ) ) THEN
1056 ILPIVT = .TRUE.
1057 TEMP = W12
1058 TEMP2 = W22
1059 W12 = W11
1060 W22 = W21
1061 W11 = TEMP
1062 W21 = TEMP2
1063 END IF
1064 *
1065 * LU-factor
1066 *
1067 TEMP = W21 / W11
1068 U2 = U2 - TEMP*U1
1069 W22 = W22 - TEMP*W12
1070 W21 = ZERO
1071 *
1072 * Compute SCALE
1073 *
1074 SCALE = ONE
1075 IF( ABS( W22 ).LT.SAFMIN ) THEN
1076 SCALE = ZERO
1077 U2 = ONE
1078 U1 = -W12 / W11
1079 GO TO 250
1080 END IF
1081 IF( ABS( W22 ).LT.ABS( U2 ) )
1082 $ SCALE = ABS( W22 / U2 )
1083 IF( ABS( W11 ).LT.ABS( U1 ) )
1084 $ SCALE = MIN( SCALE, ABS( W11 / U1 ) )
1085 *
1086 * Solve
1087 *
1088 U2 = ( SCALE*U2 ) / W22
1089 U1 = ( SCALE*U1-W12*U2 ) / W11
1090 *
1091 250 CONTINUE
1092 IF( ILPIVT ) THEN
1093 TEMP = U2
1094 U2 = U1
1095 U1 = TEMP
1096 END IF
1097 *
1098 * Compute Householder Vector
1099 *
1100 T1 = SQRT( SCALE**2+U1**2+U2**2 )
1101 TAU = ONE + SCALE / T1
1102 VS = -ONE / ( SCALE+T1 )
1103 V( 1 ) = ONE
1104 V( 2 ) = VS*U1
1105 V( 3 ) = VS*U2
1106 *
1107 * Apply transformations from the right.
1108 *
1109 DO 260 JR = IFRSTM, MIN( J+3, ILAST )
1110 TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )*
1111 $ H( JR, J+2 ) )
1112 H( JR, J ) = H( JR, J ) - TEMP
1113 H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 )
1114 H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 )
1115 260 CONTINUE
1116 DO 270 JR = IFRSTM, J + 2
1117 TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )*
1118 $ T( JR, J+2 ) )
1119 T( JR, J ) = T( JR, J ) - TEMP
1120 T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 )
1121 T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 )
1122 270 CONTINUE
1123 IF( ILZ ) THEN
1124 DO 280 JR = 1, N
1125 TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )*
1126 $ Z( JR, J+2 ) )
1127 Z( JR, J ) = Z( JR, J ) - TEMP
1128 Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 )
1129 Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 )
1130 280 CONTINUE
1131 END IF
1132 T( J+1, J ) = ZERO
1133 T( J+2, J ) = ZERO
1134 290 CONTINUE
1135 *
1136 * Last elements: Use Givens rotations
1137 *
1138 * Rotations from the left
1139 *
1140 J = ILAST - 1
1141 TEMP = H( J, J-1 )
1142 CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
1143 H( J+1, J-1 ) = ZERO
1144 *
1145 DO 300 JC = J, ILASTM
1146 TEMP = C*H( J, JC ) + S*H( J+1, JC )
1147 H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
1148 H( J, JC ) = TEMP
1149 TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
1150 T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
1151 T( J, JC ) = TEMP2
1152 300 CONTINUE
1153 IF( ILQ ) THEN
1154 DO 310 JR = 1, N
1155 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
1156 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
1157 Q( JR, J ) = TEMP
1158 310 CONTINUE
1159 END IF
1160 *
1161 * Rotations from the right.
1162 *
1163 TEMP = T( J+1, J+1 )
1164 CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
1165 T( J+1, J ) = ZERO
1166 *
1167 DO 320 JR = IFRSTM, ILAST
1168 TEMP = C*H( JR, J+1 ) + S*H( JR, J )
1169 H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
1170 H( JR, J+1 ) = TEMP
1171 320 CONTINUE
1172 DO 330 JR = IFRSTM, ILAST - 1
1173 TEMP = C*T( JR, J+1 ) + S*T( JR, J )
1174 T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
1175 T( JR, J+1 ) = TEMP
1176 330 CONTINUE
1177 IF( ILZ ) THEN
1178 DO 340 JR = 1, N
1179 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
1180 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
1181 Z( JR, J+1 ) = TEMP
1182 340 CONTINUE
1183 END IF
1184 *
1185 * End of Double-Shift code
1186 *
1187 END IF
1188 *
1189 GO TO 350
1190 *
1191 * End of iteration loop
1192 *
1193 350 CONTINUE
1194 360 CONTINUE
1195 *
1196 * Drop-through = non-convergence
1197 *
1198 INFO = ILAST
1199 GO TO 420
1200 *
1201 * Successful completion of all QZ steps
1202 *
1203 380 CONTINUE
1204 *
1205 * Set Eigenvalues 1:ILO-1
1206 *
1207 DO 410 J = 1, ILO - 1
1208 IF( T( J, J ).LT.ZERO ) THEN
1209 IF( ILSCHR ) THEN
1210 DO 390 JR = 1, J
1211 H( JR, J ) = -H( JR, J )
1212 T( JR, J ) = -T( JR, J )
1213 390 CONTINUE
1214 ELSE
1215 H( J, J ) = -H( J, J )
1216 T( J, J ) = -T( J, J )
1217 END IF
1218 IF( ILZ ) THEN
1219 DO 400 JR = 1, N
1220 Z( JR, J ) = -Z( JR, J )
1221 400 CONTINUE
1222 END IF
1223 END IF
1224 ALPHAR( J ) = H( J, J )
1225 ALPHAI( J ) = ZERO
1226 BETA( J ) = T( J, J )
1227 410 CONTINUE
1228 *
1229 * Normal Termination
1230 *
1231 INFO = 0
1232 *
1233 * Exit (other than argument error) -- return optimal workspace size
1234 *
1235 420 CONTINUE
1236 WORK( 1 ) = DBLE( N )
1237 RETURN
1238 *
1239 * End of DHGEQZ
1240 *
1241 END
2 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
3 $ LWORK, INFO )
4 *
5 * -- LAPACK routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2009 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER COMPQ, COMPZ, JOB
12 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ),
16 $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
17 $ WORK( * ), Z( LDZ, * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
24 * where H is an upper Hessenberg matrix and T is upper triangular,
25 * using the double-shift QZ method.
26 * Matrix pairs of this type are produced by the reduction to
27 * generalized upper Hessenberg form of a real matrix pair (A,B):
28 *
29 * A = Q1*H*Z1**T, B = Q1*T*Z1**T,
30 *
31 * as computed by DGGHRD.
32 *
33 * If JOB='S', then the Hessenberg-triangular pair (H,T) is
34 * also reduced to generalized Schur form,
35 *
36 * H = Q*S*Z**T, T = Q*P*Z**T,
37 *
38 * where Q and Z are orthogonal matrices, P is an upper triangular
39 * matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
40 * diagonal blocks.
41 *
42 * The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
43 * (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
44 * eigenvalues.
45 *
46 * Additionally, the 2-by-2 upper triangular diagonal blocks of P
47 * corresponding to 2-by-2 blocks of S are reduced to positive diagonal
48 * form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
49 * P(j,j) > 0, and P(j+1,j+1) > 0.
50 *
51 * Optionally, the orthogonal matrix Q from the generalized Schur
52 * factorization may be postmultiplied into an input matrix Q1, and the
53 * orthogonal matrix Z may be postmultiplied into an input matrix Z1.
54 * If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
55 * the matrix pair (A,B) to generalized upper Hessenberg form, then the
56 * output matrices Q1*Q and Z1*Z are the orthogonal factors from the
57 * generalized Schur factorization of (A,B):
58 *
59 * A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T.
60 *
61 * To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
62 * of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
63 * complex and beta real.
64 * If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
65 * generalized nonsymmetric eigenvalue problem (GNEP)
66 * A*x = lambda*B*x
67 * and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
68 * alternate form of the GNEP
69 * mu*A*y = B*y.
70 * Real eigenvalues can be read directly from the generalized Schur
71 * form:
72 * alpha = S(i,i), beta = P(i,i).
73 *
74 * Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
75 * Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
76 * pp. 241--256.
77 *
78 * Arguments
79 * =========
80 *
81 * JOB (input) CHARACTER*1
82 * = 'E': Compute eigenvalues only;
83 * = 'S': Compute eigenvalues and the Schur form.
84 *
85 * COMPQ (input) CHARACTER*1
86 * = 'N': Left Schur vectors (Q) are not computed;
87 * = 'I': Q is initialized to the unit matrix and the matrix Q
88 * of left Schur vectors of (H,T) is returned;
89 * = 'V': Q must contain an orthogonal matrix Q1 on entry and
90 * the product Q1*Q is returned.
91 *
92 * COMPZ (input) CHARACTER*1
93 * = 'N': Right Schur vectors (Z) are not computed;
94 * = 'I': Z is initialized to the unit matrix and the matrix Z
95 * of right Schur vectors of (H,T) is returned;
96 * = 'V': Z must contain an orthogonal matrix Z1 on entry and
97 * the product Z1*Z is returned.
98 *
99 * N (input) INTEGER
100 * The order of the matrices H, T, Q, and Z. N >= 0.
101 *
102 * ILO (input) INTEGER
103 * IHI (input) INTEGER
104 * ILO and IHI mark the rows and columns of H which are in
105 * Hessenberg form. It is assumed that A is already upper
106 * triangular in rows and columns 1:ILO-1 and IHI+1:N.
107 * If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
108 *
109 * H (input/output) DOUBLE PRECISION array, dimension (LDH, N)
110 * On entry, the N-by-N upper Hessenberg matrix H.
111 * On exit, if JOB = 'S', H contains the upper quasi-triangular
112 * matrix S from the generalized Schur factorization.
113 * If JOB = 'E', the diagonal blocks of H match those of S, but
114 * the rest of H is unspecified.
115 *
116 * LDH (input) INTEGER
117 * The leading dimension of the array H. LDH >= max( 1, N ).
118 *
119 * T (input/output) DOUBLE PRECISION array, dimension (LDT, N)
120 * On entry, the N-by-N upper triangular matrix T.
121 * On exit, if JOB = 'S', T contains the upper triangular
122 * matrix P from the generalized Schur factorization;
123 * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
124 * are reduced to positive diagonal form, i.e., if H(j+1,j) is
125 * non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
126 * T(j+1,j+1) > 0.
127 * If JOB = 'E', the diagonal blocks of T match those of P, but
128 * the rest of T is unspecified.
129 *
130 * LDT (input) INTEGER
131 * The leading dimension of the array T. LDT >= max( 1, N ).
132 *
133 * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
134 * The real parts of each scalar alpha defining an eigenvalue
135 * of GNEP.
136 *
137 * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
138 * The imaginary parts of each scalar alpha defining an
139 * eigenvalue of GNEP.
140 * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
141 * positive, then the j-th and (j+1)-st eigenvalues are a
142 * complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
143 *
144 * BETA (output) DOUBLE PRECISION array, dimension (N)
145 * The scalars beta that define the eigenvalues of GNEP.
146 * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
147 * beta = BETA(j) represent the j-th eigenvalue of the matrix
148 * pair (A,B), in one of the forms lambda = alpha/beta or
149 * mu = beta/alpha. Since either lambda or mu may overflow,
150 * they should not, in general, be computed.
151 *
152 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
153 * On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in
154 * the reduction of (A,B) to generalized Hessenberg form.
155 * On exit, if COMPZ = 'I', the orthogonal matrix of left Schur
156 * vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix
157 * of left Schur vectors of (A,B).
158 * Not referenced if COMPZ = 'N'.
159 *
160 * LDQ (input) INTEGER
161 * The leading dimension of the array Q. LDQ >= 1.
162 * If COMPQ='V' or 'I', then LDQ >= N.
163 *
164 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
165 * On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
166 * the reduction of (A,B) to generalized Hessenberg form.
167 * On exit, if COMPZ = 'I', the orthogonal matrix of
168 * right Schur vectors of (H,T), and if COMPZ = 'V', the
169 * orthogonal matrix of right Schur vectors of (A,B).
170 * Not referenced if COMPZ = 'N'.
171 *
172 * LDZ (input) INTEGER
173 * The leading dimension of the array Z. LDZ >= 1.
174 * If COMPZ='V' or 'I', then LDZ >= N.
175 *
176 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
177 * On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
178 *
179 * LWORK (input) INTEGER
180 * The dimension of the array WORK. LWORK >= max(1,N).
181 *
182 * If LWORK = -1, then a workspace query is assumed; the routine
183 * only calculates the optimal size of the WORK array, returns
184 * this value as the first entry of the WORK array, and no error
185 * message related to LWORK is issued by XERBLA.
186 *
187 * INFO (output) INTEGER
188 * = 0: successful exit
189 * < 0: if INFO = -i, the i-th argument had an illegal value
190 * = 1,...,N: the QZ iteration did not converge. (H,T) is not
191 * in Schur form, but ALPHAR(i), ALPHAI(i), and
192 * BETA(i), i=INFO+1,...,N should be correct.
193 * = N+1,...,2*N: the shift calculation failed. (H,T) is not
194 * in Schur form, but ALPHAR(i), ALPHAI(i), and
195 * BETA(i), i=INFO-N+1,...,N should be correct.
196 *
197 * Further Details
198 * ===============
199 *
200 * Iteration counters:
201 *
202 * JITER -- counts iterations.
203 * IITER -- counts iterations run since ILAST was last
204 * changed. This is therefore reset only when a 1-by-1 or
205 * 2-by-2 block deflates off the bottom.
206 *
207 * =====================================================================
208 *
209 * .. Parameters ..
210 * $ SAFETY = 1.0E+0 )
211 DOUBLE PRECISION HALF, ZERO, ONE, SAFETY
212 PARAMETER ( HALF = 0.5D+0, ZERO = 0.0D+0, ONE = 1.0D+0,
213 $ SAFETY = 1.0D+2 )
214 * ..
215 * .. Local Scalars ..
216 LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
217 $ LQUERY
218 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
219 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
220 $ JR, MAXIT
221 DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
222 $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
223 $ AD32L, AN, ANORM, ASCALE, ATOL, B11, B1A, B1I,
224 $ B1R, B22, B2A, B2I, B2R, BN, BNORM, BSCALE,
225 $ BTOL, C, C11I, C11R, C12, C21, C22I, C22R, CL,
226 $ CQ, CR, CZ, ESHIFT, S, S1, S1INV, S2, SAFMAX,
227 $ SAFMIN, SCALE, SL, SQI, SQR, SR, SZI, SZR, T1,
228 $ TAU, TEMP, TEMP2, TEMPI, TEMPR, U1, U12, U12L,
229 $ U2, ULP, VS, W11, W12, W21, W22, WABS, WI, WR,
230 $ WR2
231 * ..
232 * .. Local Arrays ..
233 DOUBLE PRECISION V( 3 )
234 * ..
235 * .. External Functions ..
236 LOGICAL LSAME
237 DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3
238 EXTERNAL LSAME, DLAMCH, DLANHS, DLAPY2, DLAPY3
239 * ..
240 * .. External Subroutines ..
241 EXTERNAL DLAG2, DLARFG, DLARTG, DLASET, DLASV2, DROT,
242 $ XERBLA
243 * ..
244 * .. Intrinsic Functions ..
245 INTRINSIC ABS, DBLE, MAX, MIN, SQRT
246 * ..
247 * .. Executable Statements ..
248 *
249 * Decode JOB, COMPQ, COMPZ
250 *
251 IF( LSAME( JOB, 'E' ) ) THEN
252 ILSCHR = .FALSE.
253 ISCHUR = 1
254 ELSE IF( LSAME( JOB, 'S' ) ) THEN
255 ILSCHR = .TRUE.
256 ISCHUR = 2
257 ELSE
258 ISCHUR = 0
259 END IF
260 *
261 IF( LSAME( COMPQ, 'N' ) ) THEN
262 ILQ = .FALSE.
263 ICOMPQ = 1
264 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
265 ILQ = .TRUE.
266 ICOMPQ = 2
267 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
268 ILQ = .TRUE.
269 ICOMPQ = 3
270 ELSE
271 ICOMPQ = 0
272 END IF
273 *
274 IF( LSAME( COMPZ, 'N' ) ) THEN
275 ILZ = .FALSE.
276 ICOMPZ = 1
277 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
278 ILZ = .TRUE.
279 ICOMPZ = 2
280 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
281 ILZ = .TRUE.
282 ICOMPZ = 3
283 ELSE
284 ICOMPZ = 0
285 END IF
286 *
287 * Check Argument Values
288 *
289 INFO = 0
290 WORK( 1 ) = MAX( 1, N )
291 LQUERY = ( LWORK.EQ.-1 )
292 IF( ISCHUR.EQ.0 ) THEN
293 INFO = -1
294 ELSE IF( ICOMPQ.EQ.0 ) THEN
295 INFO = -2
296 ELSE IF( ICOMPZ.EQ.0 ) THEN
297 INFO = -3
298 ELSE IF( N.LT.0 ) THEN
299 INFO = -4
300 ELSE IF( ILO.LT.1 ) THEN
301 INFO = -5
302 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
303 INFO = -6
304 ELSE IF( LDH.LT.N ) THEN
305 INFO = -8
306 ELSE IF( LDT.LT.N ) THEN
307 INFO = -10
308 ELSE IF( LDQ.LT.1 .OR. ( ILQ .AND. LDQ.LT.N ) ) THEN
309 INFO = -15
310 ELSE IF( LDZ.LT.1 .OR. ( ILZ .AND. LDZ.LT.N ) ) THEN
311 INFO = -17
312 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
313 INFO = -19
314 END IF
315 IF( INFO.NE.0 ) THEN
316 CALL XERBLA( 'DHGEQZ', -INFO )
317 RETURN
318 ELSE IF( LQUERY ) THEN
319 RETURN
320 END IF
321 *
322 * Quick return if possible
323 *
324 IF( N.LE.0 ) THEN
325 WORK( 1 ) = DBLE( 1 )
326 RETURN
327 END IF
328 *
329 * Initialize Q and Z
330 *
331 IF( ICOMPQ.EQ.3 )
332 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
333 IF( ICOMPZ.EQ.3 )
334 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
335 *
336 * Machine Constants
337 *
338 IN = IHI + 1 - ILO
339 SAFMIN = DLAMCH( 'S' )
340 SAFMAX = ONE / SAFMIN
341 ULP = DLAMCH( 'E' )*DLAMCH( 'B' )
342 ANORM = DLANHS( 'F', IN, H( ILO, ILO ), LDH, WORK )
343 BNORM = DLANHS( 'F', IN, T( ILO, ILO ), LDT, WORK )
344 ATOL = MAX( SAFMIN, ULP*ANORM )
345 BTOL = MAX( SAFMIN, ULP*BNORM )
346 ASCALE = ONE / MAX( SAFMIN, ANORM )
347 BSCALE = ONE / MAX( SAFMIN, BNORM )
348 *
349 * Set Eigenvalues IHI+1:N
350 *
351 DO 30 J = IHI + 1, N
352 IF( T( J, J ).LT.ZERO ) THEN
353 IF( ILSCHR ) THEN
354 DO 10 JR = 1, J
355 H( JR, J ) = -H( JR, J )
356 T( JR, J ) = -T( JR, J )
357 10 CONTINUE
358 ELSE
359 H( J, J ) = -H( J, J )
360 T( J, J ) = -T( J, J )
361 END IF
362 IF( ILZ ) THEN
363 DO 20 JR = 1, N
364 Z( JR, J ) = -Z( JR, J )
365 20 CONTINUE
366 END IF
367 END IF
368 ALPHAR( J ) = H( J, J )
369 ALPHAI( J ) = ZERO
370 BETA( J ) = T( J, J )
371 30 CONTINUE
372 *
373 * If IHI < ILO, skip QZ steps
374 *
375 IF( IHI.LT.ILO )
376 $ GO TO 380
377 *
378 * MAIN QZ ITERATION LOOP
379 *
380 * Initialize dynamic indices
381 *
382 * Eigenvalues ILAST+1:N have been found.
383 * Column operations modify rows IFRSTM:whatever.
384 * Row operations modify columns whatever:ILASTM.
385 *
386 * If only eigenvalues are being computed, then
387 * IFRSTM is the row of the last splitting row above row ILAST;
388 * this is always at least ILO.
389 * IITER counts iterations since the last eigenvalue was found,
390 * to tell when to use an extraordinary shift.
391 * MAXIT is the maximum number of QZ sweeps allowed.
392 *
393 ILAST = IHI
394 IF( ILSCHR ) THEN
395 IFRSTM = 1
396 ILASTM = N
397 ELSE
398 IFRSTM = ILO
399 ILASTM = IHI
400 END IF
401 IITER = 0
402 ESHIFT = ZERO
403 MAXIT = 30*( IHI-ILO+1 )
404 *
405 DO 360 JITER = 1, MAXIT
406 *
407 * Split the matrix if possible.
408 *
409 * Two tests:
410 * 1: H(j,j-1)=0 or j=ILO
411 * 2: T(j,j)=0
412 *
413 IF( ILAST.EQ.ILO ) THEN
414 *
415 * Special case: j=ILAST
416 *
417 GO TO 80
418 ELSE
419 IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN
420 H( ILAST, ILAST-1 ) = ZERO
421 GO TO 80
422 END IF
423 END IF
424 *
425 IF( ABS( T( ILAST, ILAST ) ).LE.BTOL ) THEN
426 T( ILAST, ILAST ) = ZERO
427 GO TO 70
428 END IF
429 *
430 * General case: j<ILAST
431 *
432 DO 60 J = ILAST - 1, ILO, -1
433 *
434 * Test 1: for H(j,j-1)=0 or j=ILO
435 *
436 IF( J.EQ.ILO ) THEN
437 ILAZRO = .TRUE.
438 ELSE
439 IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN
440 H( J, J-1 ) = ZERO
441 ILAZRO = .TRUE.
442 ELSE
443 ILAZRO = .FALSE.
444 END IF
445 END IF
446 *
447 * Test 2: for T(j,j)=0
448 *
449 IF( ABS( T( J, J ) ).LT.BTOL ) THEN
450 T( J, J ) = ZERO
451 *
452 * Test 1a: Check for 2 consecutive small subdiagonals in A
453 *
454 ILAZR2 = .FALSE.
455 IF( .NOT.ILAZRO ) THEN
456 TEMP = ABS( H( J, J-1 ) )
457 TEMP2 = ABS( H( J, J ) )
458 TEMPR = MAX( TEMP, TEMP2 )
459 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
460 TEMP = TEMP / TEMPR
461 TEMP2 = TEMP2 / TEMPR
462 END IF
463 IF( TEMP*( ASCALE*ABS( H( J+1, J ) ) ).LE.TEMP2*
464 $ ( ASCALE*ATOL ) )ILAZR2 = .TRUE.
465 END IF
466 *
467 * If both tests pass (1 & 2), i.e., the leading diagonal
468 * element of B in the block is zero, split a 1x1 block off
469 * at the top. (I.e., at the J-th row/column) The leading
470 * diagonal element of the remainder can also be zero, so
471 * this may have to be done repeatedly.
472 *
473 IF( ILAZRO .OR. ILAZR2 ) THEN
474 DO 40 JCH = J, ILAST - 1
475 TEMP = H( JCH, JCH )
476 CALL DLARTG( TEMP, H( JCH+1, JCH ), C, S,
477 $ H( JCH, JCH ) )
478 H( JCH+1, JCH ) = ZERO
479 CALL DROT( ILASTM-JCH, H( JCH, JCH+1 ), LDH,
480 $ H( JCH+1, JCH+1 ), LDH, C, S )
481 CALL DROT( ILASTM-JCH, T( JCH, JCH+1 ), LDT,
482 $ T( JCH+1, JCH+1 ), LDT, C, S )
483 IF( ILQ )
484 $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
485 $ C, S )
486 IF( ILAZR2 )
487 $ H( JCH, JCH-1 ) = H( JCH, JCH-1 )*C
488 ILAZR2 = .FALSE.
489 IF( ABS( T( JCH+1, JCH+1 ) ).GE.BTOL ) THEN
490 IF( JCH+1.GE.ILAST ) THEN
491 GO TO 80
492 ELSE
493 IFIRST = JCH + 1
494 GO TO 110
495 END IF
496 END IF
497 T( JCH+1, JCH+1 ) = ZERO
498 40 CONTINUE
499 GO TO 70
500 ELSE
501 *
502 * Only test 2 passed -- chase the zero to T(ILAST,ILAST)
503 * Then process as in the case T(ILAST,ILAST)=0
504 *
505 DO 50 JCH = J, ILAST - 1
506 TEMP = T( JCH, JCH+1 )
507 CALL DLARTG( TEMP, T( JCH+1, JCH+1 ), C, S,
508 $ T( JCH, JCH+1 ) )
509 T( JCH+1, JCH+1 ) = ZERO
510 IF( JCH.LT.ILASTM-1 )
511 $ CALL DROT( ILASTM-JCH-1, T( JCH, JCH+2 ), LDT,
512 $ T( JCH+1, JCH+2 ), LDT, C, S )
513 CALL DROT( ILASTM-JCH+2, H( JCH, JCH-1 ), LDH,
514 $ H( JCH+1, JCH-1 ), LDH, C, S )
515 IF( ILQ )
516 $ CALL DROT( N, Q( 1, JCH ), 1, Q( 1, JCH+1 ), 1,
517 $ C, S )
518 TEMP = H( JCH+1, JCH )
519 CALL DLARTG( TEMP, H( JCH+1, JCH-1 ), C, S,
520 $ H( JCH+1, JCH ) )
521 H( JCH+1, JCH-1 ) = ZERO
522 CALL DROT( JCH+1-IFRSTM, H( IFRSTM, JCH ), 1,
523 $ H( IFRSTM, JCH-1 ), 1, C, S )
524 CALL DROT( JCH-IFRSTM, T( IFRSTM, JCH ), 1,
525 $ T( IFRSTM, JCH-1 ), 1, C, S )
526 IF( ILZ )
527 $ CALL DROT( N, Z( 1, JCH ), 1, Z( 1, JCH-1 ), 1,
528 $ C, S )
529 50 CONTINUE
530 GO TO 70
531 END IF
532 ELSE IF( ILAZRO ) THEN
533 *
534 * Only test 1 passed -- work on J:ILAST
535 *
536 IFIRST = J
537 GO TO 110
538 END IF
539 *
540 * Neither test passed -- try next J
541 *
542 60 CONTINUE
543 *
544 * (Drop-through is "impossible")
545 *
546 INFO = N + 1
547 GO TO 420
548 *
549 * T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
550 * 1x1 block.
551 *
552 70 CONTINUE
553 TEMP = H( ILAST, ILAST )
554 CALL DLARTG( TEMP, H( ILAST, ILAST-1 ), C, S,
555 $ H( ILAST, ILAST ) )
556 H( ILAST, ILAST-1 ) = ZERO
557 CALL DROT( ILAST-IFRSTM, H( IFRSTM, ILAST ), 1,
558 $ H( IFRSTM, ILAST-1 ), 1, C, S )
559 CALL DROT( ILAST-IFRSTM, T( IFRSTM, ILAST ), 1,
560 $ T( IFRSTM, ILAST-1 ), 1, C, S )
561 IF( ILZ )
562 $ CALL DROT( N, Z( 1, ILAST ), 1, Z( 1, ILAST-1 ), 1, C, S )
563 *
564 * H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
565 * and BETA
566 *
567 80 CONTINUE
568 IF( T( ILAST, ILAST ).LT.ZERO ) THEN
569 IF( ILSCHR ) THEN
570 DO 90 J = IFRSTM, ILAST
571 H( J, ILAST ) = -H( J, ILAST )
572 T( J, ILAST ) = -T( J, ILAST )
573 90 CONTINUE
574 ELSE
575 H( ILAST, ILAST ) = -H( ILAST, ILAST )
576 T( ILAST, ILAST ) = -T( ILAST, ILAST )
577 END IF
578 IF( ILZ ) THEN
579 DO 100 J = 1, N
580 Z( J, ILAST ) = -Z( J, ILAST )
581 100 CONTINUE
582 END IF
583 END IF
584 ALPHAR( ILAST ) = H( ILAST, ILAST )
585 ALPHAI( ILAST ) = ZERO
586 BETA( ILAST ) = T( ILAST, ILAST )
587 *
588 * Go to next block -- exit if finished.
589 *
590 ILAST = ILAST - 1
591 IF( ILAST.LT.ILO )
592 $ GO TO 380
593 *
594 * Reset counters
595 *
596 IITER = 0
597 ESHIFT = ZERO
598 IF( .NOT.ILSCHR ) THEN
599 ILASTM = ILAST
600 IF( IFRSTM.GT.ILAST )
601 $ IFRSTM = ILO
602 END IF
603 GO TO 350
604 *
605 * QZ step
606 *
607 * This iteration only involves rows/columns IFIRST:ILAST. We
608 * assume IFIRST < ILAST, and that the diagonal of B is non-zero.
609 *
610 110 CONTINUE
611 IITER = IITER + 1
612 IF( .NOT.ILSCHR ) THEN
613 IFRSTM = IFIRST
614 END IF
615 *
616 * Compute single shifts.
617 *
618 * At this point, IFIRST < ILAST, and the diagonal elements of
619 * T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
620 * magnitude)
621 *
622 IF( ( IITER / 10 )*10.EQ.IITER ) THEN
623 *
624 * Exceptional shift. Chosen for no particularly good reason.
625 * (Single shift only.)
626 *
627 IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.
628 $ ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
629 ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /
630 $ T( ILAST-1, ILAST-1 )
631 ELSE
632 ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
633 END IF
634 S1 = ONE
635 WR = ESHIFT
636 *
637 ELSE
638 *
639 * Shifts based on the generalized eigenvalues of the
640 * bottom-right 2x2 block of A and B. The first eigenvalue
641 * returned by DLAG2 is the Wilkinson shift (AEP p.512),
642 *
643 CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
644 $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
645 $ S2, WR, WR2, WI )
646 *
647 TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
648 IF( WI.NE.ZERO )
649 $ GO TO 200
650 END IF
651 *
652 * Fiddle with shift to avoid overflow
653 *
654 TEMP = MIN( ASCALE, ONE )*( HALF*SAFMAX )
655 IF( S1.GT.TEMP ) THEN
656 SCALE = TEMP / S1
657 ELSE
658 SCALE = ONE
659 END IF
660 *
661 TEMP = MIN( BSCALE, ONE )*( HALF*SAFMAX )
662 IF( ABS( WR ).GT.TEMP )
663 $ SCALE = MIN( SCALE, TEMP / ABS( WR ) )
664 S1 = SCALE*S1
665 WR = SCALE*WR
666 *
667 * Now check for two consecutive small subdiagonals.
668 *
669 DO 120 J = ILAST - 1, IFIRST + 1, -1
670 ISTART = J
671 TEMP = ABS( S1*H( J, J-1 ) )
672 TEMP2 = ABS( S1*H( J, J )-WR*T( J, J ) )
673 TEMPR = MAX( TEMP, TEMP2 )
674 IF( TEMPR.LT.ONE .AND. TEMPR.NE.ZERO ) THEN
675 TEMP = TEMP / TEMPR
676 TEMP2 = TEMP2 / TEMPR
677 END IF
678 IF( ABS( ( ASCALE*H( J+1, J ) )*TEMP ).LE.( ASCALE*ATOL )*
679 $ TEMP2 )GO TO 130
680 120 CONTINUE
681 *
682 ISTART = IFIRST
683 130 CONTINUE
684 *
685 * Do an implicit single-shift QZ sweep.
686 *
687 * Initial Q
688 *
689 TEMP = S1*H( ISTART, ISTART ) - WR*T( ISTART, ISTART )
690 TEMP2 = S1*H( ISTART+1, ISTART )
691 CALL DLARTG( TEMP, TEMP2, C, S, TEMPR )
692 *
693 * Sweep
694 *
695 DO 190 J = ISTART, ILAST - 1
696 IF( J.GT.ISTART ) THEN
697 TEMP = H( J, J-1 )
698 CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
699 H( J+1, J-1 ) = ZERO
700 END IF
701 *
702 DO 140 JC = J, ILASTM
703 TEMP = C*H( J, JC ) + S*H( J+1, JC )
704 H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
705 H( J, JC ) = TEMP
706 TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
707 T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
708 T( J, JC ) = TEMP2
709 140 CONTINUE
710 IF( ILQ ) THEN
711 DO 150 JR = 1, N
712 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
713 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
714 Q( JR, J ) = TEMP
715 150 CONTINUE
716 END IF
717 *
718 TEMP = T( J+1, J+1 )
719 CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
720 T( J+1, J ) = ZERO
721 *
722 DO 160 JR = IFRSTM, MIN( J+2, ILAST )
723 TEMP = C*H( JR, J+1 ) + S*H( JR, J )
724 H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
725 H( JR, J+1 ) = TEMP
726 160 CONTINUE
727 DO 170 JR = IFRSTM, J
728 TEMP = C*T( JR, J+1 ) + S*T( JR, J )
729 T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
730 T( JR, J+1 ) = TEMP
731 170 CONTINUE
732 IF( ILZ ) THEN
733 DO 180 JR = 1, N
734 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
735 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
736 Z( JR, J+1 ) = TEMP
737 180 CONTINUE
738 END IF
739 190 CONTINUE
740 *
741 GO TO 350
742 *
743 * Use Francis double-shift
744 *
745 * Note: the Francis double-shift should work with real shifts,
746 * but only if the block is at least 3x3.
747 * This code may break if this point is reached with
748 * a 2x2 block with real eigenvalues.
749 *
750 200 CONTINUE
751 IF( IFIRST+1.EQ.ILAST ) THEN
752 *
753 * Special case -- 2x2 block with complex eigenvectors
754 *
755 * Step 1: Standardize, that is, rotate so that
756 *
757 * ( B11 0 )
758 * B = ( ) with B11 non-negative.
759 * ( 0 B22 )
760 *
761 CALL DLASV2( T( ILAST-1, ILAST-1 ), T( ILAST-1, ILAST ),
762 $ T( ILAST, ILAST ), B22, B11, SR, CR, SL, CL )
763 *
764 IF( B11.LT.ZERO ) THEN
765 CR = -CR
766 SR = -SR
767 B11 = -B11
768 B22 = -B22
769 END IF
770 *
771 CALL DROT( ILASTM+1-IFIRST, H( ILAST-1, ILAST-1 ), LDH,
772 $ H( ILAST, ILAST-1 ), LDH, CL, SL )
773 CALL DROT( ILAST+1-IFRSTM, H( IFRSTM, ILAST-1 ), 1,
774 $ H( IFRSTM, ILAST ), 1, CR, SR )
775 *
776 IF( ILAST.LT.ILASTM )
777 $ CALL DROT( ILASTM-ILAST, T( ILAST-1, ILAST+1 ), LDT,
778 $ T( ILAST, ILAST+1 ), LDT, CL, SL )
779 IF( IFRSTM.LT.ILAST-1 )
780 $ CALL DROT( IFIRST-IFRSTM, T( IFRSTM, ILAST-1 ), 1,
781 $ T( IFRSTM, ILAST ), 1, CR, SR )
782 *
783 IF( ILQ )
784 $ CALL DROT( N, Q( 1, ILAST-1 ), 1, Q( 1, ILAST ), 1, CL,
785 $ SL )
786 IF( ILZ )
787 $ CALL DROT( N, Z( 1, ILAST-1 ), 1, Z( 1, ILAST ), 1, CR,
788 $ SR )
789 *
790 T( ILAST-1, ILAST-1 ) = B11
791 T( ILAST-1, ILAST ) = ZERO
792 T( ILAST, ILAST-1 ) = ZERO
793 T( ILAST, ILAST ) = B22
794 *
795 * If B22 is negative, negate column ILAST
796 *
797 IF( B22.LT.ZERO ) THEN
798 DO 210 J = IFRSTM, ILAST
799 H( J, ILAST ) = -H( J, ILAST )
800 T( J, ILAST ) = -T( J, ILAST )
801 210 CONTINUE
802 *
803 IF( ILZ ) THEN
804 DO 220 J = 1, N
805 Z( J, ILAST ) = -Z( J, ILAST )
806 220 CONTINUE
807 END IF
808 END IF
809 *
810 * Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
811 *
812 * Recompute shift
813 *
814 CALL DLAG2( H( ILAST-1, ILAST-1 ), LDH,
815 $ T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
816 $ TEMP, WR, TEMP2, WI )
817 *
818 * If standardization has perturbed the shift onto real line,
819 * do another (real single-shift) QR step.
820 *
821 IF( WI.EQ.ZERO )
822 $ GO TO 350
823 S1INV = ONE / S1
824 *
825 * Do EISPACK (QZVAL) computation of alpha and beta
826 *
827 A11 = H( ILAST-1, ILAST-1 )
828 A21 = H( ILAST, ILAST-1 )
829 A12 = H( ILAST-1, ILAST )
830 A22 = H( ILAST, ILAST )
831 *
832 * Compute complex Givens rotation on right
833 * (Assume some element of C = (sA - wB) > unfl )
834 * __
835 * (sA - wB) ( CZ -SZ )
836 * ( SZ CZ )
837 *
838 C11R = S1*A11 - WR*B11
839 C11I = -WI*B11
840 C12 = S1*A12
841 C21 = S1*A21
842 C22R = S1*A22 - WR*B22
843 C22I = -WI*B22
844 *
845 IF( ABS( C11R )+ABS( C11I )+ABS( C12 ).GT.ABS( C21 )+
846 $ ABS( C22R )+ABS( C22I ) ) THEN
847 T1 = DLAPY3( C12, C11R, C11I )
848 CZ = C12 / T1
849 SZR = -C11R / T1
850 SZI = -C11I / T1
851 ELSE
852 CZ = DLAPY2( C22R, C22I )
853 IF( CZ.LE.SAFMIN ) THEN
854 CZ = ZERO
855 SZR = ONE
856 SZI = ZERO
857 ELSE
858 TEMPR = C22R / CZ
859 TEMPI = C22I / CZ
860 T1 = DLAPY2( CZ, C21 )
861 CZ = CZ / T1
862 SZR = -C21*TEMPR / T1
863 SZI = C21*TEMPI / T1
864 END IF
865 END IF
866 *
867 * Compute Givens rotation on left
868 *
869 * ( CQ SQ )
870 * ( __ ) A or B
871 * ( -SQ CQ )
872 *
873 AN = ABS( A11 ) + ABS( A12 ) + ABS( A21 ) + ABS( A22 )
874 BN = ABS( B11 ) + ABS( B22 )
875 WABS = ABS( WR ) + ABS( WI )
876 IF( S1*AN.GT.WABS*BN ) THEN
877 CQ = CZ*B11
878 SQR = SZR*B22
879 SQI = -SZI*B22
880 ELSE
881 A1R = CZ*A11 + SZR*A12
882 A1I = SZI*A12
883 A2R = CZ*A21 + SZR*A22
884 A2I = SZI*A22
885 CQ = DLAPY2( A1R, A1I )
886 IF( CQ.LE.SAFMIN ) THEN
887 CQ = ZERO
888 SQR = ONE
889 SQI = ZERO
890 ELSE
891 TEMPR = A1R / CQ
892 TEMPI = A1I / CQ
893 SQR = TEMPR*A2R + TEMPI*A2I
894 SQI = TEMPI*A2R - TEMPR*A2I
895 END IF
896 END IF
897 T1 = DLAPY3( CQ, SQR, SQI )
898 CQ = CQ / T1
899 SQR = SQR / T1
900 SQI = SQI / T1
901 *
902 * Compute diagonal elements of QBZ
903 *
904 TEMPR = SQR*SZR - SQI*SZI
905 TEMPI = SQR*SZI + SQI*SZR
906 B1R = CQ*CZ*B11 + TEMPR*B22
907 B1I = TEMPI*B22
908 B1A = DLAPY2( B1R, B1I )
909 B2R = CQ*CZ*B22 + TEMPR*B11
910 B2I = -TEMPI*B11
911 B2A = DLAPY2( B2R, B2I )
912 *
913 * Normalize so beta > 0, and Im( alpha1 ) > 0
914 *
915 BETA( ILAST-1 ) = B1A
916 BETA( ILAST ) = B2A
917 ALPHAR( ILAST-1 ) = ( WR*B1A )*S1INV
918 ALPHAI( ILAST-1 ) = ( WI*B1A )*S1INV
919 ALPHAR( ILAST ) = ( WR*B2A )*S1INV
920 ALPHAI( ILAST ) = -( WI*B2A )*S1INV
921 *
922 * Step 3: Go to next block -- exit if finished.
923 *
924 ILAST = IFIRST - 1
925 IF( ILAST.LT.ILO )
926 $ GO TO 380
927 *
928 * Reset counters
929 *
930 IITER = 0
931 ESHIFT = ZERO
932 IF( .NOT.ILSCHR ) THEN
933 ILASTM = ILAST
934 IF( IFRSTM.GT.ILAST )
935 $ IFRSTM = ILO
936 END IF
937 GO TO 350
938 ELSE
939 *
940 * Usual case: 3x3 or larger block, using Francis implicit
941 * double-shift
942 *
943 * 2
944 * Eigenvalue equation is w - c w + d = 0,
945 *
946 * -1 2 -1
947 * so compute 1st column of (A B ) - c A B + d
948 * using the formula in QZIT (from EISPACK)
949 *
950 * We assume that the block is at least 3x3
951 *
952 AD11 = ( ASCALE*H( ILAST-1, ILAST-1 ) ) /
953 $ ( BSCALE*T( ILAST-1, ILAST-1 ) )
954 AD21 = ( ASCALE*H( ILAST, ILAST-1 ) ) /
955 $ ( BSCALE*T( ILAST-1, ILAST-1 ) )
956 AD12 = ( ASCALE*H( ILAST-1, ILAST ) ) /
957 $ ( BSCALE*T( ILAST, ILAST ) )
958 AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
959 $ ( BSCALE*T( ILAST, ILAST ) )
960 U12 = T( ILAST-1, ILAST ) / T( ILAST, ILAST )
961 AD11L = ( ASCALE*H( IFIRST, IFIRST ) ) /
962 $ ( BSCALE*T( IFIRST, IFIRST ) )
963 AD21L = ( ASCALE*H( IFIRST+1, IFIRST ) ) /
964 $ ( BSCALE*T( IFIRST, IFIRST ) )
965 AD12L = ( ASCALE*H( IFIRST, IFIRST+1 ) ) /
966 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
967 AD22L = ( ASCALE*H( IFIRST+1, IFIRST+1 ) ) /
968 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
969 AD32L = ( ASCALE*H( IFIRST+2, IFIRST+1 ) ) /
970 $ ( BSCALE*T( IFIRST+1, IFIRST+1 ) )
971 U12L = T( IFIRST, IFIRST+1 ) / T( IFIRST+1, IFIRST+1 )
972 *
973 V( 1 ) = ( AD11-AD11L )*( AD22-AD11L ) - AD12*AD21 +
974 $ AD21*U12*AD11L + ( AD12L-AD11L*U12L )*AD21L
975 V( 2 ) = ( ( AD22L-AD11L )-AD21L*U12L-( AD11-AD11L )-
976 $ ( AD22-AD11L )+AD21*U12 )*AD21L
977 V( 3 ) = AD32L*AD21L
978 *
979 ISTART = IFIRST
980 *
981 CALL DLARFG( 3, V( 1 ), V( 2 ), 1, TAU )
982 V( 1 ) = ONE
983 *
984 * Sweep
985 *
986 DO 290 J = ISTART, ILAST - 2
987 *
988 * All but last elements: use 3x3 Householder transforms.
989 *
990 * Zero (j-1)st column of A
991 *
992 IF( J.GT.ISTART ) THEN
993 V( 1 ) = H( J, J-1 )
994 V( 2 ) = H( J+1, J-1 )
995 V( 3 ) = H( J+2, J-1 )
996 *
997 CALL DLARFG( 3, H( J, J-1 ), V( 2 ), 1, TAU )
998 V( 1 ) = ONE
999 H( J+1, J-1 ) = ZERO
1000 H( J+2, J-1 ) = ZERO
1001 END IF
1002 *
1003 DO 230 JC = J, ILASTM
1004 TEMP = TAU*( H( J, JC )+V( 2 )*H( J+1, JC )+V( 3 )*
1005 $ H( J+2, JC ) )
1006 H( J, JC ) = H( J, JC ) - TEMP
1007 H( J+1, JC ) = H( J+1, JC ) - TEMP*V( 2 )
1008 H( J+2, JC ) = H( J+2, JC ) - TEMP*V( 3 )
1009 TEMP2 = TAU*( T( J, JC )+V( 2 )*T( J+1, JC )+V( 3 )*
1010 $ T( J+2, JC ) )
1011 T( J, JC ) = T( J, JC ) - TEMP2
1012 T( J+1, JC ) = T( J+1, JC ) - TEMP2*V( 2 )
1013 T( J+2, JC ) = T( J+2, JC ) - TEMP2*V( 3 )
1014 230 CONTINUE
1015 IF( ILQ ) THEN
1016 DO 240 JR = 1, N
1017 TEMP = TAU*( Q( JR, J )+V( 2 )*Q( JR, J+1 )+V( 3 )*
1018 $ Q( JR, J+2 ) )
1019 Q( JR, J ) = Q( JR, J ) - TEMP
1020 Q( JR, J+1 ) = Q( JR, J+1 ) - TEMP*V( 2 )
1021 Q( JR, J+2 ) = Q( JR, J+2 ) - TEMP*V( 3 )
1022 240 CONTINUE
1023 END IF
1024 *
1025 * Zero j-th column of B (see DLAGBC for details)
1026 *
1027 * Swap rows to pivot
1028 *
1029 ILPIVT = .FALSE.
1030 TEMP = MAX( ABS( T( J+1, J+1 ) ), ABS( T( J+1, J+2 ) ) )
1031 TEMP2 = MAX( ABS( T( J+2, J+1 ) ), ABS( T( J+2, J+2 ) ) )
1032 IF( MAX( TEMP, TEMP2 ).LT.SAFMIN ) THEN
1033 SCALE = ZERO
1034 U1 = ONE
1035 U2 = ZERO
1036 GO TO 250
1037 ELSE IF( TEMP.GE.TEMP2 ) THEN
1038 W11 = T( J+1, J+1 )
1039 W21 = T( J+2, J+1 )
1040 W12 = T( J+1, J+2 )
1041 W22 = T( J+2, J+2 )
1042 U1 = T( J+1, J )
1043 U2 = T( J+2, J )
1044 ELSE
1045 W21 = T( J+1, J+1 )
1046 W11 = T( J+2, J+1 )
1047 W22 = T( J+1, J+2 )
1048 W12 = T( J+2, J+2 )
1049 U2 = T( J+1, J )
1050 U1 = T( J+2, J )
1051 END IF
1052 *
1053 * Swap columns if nec.
1054 *
1055 IF( ABS( W12 ).GT.ABS( W11 ) ) THEN
1056 ILPIVT = .TRUE.
1057 TEMP = W12
1058 TEMP2 = W22
1059 W12 = W11
1060 W22 = W21
1061 W11 = TEMP
1062 W21 = TEMP2
1063 END IF
1064 *
1065 * LU-factor
1066 *
1067 TEMP = W21 / W11
1068 U2 = U2 - TEMP*U1
1069 W22 = W22 - TEMP*W12
1070 W21 = ZERO
1071 *
1072 * Compute SCALE
1073 *
1074 SCALE = ONE
1075 IF( ABS( W22 ).LT.SAFMIN ) THEN
1076 SCALE = ZERO
1077 U2 = ONE
1078 U1 = -W12 / W11
1079 GO TO 250
1080 END IF
1081 IF( ABS( W22 ).LT.ABS( U2 ) )
1082 $ SCALE = ABS( W22 / U2 )
1083 IF( ABS( W11 ).LT.ABS( U1 ) )
1084 $ SCALE = MIN( SCALE, ABS( W11 / U1 ) )
1085 *
1086 * Solve
1087 *
1088 U2 = ( SCALE*U2 ) / W22
1089 U1 = ( SCALE*U1-W12*U2 ) / W11
1090 *
1091 250 CONTINUE
1092 IF( ILPIVT ) THEN
1093 TEMP = U2
1094 U2 = U1
1095 U1 = TEMP
1096 END IF
1097 *
1098 * Compute Householder Vector
1099 *
1100 T1 = SQRT( SCALE**2+U1**2+U2**2 )
1101 TAU = ONE + SCALE / T1
1102 VS = -ONE / ( SCALE+T1 )
1103 V( 1 ) = ONE
1104 V( 2 ) = VS*U1
1105 V( 3 ) = VS*U2
1106 *
1107 * Apply transformations from the right.
1108 *
1109 DO 260 JR = IFRSTM, MIN( J+3, ILAST )
1110 TEMP = TAU*( H( JR, J )+V( 2 )*H( JR, J+1 )+V( 3 )*
1111 $ H( JR, J+2 ) )
1112 H( JR, J ) = H( JR, J ) - TEMP
1113 H( JR, J+1 ) = H( JR, J+1 ) - TEMP*V( 2 )
1114 H( JR, J+2 ) = H( JR, J+2 ) - TEMP*V( 3 )
1115 260 CONTINUE
1116 DO 270 JR = IFRSTM, J + 2
1117 TEMP = TAU*( T( JR, J )+V( 2 )*T( JR, J+1 )+V( 3 )*
1118 $ T( JR, J+2 ) )
1119 T( JR, J ) = T( JR, J ) - TEMP
1120 T( JR, J+1 ) = T( JR, J+1 ) - TEMP*V( 2 )
1121 T( JR, J+2 ) = T( JR, J+2 ) - TEMP*V( 3 )
1122 270 CONTINUE
1123 IF( ILZ ) THEN
1124 DO 280 JR = 1, N
1125 TEMP = TAU*( Z( JR, J )+V( 2 )*Z( JR, J+1 )+V( 3 )*
1126 $ Z( JR, J+2 ) )
1127 Z( JR, J ) = Z( JR, J ) - TEMP
1128 Z( JR, J+1 ) = Z( JR, J+1 ) - TEMP*V( 2 )
1129 Z( JR, J+2 ) = Z( JR, J+2 ) - TEMP*V( 3 )
1130 280 CONTINUE
1131 END IF
1132 T( J+1, J ) = ZERO
1133 T( J+2, J ) = ZERO
1134 290 CONTINUE
1135 *
1136 * Last elements: Use Givens rotations
1137 *
1138 * Rotations from the left
1139 *
1140 J = ILAST - 1
1141 TEMP = H( J, J-1 )
1142 CALL DLARTG( TEMP, H( J+1, J-1 ), C, S, H( J, J-1 ) )
1143 H( J+1, J-1 ) = ZERO
1144 *
1145 DO 300 JC = J, ILASTM
1146 TEMP = C*H( J, JC ) + S*H( J+1, JC )
1147 H( J+1, JC ) = -S*H( J, JC ) + C*H( J+1, JC )
1148 H( J, JC ) = TEMP
1149 TEMP2 = C*T( J, JC ) + S*T( J+1, JC )
1150 T( J+1, JC ) = -S*T( J, JC ) + C*T( J+1, JC )
1151 T( J, JC ) = TEMP2
1152 300 CONTINUE
1153 IF( ILQ ) THEN
1154 DO 310 JR = 1, N
1155 TEMP = C*Q( JR, J ) + S*Q( JR, J+1 )
1156 Q( JR, J+1 ) = -S*Q( JR, J ) + C*Q( JR, J+1 )
1157 Q( JR, J ) = TEMP
1158 310 CONTINUE
1159 END IF
1160 *
1161 * Rotations from the right.
1162 *
1163 TEMP = T( J+1, J+1 )
1164 CALL DLARTG( TEMP, T( J+1, J ), C, S, T( J+1, J+1 ) )
1165 T( J+1, J ) = ZERO
1166 *
1167 DO 320 JR = IFRSTM, ILAST
1168 TEMP = C*H( JR, J+1 ) + S*H( JR, J )
1169 H( JR, J ) = -S*H( JR, J+1 ) + C*H( JR, J )
1170 H( JR, J+1 ) = TEMP
1171 320 CONTINUE
1172 DO 330 JR = IFRSTM, ILAST - 1
1173 TEMP = C*T( JR, J+1 ) + S*T( JR, J )
1174 T( JR, J ) = -S*T( JR, J+1 ) + C*T( JR, J )
1175 T( JR, J+1 ) = TEMP
1176 330 CONTINUE
1177 IF( ILZ ) THEN
1178 DO 340 JR = 1, N
1179 TEMP = C*Z( JR, J+1 ) + S*Z( JR, J )
1180 Z( JR, J ) = -S*Z( JR, J+1 ) + C*Z( JR, J )
1181 Z( JR, J+1 ) = TEMP
1182 340 CONTINUE
1183 END IF
1184 *
1185 * End of Double-Shift code
1186 *
1187 END IF
1188 *
1189 GO TO 350
1190 *
1191 * End of iteration loop
1192 *
1193 350 CONTINUE
1194 360 CONTINUE
1195 *
1196 * Drop-through = non-convergence
1197 *
1198 INFO = ILAST
1199 GO TO 420
1200 *
1201 * Successful completion of all QZ steps
1202 *
1203 380 CONTINUE
1204 *
1205 * Set Eigenvalues 1:ILO-1
1206 *
1207 DO 410 J = 1, ILO - 1
1208 IF( T( J, J ).LT.ZERO ) THEN
1209 IF( ILSCHR ) THEN
1210 DO 390 JR = 1, J
1211 H( JR, J ) = -H( JR, J )
1212 T( JR, J ) = -T( JR, J )
1213 390 CONTINUE
1214 ELSE
1215 H( J, J ) = -H( J, J )
1216 T( J, J ) = -T( J, J )
1217 END IF
1218 IF( ILZ ) THEN
1219 DO 400 JR = 1, N
1220 Z( JR, J ) = -Z( JR, J )
1221 400 CONTINUE
1222 END IF
1223 END IF
1224 ALPHAR( J ) = H( J, J )
1225 ALPHAI( J ) = ZERO
1226 BETA( J ) = T( J, J )
1227 410 CONTINUE
1228 *
1229 * Normal Termination
1230 *
1231 INFO = 0
1232 *
1233 * Exit (other than argument error) -- return optimal workspace size
1234 *
1235 420 CONTINUE
1236 WORK( 1 ) = DBLE( N )
1237 RETURN
1238 *
1239 * End of DHGEQZ
1240 *
1241 END