1      SUBROUTINE QZIT(NM,N,A,B,EPS1,MATZ,Z,IERR)
2C***BEGIN PROLOGUE  QZIT
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4C1B3
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  The second step of the QZ algorithm for generalized
9C            eigenproblems. Accepts an upper Hessenberg and an upper
10C            triangular matrix and reduces the former to
11C            quasi-triangular form while preserving the form of the
12C            latter. Usually preceeded by QZHES and followed by QZVAL
13C            and QZVEC.
14C***DESCRIPTION
15C
16C     This subroutine is the second step of the QZ algorithm
17C     for solving generalized matrix eigenvalue problems,
18C     SIAM J. NUMER. ANAL. 10, 241-256(1973) by MOLER and STEWART,
19C     as modified in technical note NASA TN D-7305(1973) by WARD.
20C
21C     This subroutine accepts a pair of REAL matrices, one of them
22C     in upper Hessenberg form and the other in upper triangular form.
23C     It reduces the Hessenberg matrix to quasi-triangular form using
24C     orthogonal transformations while maintaining the triangular form
25C     of the other matrix.  It is usually preceded by  QZHES  and
26C     followed by  QZVAL  and, possibly,  QZVEC.
27C
28C     On Input
29C
30C        NM must be set to the row dimension of two-dimensional
31C          array parameters as declared in the calling program
32C          dimension statement.
33C
34C        N is the order of the matrices.
35C
36C        A contains a real upper Hessenberg matrix.
37C
38C        B contains a real upper triangular matrix.
39C
40C        EPS1 is a tolerance used to determine negligible elements.
41C          EPS1 = 0.0 (or negative) may be input, in which case an
42C          element will be neglected only if it is less than roundoff
43C          error times the norm of its matrix.  If the input EPS1 is
44C          positive, then an element will be considered negligible
45C          if it is less than EPS1 times the norm of its matrix.  A
46C          positive value of EPS1 may result in faster execution,
47C          but less accurate results.
48C
49C        MATZ should be set to .TRUE. If the right hand transformations
50C          are to be accumulated for later use in computing
51C          eigenvectors, and to .FALSE. otherwise.
52C
53C        Z contains, if MATZ has been set to .TRUE., the
54C          transformation matrix produced in the reduction
55C          by  QZHES, if performed, or else the identity matrix.
56C          If MATZ has been set to .FALSE., Z is not referenced.
57C
58C     On Output
59C
60C        A has been reduced to quasi-triangular form.  The elements
61C          below the first subdiagonal are still zero and no two
62C          consecutive subdiagonal elements are nonzero.
63C
64C        B is still in upper triangular form, although its elements
65C          have been altered.  The location B(N,1) is used to store
66C          EPS1 times the norm of B for later use by  QZVAL  and  QZVEC.
67C
68C        Z contains the product of the right hand transformations
69C          (for both steps) if MATZ has been set to .TRUE.
70C
71C        IERR is set to
72C          ZERO       for normal return,
73C          J          if neither A(J,J-1) nor A(J-1,J-2) has become
74C                     zero after a total of 30*N iterations.
75C
76C     Questions and comments should be directed to B. S. Garbow,
77C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
78C     ------------------------------------------------------------------
79C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
80C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
81C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
82C                 1976.
83C***ROUTINES CALLED  (NONE)
84C***END PROLOGUE  QZIT
85C
86      INTEGER I,J,K,L,N,EN,K1,K2,LD,LL,L1,NA,NM,ISH,ITN,ITS,KM1,LM1
87      INTEGER ENM2,IERR,LOR1,ENORN
88      REAL A(NM,N),B(NM,N),Z(NM,N)
89      REAL R,S,T,A1,A2,A3,EP,SH,U1,U2,U3,V1,V2,V3,ANI
90      REAL A11,A12,A21,A22,A33,A34,A43,A44,BNI,B11
91      REAL B12,B22,B33,B34,B44,EPSA,EPSB,EPS1,ANORM,BNORM
92      LOGICAL MATZ,NOTLAS
93C
94C***FIRST EXECUTABLE STATEMENT  QZIT
95      IERR = 0
96C     .......... COMPUTE EPSA,EPSB ..........
97      ANORM = 0.0E0
98      BNORM = 0.0E0
99C
100      DO 30 I = 1, N
101         ANI = 0.0E0
102         IF (I .NE. 1) ANI = ABS(A(I,I-1))
103         BNI = 0.0E0
104C
105         DO 20 J = I, N
106            ANI = ANI + ABS(A(I,J))
107            BNI = BNI + ABS(B(I,J))
108   20    CONTINUE
109C
110         IF (ANI .GT. ANORM) ANORM = ANI
111         IF (BNI .GT. BNORM) BNORM = BNI
112   30 CONTINUE
113C
114      IF (ANORM .EQ. 0.0E0) ANORM = 1.0E0
115      IF (BNORM .EQ. 0.0E0) BNORM = 1.0E0
116      EP = EPS1
117      IF (EP .GT. 0.0E0) GO TO 50
118C     .......... COMPUTE ROUNDOFF LEVEL IF EPS1 IS ZERO ..........
119      EP = 1.0E0
120   40 EP = EP / 2.0E0
121      IF (1.0E0 + EP .GT. 1.0E0) GO TO 40
122   50 EPSA = EP * ANORM
123      EPSB = EP * BNORM
124C     .......... REDUCE A TO QUASI-TRIANGULAR FORM, WHILE
125C                KEEPING B TRIANGULAR ..........
126      LOR1 = 1
127      ENORN = N
128      EN = N
129      ITN = 30*N
130C     .......... BEGIN QZ STEP ..........
131   60 IF (EN .LE. 2) GO TO 1001
132      IF (.NOT. MATZ) ENORN = EN
133      ITS = 0
134      NA = EN - 1
135      ENM2 = NA - 1
136   70 ISH = 2
137C     .......... CHECK FOR CONVERGENCE OR REDUCIBILITY.
138C                FOR L=EN STEP -1 UNTIL 1 DO -- ..........
139      DO 80 LL = 1, EN
140         LM1 = EN - LL
141         L = LM1 + 1
142         IF (L .EQ. 1) GO TO 95
143         IF (ABS(A(L,LM1)) .LE. EPSA) GO TO 90
144   80 CONTINUE
145C
146   90 A(L,LM1) = 0.0E0
147      IF (L .LT. NA) GO TO 95
148C     .......... 1-BY-1 OR 2-BY-2 BLOCK ISOLATED ..........
149      EN = LM1
150      GO TO 60
151C     .......... CHECK FOR SMALL TOP OF B ..........
152   95 LD = L
153  100 L1 = L + 1
154      B11 = B(L,L)
155      IF (ABS(B11) .GT. EPSB) GO TO 120
156      B(L,L) = 0.0E0
157      S = ABS(A(L,L)) + ABS(A(L1,L))
158      U1 = A(L,L) / S
159      U2 = A(L1,L) / S
160      R = SIGN(SQRT(U1*U1+U2*U2),U1)
161      V1 = -(U1 + R) / R
162      V2 = -U2 / R
163      U2 = V2 / V1
164C
165      DO 110 J = L, ENORN
166         T = A(L,J) + U2 * A(L1,J)
167         A(L,J) = A(L,J) + T * V1
168         A(L1,J) = A(L1,J) + T * V2
169         T = B(L,J) + U2 * B(L1,J)
170         B(L,J) = B(L,J) + T * V1
171         B(L1,J) = B(L1,J) + T * V2
172  110 CONTINUE
173C
174      IF (L .NE. 1) A(L,LM1) = -A(L,LM1)
175      LM1 = L
176      L = L1
177      GO TO 90
178  120 A11 = A(L,L) / B11
179      A21 = A(L1,L) / B11
180      IF (ISH .EQ. 1) GO TO 140
181C     .......... ITERATION STRATEGY ..........
182      IF (ITN .EQ. 0) GO TO 1000
183      IF (ITS .EQ. 10) GO TO 155
184C     .......... DETERMINE TYPE OF SHIFT ..........
185      B22 = B(L1,L1)
186      IF (ABS(B22) .LT. EPSB) B22 = EPSB
187      B33 = B(NA,NA)
188      IF (ABS(B33) .LT. EPSB) B33 = EPSB
189      B44 = B(EN,EN)
190      IF (ABS(B44) .LT. EPSB) B44 = EPSB
191      A33 = A(NA,NA) / B33
192      A34 = A(NA,EN) / B44
193      A43 = A(EN,NA) / B33
194      A44 = A(EN,EN) / B44
195      B34 = B(NA,EN) / B44
196      T = 0.5E0 * (A43 * B34 - A33 - A44)
197      R = T * T + A34 * A43 - A33 * A44
198      IF (R .LT. 0.0E0) GO TO 150
199C     .......... DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A ..........
200      ISH = 1
201      R = SQRT(R)
202      SH = -T + R
203      S = -T - R
204      IF (ABS(S-A44) .LT. ABS(SH-A44)) SH = S
205C     .......... LOOK FOR TWO CONSECUTIVE SMALL
206C                SUB-DIAGONAL ELEMENTS OF A.
207C                FOR L=EN-2 STEP -1 UNTIL LD DO -- ..........
208      DO 130 LL = LD, ENM2
209         L = ENM2 + LD - LL
210         IF (L .EQ. LD) GO TO 140
211         LM1 = L - 1
212         L1 = L + 1
213         T = A(L,L)
214         IF (ABS(B(L,L)) .GT. EPSB) T = T - SH * B(L,L)
215         IF (ABS(A(L,LM1)) .LE. ABS(T/A(L1,L)) * EPSA) GO TO 100
216  130 CONTINUE
217C
218  140 A1 = A11 - SH
219      A2 = A21
220      IF (L .NE. LD) A(L,LM1) = -A(L,LM1)
221      GO TO 160
222C     .......... DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A ..........
223  150 A12 = A(L,L1) / B22
224      A22 = A(L1,L1) / B22
225      B12 = B(L,L1) / B22
226      A1 = ((A33 - A11) * (A44 - A11) - A34 * A43 + A43 * B34 * A11)
227     1     / A21 + A12 - A11 * B12
228      A2 = (A22 - A11) - A21 * B12 - (A33 - A11) - (A44 - A11)
229     1     + A43 * B34
230      A3 = A(L1+1,L1) / B22
231      GO TO 160
232C     .......... AD HOC SHIFT ..........
233  155 A1 = 0.0E0
234      A2 = 1.0E0
235      A3 = 1.1605E0
236  160 ITS = ITS + 1
237      ITN = ITN - 1
238      IF (.NOT. MATZ) LOR1 = LD
239C     .......... MAIN LOOP ..........
240      DO 260 K = L, NA
241         NOTLAS = K .NE. NA .AND. ISH .EQ. 2
242         K1 = K + 1
243         K2 = K + 2
244         KM1 = MAX0(K-1,L)
245         LL = MIN0(EN,K1+ISH)
246         IF (NOTLAS) GO TO 190
247C     .......... ZERO A(K+1,K-1) ..........
248         IF (K .EQ. L) GO TO 170
249         A1 = A(K,KM1)
250         A2 = A(K1,KM1)
251  170    S = ABS(A1) + ABS(A2)
252         IF (S .EQ. 0.0E0) GO TO 70
253         U1 = A1 / S
254         U2 = A2 / S
255         R = SIGN(SQRT(U1*U1+U2*U2),U1)
256         V1 = -(U1 + R) / R
257         V2 = -U2 / R
258         U2 = V2 / V1
259C
260         DO 180 J = KM1, ENORN
261            T = A(K,J) + U2 * A(K1,J)
262            A(K,J) = A(K,J) + T * V1
263            A(K1,J) = A(K1,J) + T * V2
264            T = B(K,J) + U2 * B(K1,J)
265            B(K,J) = B(K,J) + T * V1
266            B(K1,J) = B(K1,J) + T * V2
267  180    CONTINUE
268C
269         IF (K .NE. L) A(K1,KM1) = 0.0E0
270         GO TO 240
271C     .......... ZERO A(K+1,K-1) AND A(K+2,K-1) ..........
272  190    IF (K .EQ. L) GO TO 200
273         A1 = A(K,KM1)
274         A2 = A(K1,KM1)
275         A3 = A(K2,KM1)
276  200    S = ABS(A1) + ABS(A2) + ABS(A3)
277         IF (S .EQ. 0.0E0) GO TO 260
278         U1 = A1 / S
279         U2 = A2 / S
280         U3 = A3 / S
281         R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1)
282         V1 = -(U1 + R) / R
283         V2 = -U2 / R
284         V3 = -U3 / R
285         U2 = V2 / V1
286         U3 = V3 / V1
287C
288         DO 210 J = KM1, ENORN
289            T = A(K,J) + U2 * A(K1,J) + U3 * A(K2,J)
290            A(K,J) = A(K,J) + T * V1
291            A(K1,J) = A(K1,J) + T * V2
292            A(K2,J) = A(K2,J) + T * V3
293            T = B(K,J) + U2 * B(K1,J) + U3 * B(K2,J)
294            B(K,J) = B(K,J) + T * V1
295            B(K1,J) = B(K1,J) + T * V2
296            B(K2,J) = B(K2,J) + T * V3
297  210    CONTINUE
298C
299         IF (K .EQ. L) GO TO 220
300         A(K1,KM1) = 0.0E0
301         A(K2,KM1) = 0.0E0
302C     .......... ZERO B(K+2,K+1) AND B(K+2,K) ..........
303  220    S = ABS(B(K2,K2)) + ABS(B(K2,K1)) + ABS(B(K2,K))
304         IF (S .EQ. 0.0E0) GO TO 240
305         U1 = B(K2,K2) / S
306         U2 = B(K2,K1) / S
307         U3 = B(K2,K) / S
308         R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1)
309         V1 = -(U1 + R) / R
310         V2 = -U2 / R
311         V3 = -U3 / R
312         U2 = V2 / V1
313         U3 = V3 / V1
314C
315         DO 230 I = LOR1, LL
316            T = A(I,K2) + U2 * A(I,K1) + U3 * A(I,K)
317            A(I,K2) = A(I,K2) + T * V1
318            A(I,K1) = A(I,K1) + T * V2
319            A(I,K) = A(I,K) + T * V3
320            T = B(I,K2) + U2 * B(I,K1) + U3 * B(I,K)
321            B(I,K2) = B(I,K2) + T * V1
322            B(I,K1) = B(I,K1) + T * V2
323            B(I,K) = B(I,K) + T * V3
324  230    CONTINUE
325C
326         B(K2,K) = 0.0E0
327         B(K2,K1) = 0.0E0
328         IF (.NOT. MATZ) GO TO 240
329C
330         DO 235 I = 1, N
331            T = Z(I,K2) + U2 * Z(I,K1) + U3 * Z(I,K)
332            Z(I,K2) = Z(I,K2) + T * V1
333            Z(I,K1) = Z(I,K1) + T * V2
334            Z(I,K) = Z(I,K) + T * V3
335  235    CONTINUE
336C     .......... ZERO B(K+1,K) ..........
337  240    S = ABS(B(K1,K1)) + ABS(B(K1,K))
338         IF (S .EQ. 0.0E0) GO TO 260
339         U1 = B(K1,K1) / S
340         U2 = B(K1,K) / S
341         R = SIGN(SQRT(U1*U1+U2*U2),U1)
342         V1 = -(U1 + R) / R
343         V2 = -U2 / R
344         U2 = V2 / V1
345C
346         DO 250 I = LOR1, LL
347            T = A(I,K1) + U2 * A(I,K)
348            A(I,K1) = A(I,K1) + T * V1
349            A(I,K) = A(I,K) + T * V2
350            T = B(I,K1) + U2 * B(I,K)
351            B(I,K1) = B(I,K1) + T * V1
352            B(I,K) = B(I,K) + T * V2
353  250    CONTINUE
354C
355         B(K1,K) = 0.0E0
356         IF (.NOT. MATZ) GO TO 260
357C
358         DO 255 I = 1, N
359            T = Z(I,K1) + U2 * Z(I,K)
360            Z(I,K1) = Z(I,K1) + T * V1
361            Z(I,K) = Z(I,K) + T * V2
362  255    CONTINUE
363C
364  260 CONTINUE
365C     .......... END QZ STEP ..........
366      GO TO 70
367C     .......... SET ERROR -- NEITHER BOTTOM SUBDIAGONAL ELEMENT
368C                HAS BECOME NEGLIGIBLE AFTER 30*N ITERATIONS ..........
369 1000 IERR = EN
370C     .......... SAVE EPSB FOR USE BY QZVAL AND QZVEC ..........
371 1001 IF (N .GT. 1) B(N,1) = EPSB
372      RETURN
373      END
374