1      SUBROUTINE QZVAL(NM,N,A,B,ALFR,ALFI,BETA,MATZ,Z)
2C***BEGIN PROLOGUE  QZVAL
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4C2C
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  The third step of the QZ algorithm for generalized
9C            eigenproblems. Accepts a pair of real matrices, one
10C            quasi-triangular form and the other in upper triangular
11C            form and computes the eigenvalues of the associated
12C            eigenproblem. Usually preceded by QZHES, QZIT, and
13C            followed by QZVEC.
14C***DESCRIPTION
15C
16C     This subroutine is the third 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
20C     This subroutine accepts a pair of REAL matrices, one of them
21C     in quasi-triangular form and the other in upper triangular form.
22C     It reduces the quasi-triangular matrix further, so that any
23C     remaining 2-by-2 blocks correspond to pairs of complex
24C     eigenvalues, and returns quantities whose ratios give the
25C     generalized eigenvalues.  It is usually preceded by  QZHES
26C     and  QZIT  and may be followed by  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 quasi-triangular matrix.
37C
38C        B contains a real upper triangular matrix.  In addition,
39C          location B(N,1) contains the tolerance quantity (EPSB)
40C          computed and saved in  QZIT.
41C
42C        MATZ should be set to .TRUE. If the right hand transformations
43C          are to be accumulated for later use in computing
44C          eigenvectors, and to .FALSE. otherwise.
45C
46C        Z contains, if MATZ has been set to .TRUE., the
47C          transformation matrix produced in the reductions by QZHES
48C          and QZIT, if performed, or else the identity matrix.
49C          If MATZ has been set to .FALSE., Z is not referenced.
50C
51C     On Output
52C
53C        A has been reduced further to a quasi-triangular matrix
54C          in which all nonzero subdiagonal elements correspond to
55C          pairs of complex eigenvalues.
56C
57C        B is still in upper triangular form, although its elements
58C          have been altered.  B(N,1) is unaltered.
59C
60C        ALFR and ALFI contain the real and imaginary parts of the
61C          diagonal elements of the triangular matrix that would be
62C          obtained if a were reduced completely to triangular form
63C          by unitary transformations.  Non-zero values of ALFI occur
64C          in pairs, the first member positive and the second negative.
65C
66C        BETA contains the diagonal elements of the corresponding B,
67C          normalized to be real and non-negative.  The generalized
68C          eigenvalues are then the ratios ((ALFR+I*ALFI)/BETA).
69C
70C        Z contains the product of the right hand transformations
71C          (for all three steps) if MATZ has been set to .TRUE.
72C
73C     Questions and comments should be directed to B. S. Garbow,
74C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
75C     ------------------------------------------------------------------
76C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
77C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
78C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
79C                 1976.
80C***ROUTINES CALLED  (NONE)
81C***END PROLOGUE  QZVAL
82C
83      INTEGER I,J,N,EN,NA,NM,NN,ISW
84      REAL A(NM,N),B(NM,N),ALFR(N),ALFI(N),BETA(N),Z(NM,N)
85      REAL C,D,E,R,S,T,AN,A1,A2,BN,CQ,CZ,DI,DR,EI,TI,TR
86      REAL U1,U2,V1,V2,A1I,A11,A12,A2I,A21,A22,B11,B12,B22
87      REAL SQI,SQR,SSI,SSR,SZI,SZR,A11I,A11R,A12I,A12R
88      REAL A22I,A22R,EPSB
89      LOGICAL MATZ
90C
91C***FIRST EXECUTABLE STATEMENT  QZVAL
92      EPSB = B(N,1)
93      ISW = 1
94C     .......... FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES.
95C                FOR EN=N STEP -1 UNTIL 1 DO -- ..........
96      DO 510 NN = 1, N
97         EN = N + 1 - NN
98         NA = EN - 1
99         IF (ISW .EQ. 2) GO TO 505
100         IF (EN .EQ. 1) GO TO 410
101         IF (A(EN,NA) .NE. 0.0E0) GO TO 420
102C     .......... 1-BY-1 BLOCK, ONE REAL ROOT ..........
103  410    ALFR(EN) = A(EN,EN)
104         IF (B(EN,EN) .LT. 0.0E0) ALFR(EN) = -ALFR(EN)
105         BETA(EN) = ABS(B(EN,EN))
106         ALFI(EN) = 0.0E0
107         GO TO 510
108C     .......... 2-BY-2 BLOCK ..........
109  420    IF (ABS(B(NA,NA)) .LE. EPSB) GO TO 455
110         IF (ABS(B(EN,EN)) .GT. EPSB) GO TO 430
111         A1 = A(EN,EN)
112         A2 = A(EN,NA)
113         BN = 0.0E0
114         GO TO 435
115  430    AN = ABS(A(NA,NA)) + ABS(A(NA,EN)) + ABS(A(EN,NA))
116     1      + ABS(A(EN,EN))
117         BN = ABS(B(NA,NA)) + ABS(B(NA,EN)) + ABS(B(EN,EN))
118         A11 = A(NA,NA) / AN
119         A12 = A(NA,EN) / AN
120         A21 = A(EN,NA) / AN
121         A22 = A(EN,EN) / AN
122         B11 = B(NA,NA) / BN
123         B12 = B(NA,EN) / BN
124         B22 = B(EN,EN) / BN
125         E = A11 / B11
126         EI = A22 / B22
127         S = A21 / (B11 * B22)
128         T = (A22 - E * B22) / B22
129         IF (ABS(E) .LE. ABS(EI)) GO TO 431
130         E = EI
131         T = (A11 - E * B11) / B11
132  431    C = 0.5E0 * (T - S * B12)
133         D = C * C + S * (A12 - E * B12)
134         IF (D .LT. 0.0E0) GO TO 480
135C     .......... TWO REAL ROOTS.
136C                ZERO BOTH A(EN,NA) AND B(EN,NA) ..........
137         E = E + (C + SIGN(SQRT(D),C))
138         A11 = A11 - E * B11
139         A12 = A12 - E * B12
140         A22 = A22 - E * B22
141         IF (ABS(A11) + ABS(A12) .LT.
142     1       ABS(A21) + ABS(A22)) GO TO 432
143         A1 = A12
144         A2 = A11
145         GO TO 435
146  432    A1 = A22
147         A2 = A21
148C     .......... CHOOSE AND APPLY REAL Z ..........
149  435    S = ABS(A1) + ABS(A2)
150         U1 = A1 / S
151         U2 = A2 / S
152         R = SIGN(SQRT(U1*U1+U2*U2),U1)
153         V1 = -(U1 + R) / R
154         V2 = -U2 / R
155         U2 = V2 / V1
156C
157         DO 440 I = 1, EN
158            T = A(I,EN) + U2 * A(I,NA)
159            A(I,EN) = A(I,EN) + T * V1
160            A(I,NA) = A(I,NA) + T * V2
161            T = B(I,EN) + U2 * B(I,NA)
162            B(I,EN) = B(I,EN) + T * V1
163            B(I,NA) = B(I,NA) + T * V2
164  440    CONTINUE
165C
166         IF (.NOT. MATZ) GO TO 450
167C
168         DO 445 I = 1, N
169            T = Z(I,EN) + U2 * Z(I,NA)
170            Z(I,EN) = Z(I,EN) + T * V1
171            Z(I,NA) = Z(I,NA) + T * V2
172  445    CONTINUE
173C
174  450    IF (BN .EQ. 0.0E0) GO TO 475
175         IF (AN .LT. ABS(E) * BN) GO TO 455
176         A1 = B(NA,NA)
177         A2 = B(EN,NA)
178         GO TO 460
179  455    A1 = A(NA,NA)
180         A2 = A(EN,NA)
181C     .......... CHOOSE AND APPLY REAL Q ..........
182  460    S = ABS(A1) + ABS(A2)
183         IF (S .EQ. 0.0E0) GO TO 475
184         U1 = A1 / S
185         U2 = A2 / S
186         R = SIGN(SQRT(U1*U1+U2*U2),U1)
187         V1 = -(U1 + R) / R
188         V2 = -U2 / R
189         U2 = V2 / V1
190C
191         DO 470 J = NA, N
192            T = A(NA,J) + U2 * A(EN,J)
193            A(NA,J) = A(NA,J) + T * V1
194            A(EN,J) = A(EN,J) + T * V2
195            T = B(NA,J) + U2 * B(EN,J)
196            B(NA,J) = B(NA,J) + T * V1
197            B(EN,J) = B(EN,J) + T * V2
198  470    CONTINUE
199C
200  475    A(EN,NA) = 0.0E0
201         B(EN,NA) = 0.0E0
202         ALFR(NA) = A(NA,NA)
203         ALFR(EN) = A(EN,EN)
204         IF (B(NA,NA) .LT. 0.0E0) ALFR(NA) = -ALFR(NA)
205         IF (B(EN,EN) .LT. 0.0E0) ALFR(EN) = -ALFR(EN)
206         BETA(NA) = ABS(B(NA,NA))
207         BETA(EN) = ABS(B(EN,EN))
208         ALFI(EN) = 0.0E0
209         ALFI(NA) = 0.0E0
210         GO TO 505
211C     .......... TWO COMPLEX ROOTS ..........
212  480    E = E + C
213         EI = SQRT(-D)
214         A11R = A11 - E * B11
215         A11I = EI * B11
216         A12R = A12 - E * B12
217         A12I = EI * B12
218         A22R = A22 - E * B22
219         A22I = EI * B22
220         IF (ABS(A11R) + ABS(A11I) + ABS(A12R) + ABS(A12I) .LT.
221     1       ABS(A21) + ABS(A22R) + ABS(A22I)) GO TO 482
222         A1 = A12R
223         A1I = A12I
224         A2 = -A11R
225         A2I = -A11I
226         GO TO 485
227  482    A1 = A22R
228         A1I = A22I
229         A2 = -A21
230         A2I = 0.0E0
231C     .......... CHOOSE COMPLEX Z ..........
232  485    CZ = SQRT(A1*A1+A1I*A1I)
233         IF (CZ .EQ. 0.0E0) GO TO 487
234         SZR = (A1 * A2 + A1I * A2I) / CZ
235         SZI = (A1 * A2I - A1I * A2) / CZ
236         R = SQRT(CZ*CZ+SZR*SZR+SZI*SZI)
237         CZ = CZ / R
238         SZR = SZR / R
239         SZI = SZI / R
240         GO TO 490
241  487    SZR = 1.0E0
242         SZI = 0.0E0
243  490    IF (AN .LT. (ABS(E) + EI) * BN) GO TO 492
244         A1 = CZ * B11 + SZR * B12
245         A1I = SZI * B12
246         A2 = SZR * B22
247         A2I = SZI * B22
248         GO TO 495
249  492    A1 = CZ * A11 + SZR * A12
250         A1I = SZI * A12
251         A2 = CZ * A21 + SZR * A22
252         A2I = SZI * A22
253C     .......... CHOOSE COMPLEX Q ..........
254  495    CQ = SQRT(A1*A1+A1I*A1I)
255         IF (CQ .EQ. 0.0E0) GO TO 497
256         SQR = (A1 * A2 + A1I * A2I) / CQ
257         SQI = (A1 * A2I - A1I * A2) / CQ
258         R = SQRT(CQ*CQ+SQR*SQR+SQI*SQI)
259         CQ = CQ / R
260         SQR = SQR / R
261         SQI = SQI / R
262         GO TO 500
263  497    SQR = 1.0E0
264         SQI = 0.0E0
265C     .......... COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT
266C                IF TRANSFORMATIONS WERE APPLIED ..........
267  500    SSR = SQR * SZR + SQI * SZI
268         SSI = SQR * SZI - SQI * SZR
269         I = 1
270         TR = CQ * CZ * A11 + CQ * SZR * A12 + SQR * CZ * A21
271     1      + SSR * A22
272         TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22
273         DR = CQ * CZ * B11 + CQ * SZR * B12 + SSR * B22
274         DI = CQ * SZI * B12 + SSI * B22
275         GO TO 503
276  502    I = 2
277         TR = SSR * A11 - SQR * CZ * A12 - CQ * SZR * A21
278     1      + CQ * CZ * A22
279         TI = -SSI * A11 - SQI * CZ * A12 + CQ * SZI * A21
280         DR = SSR * B11 - SQR * CZ * B12 + CQ * CZ * B22
281         DI = -SSI * B11 - SQI * CZ * B12
282  503    T = TI * DR - TR * DI
283         J = NA
284         IF (T .LT. 0.0E0) J = EN
285         R = SQRT(DR*DR+DI*DI)
286         BETA(J) = BN * R
287         ALFR(J) = AN * (TR * DR + TI * DI) / R
288         ALFI(J) = AN * T / R
289         IF (I .EQ. 1) GO TO 502
290  505    ISW = 3 - ISW
291  510 CONTINUE
292C
293      RETURN
294      END
295