2
SUBROUTINE QZIT (NM, N, A, B, EPS1, MATZ, Z, IERR)
3
C***BEGIN PROLOGUE QZIT
4
C***PURPOSE The second step of the QZ algorithm for generalized
5
C eigenproblems. Accepts an upper Hessenberg and an upper
6
C triangular matrix and reduces the former to
7
C quasi-triangular form while preserving the form of the
8
C latter. Usually preceded by QZHES and followed by QZVAL
10
C***LIBRARY SLATEC (EISPACK)
12
C***TYPE SINGLE PRECISION (QZIT-S)
13
C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
14
C***AUTHOR Smith, B. T., et al.
17
C This subroutine is the second step of the QZ algorithm
18
C for solving generalized matrix eigenvalue problems,
19
C SIAM J. NUMER. ANAL. 10, 241-256(1973) by MOLER and STEWART,
20
C as modified in technical note NASA TN D-7305(1973) by WARD.
22
C This subroutine accepts a pair of REAL matrices, one of them
23
C in upper Hessenberg form and the other in upper triangular form.
24
C It reduces the Hessenberg matrix to quasi-triangular form using
25
C orthogonal transformations while maintaining the triangular form
26
C of the other matrix. It is usually preceded by QZHES and
27
C followed by QZVAL and, possibly, QZVEC.
31
C NM must be set to the row dimension of the two-dimensional
32
C array parameters, A, B, and Z, as declared in the calling
33
C program dimension statement. NM is an INTEGER variable.
35
C N is the order of the matrices A and B. N is an INTEGER
36
C variable. N must be less than or equal to NM.
38
C A contains a real upper Hessenberg matrix. A is a two-
39
C dimensional REAL array, dimensioned A(NM,N).
41
C B contains a real upper triangular matrix. B is a two-
42
C dimensional REAL array, dimensioned B(NM,N).
44
C EPS1 is a tolerance used to determine negligible elements.
45
C EPS1 = 0.0 (or negative) may be input, in which case an
46
C element will be neglected only if it is less than roundoff
47
C error times the norm of its matrix. If the input EPS1 is
48
C positive, then an element will be considered negligible
49
C if it is less than EPS1 times the norm of its matrix. A
50
C positive value of EPS1 may result in faster execution,
51
C but less accurate results. EPS1 is a REAL variable.
53
C MATZ should be set to .TRUE. if the right hand transformations
54
C are to be accumulated for later use in computing
55
C eigenvectors, and to .FALSE. otherwise. MATZ is a LOGICAL
58
C Z contains, if MATZ has been set to .TRUE., the transformation
59
C matrix produced in the reduction by QZHES, if performed, or
60
C else the identity matrix. If MATZ has been set to .FALSE.,
61
C Z is not referenced. Z is a two-dimensional REAL array,
62
C dimensioned Z(NM,N).
66
C A has been reduced to quasi-triangular form. The elements
67
C below the first subdiagonal are still zero, and no two
68
C consecutive subdiagonal elements are nonzero.
70
C B is still in upper triangular form, although its elements
71
C have been altered. The location B(N,1) is used to store
72
C EPS1 times the norm of B for later use by QZVAL and QZVEC.
74
C Z contains the product of the right hand transformations
75
C (for both steps) if MATZ has been set to .TRUE.
77
C IERR is an INTEGER flag set to
78
C Zero for normal return,
79
C J if neither A(J,J-1) nor A(J-1,J-2) has become
80
C zero after a total of 30*N iterations.
82
C Questions and comments should be directed to B. S. Garbow,
83
C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
84
C ------------------------------------------------------------------
86
C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
87
C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
88
C system Routines - EISPACK Guide, Springer-Verlag,
90
C***ROUTINES CALLED (NONE)
91
C***REVISION HISTORY (YYMMDD)
93
C 890531 Changed all specific intrinsics to generic. (WRB)
94
C 890831 Modified array declarations. (WRB)
95
C 890831 REVISION DATE from Version 3.2
96
C 891214 Prologue converted to Version 4.0 format. (BAB)
97
C 920501 Reformatted the REFERENCES section. (WRB)
100
INTEGER I,J,K,L,N,EN,K1,K2,LD,LL,L1,NA,NM,ISH,ITN,ITS,KM1,LM1
101
INTEGER ENM2,IERR,LOR1,ENORN
102
REAL A(NM,*),B(NM,*),Z(NM,*)
103
REAL R,S,T,A1,A2,A3,EP,SH,U1,U2,U3,V1,V2,V3,ANI
104
REAL A11,A12,A21,A22,A33,A34,A43,A44,BNI,B11
105
REAL B12,B22,B33,B34,B44,EPSA,EPSB,EPS1,ANORM,BNORM
108
C***FIRST EXECUTABLE STATEMENT QZIT
110
C .......... COMPUTE EPSA,EPSB ..........
116
IF (I .NE. 1) ANI = ABS(A(I,I-1))
120
ANI = ANI + ABS(A(I,J))
121
BNI = BNI + ABS(B(I,J))
124
IF (ANI .GT. ANORM) ANORM = ANI
125
IF (BNI .GT. BNORM) BNORM = BNI
128
IF (ANORM .EQ. 0.0E0) ANORM = 1.0E0
129
IF (BNORM .EQ. 0.0E0) BNORM = 1.0E0
131
IF (EP .GT. 0.0E0) GO TO 50
132
C .......... COMPUTE ROUNDOFF LEVEL IF EPS1 IS ZERO ..........
135
IF (1.0E0 + EP .GT. 1.0E0) GO TO 40
138
C .......... REDUCE A TO QUASI-TRIANGULAR FORM, WHILE
139
C KEEPING B TRIANGULAR ..........
144
C .......... BEGIN QZ STEP ..........
145
60 IF (EN .LE. 2) GO TO 1001
146
IF (.NOT. MATZ) ENORN = EN
151
C .......... CHECK FOR CONVERGENCE OR REDUCIBILITY.
152
C FOR L=EN STEP -1 UNTIL 1 DO -- ..........
156
IF (L .EQ. 1) GO TO 95
157
IF (ABS(A(L,LM1)) .LE. EPSA) GO TO 90
161
IF (L .LT. NA) GO TO 95
162
C .......... 1-BY-1 OR 2-BY-2 BLOCK ISOLATED ..........
165
C .......... CHECK FOR SMALL TOP OF B ..........
169
IF (ABS(B11) .GT. EPSB) GO TO 120
171
S = ABS(A(L,L)) + ABS(A(L1,L))
174
R = SIGN(SQRT(U1*U1+U2*U2),U1)
180
T = A(L,J) + U2 * A(L1,J)
181
A(L,J) = A(L,J) + T * V1
182
A(L1,J) = A(L1,J) + T * V2
183
T = B(L,J) + U2 * B(L1,J)
184
B(L,J) = B(L,J) + T * V1
185
B(L1,J) = B(L1,J) + T * V2
188
IF (L .NE. 1) A(L,LM1) = -A(L,LM1)
192
120 A11 = A(L,L) / B11
194
IF (ISH .EQ. 1) GO TO 140
195
C .......... ITERATION STRATEGY ..........
196
IF (ITN .EQ. 0) GO TO 1000
197
IF (ITS .EQ. 10) GO TO 155
198
C .......... DETERMINE TYPE OF SHIFT ..........
200
IF (ABS(B22) .LT. EPSB) B22 = EPSB
202
IF (ABS(B33) .LT. EPSB) B33 = EPSB
204
IF (ABS(B44) .LT. EPSB) B44 = EPSB
210
T = 0.5E0 * (A43 * B34 - A33 - A44)
211
R = T * T + A34 * A43 - A33 * A44
212
IF (R .LT. 0.0E0) GO TO 150
213
C .......... DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A ..........
218
IF (ABS(S-A44) .LT. ABS(SH-A44)) SH = S
219
C .......... LOOK FOR TWO CONSECUTIVE SMALL
220
C SUB-DIAGONAL ELEMENTS OF A.
221
C FOR L=EN-2 STEP -1 UNTIL LD DO -- ..........
224
IF (L .EQ. LD) GO TO 140
228
IF (ABS(B(L,L)) .GT. EPSB) T = T - SH * B(L,L)
229
IF (ABS(A(L,LM1)) .LE. ABS(T/A(L1,L)) * EPSA) GO TO 100
234
IF (L .NE. LD) A(L,LM1) = -A(L,LM1)
236
C .......... DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A ..........
237
150 A12 = A(L,L1) / B22
240
A1 = ((A33 - A11) * (A44 - A11) - A34 * A43 + A43 * B34 * A11)
241
1 / A21 + A12 - A11 * B12
242
A2 = (A22 - A11) - A21 * B12 - (A33 - A11) - (A44 - A11)
244
A3 = A(L1+1,L1) / B22
246
C .......... AD HOC SHIFT ..........
252
IF (.NOT. MATZ) LOR1 = LD
253
C .......... MAIN LOOP ..........
255
NOTLAS = K .NE. NA .AND. ISH .EQ. 2
260
IF (NOTLAS) GO TO 190
261
C .......... ZERO A(K+1,K-1) ..........
262
IF (K .EQ. L) GO TO 170
265
170 S = ABS(A1) + ABS(A2)
266
IF (S .EQ. 0.0E0) GO TO 70
269
R = SIGN(SQRT(U1*U1+U2*U2),U1)
274
DO 180 J = KM1, ENORN
275
T = A(K,J) + U2 * A(K1,J)
276
A(K,J) = A(K,J) + T * V1
277
A(K1,J) = A(K1,J) + T * V2
278
T = B(K,J) + U2 * B(K1,J)
279
B(K,J) = B(K,J) + T * V1
280
B(K1,J) = B(K1,J) + T * V2
283
IF (K .NE. L) A(K1,KM1) = 0.0E0
285
C .......... ZERO A(K+1,K-1) AND A(K+2,K-1) ..........
286
190 IF (K .EQ. L) GO TO 200
290
200 S = ABS(A1) + ABS(A2) + ABS(A3)
291
IF (S .EQ. 0.0E0) GO TO 260
295
R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1)
302
DO 210 J = KM1, ENORN
303
T = A(K,J) + U2 * A(K1,J) + U3 * A(K2,J)
304
A(K,J) = A(K,J) + T * V1
305
A(K1,J) = A(K1,J) + T * V2
306
A(K2,J) = A(K2,J) + T * V3
307
T = B(K,J) + U2 * B(K1,J) + U3 * B(K2,J)
308
B(K,J) = B(K,J) + T * V1
309
B(K1,J) = B(K1,J) + T * V2
310
B(K2,J) = B(K2,J) + T * V3
313
IF (K .EQ. L) GO TO 220
316
C .......... ZERO B(K+2,K+1) AND B(K+2,K) ..........
317
220 S = ABS(B(K2,K2)) + ABS(B(K2,K1)) + ABS(B(K2,K))
318
IF (S .EQ. 0.0E0) GO TO 240
322
R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1)
330
T = A(I,K2) + U2 * A(I,K1) + U3 * A(I,K)
331
A(I,K2) = A(I,K2) + T * V1
332
A(I,K1) = A(I,K1) + T * V2
333
A(I,K) = A(I,K) + T * V3
334
T = B(I,K2) + U2 * B(I,K1) + U3 * B(I,K)
335
B(I,K2) = B(I,K2) + T * V1
336
B(I,K1) = B(I,K1) + T * V2
337
B(I,K) = B(I,K) + T * V3
342
IF (.NOT. MATZ) GO TO 240
345
T = Z(I,K2) + U2 * Z(I,K1) + U3 * Z(I,K)
346
Z(I,K2) = Z(I,K2) + T * V1
347
Z(I,K1) = Z(I,K1) + T * V2
348
Z(I,K) = Z(I,K) + T * V3
350
C .......... ZERO B(K+1,K) ..........
351
240 S = ABS(B(K1,K1)) + ABS(B(K1,K))
352
IF (S .EQ. 0.0E0) GO TO 260
355
R = SIGN(SQRT(U1*U1+U2*U2),U1)
361
T = A(I,K1) + U2 * A(I,K)
362
A(I,K1) = A(I,K1) + T * V1
363
A(I,K) = A(I,K) + T * V2
364
T = B(I,K1) + U2 * B(I,K)
365
B(I,K1) = B(I,K1) + T * V1
366
B(I,K) = B(I,K) + T * V2
370
IF (.NOT. MATZ) GO TO 260
373
T = Z(I,K1) + U2 * Z(I,K)
374
Z(I,K1) = Z(I,K1) + T * V1
375
Z(I,K) = Z(I,K) + T * V2
379
C .......... END QZ STEP ..........
381
C .......... SET ERROR -- NEITHER BOTTOM SUBDIAGONAL ELEMENT
382
C HAS BECOME NEGLIGIBLE AFTER 30*N ITERATIONS ..........
384
C .......... SAVE EPSB FOR USE BY QZVAL AND QZVEC ..........
385
1001 IF (N .GT. 1) B(N,1) = EPSB