![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dhgeqz.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dhgeqz.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 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: *
1.8 ! bertrand 5: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 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
1.8 ! bertrand 112: * matrix S from the generalized Schur factorization.
1.1 bertrand 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