2
SUBROUTINE QZVAL (NM, N, A, B, ALFR, ALFI, BETA, MATZ, Z)
3
C***BEGIN PROLOGUE QZVAL
4
C***PURPOSE The third step of the QZ algorithm for generalized
5
C eigenproblems. Accepts a pair of real matrices, one in
6
C quasi-triangular form and the other in upper triangular
7
C form and computes the eigenvalues of the associated
8
C eigenproblem. Usually preceded by QZHES, QZIT, and
10
C***LIBRARY SLATEC (EISPACK)
12
C***TYPE SINGLE PRECISION (QZVAL-S)
13
C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
14
C***AUTHOR Smith, B. T., et al.
17
C This subroutine is the third 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.
21
C This subroutine accepts a pair of REAL matrices, one of them
22
C in quasi-triangular form and the other in upper triangular form.
23
C It reduces the quasi-triangular matrix further, so that any
24
C remaining 2-by-2 blocks correspond to pairs of complex
25
C eigenvalues, and returns quantities whose ratios give the
26
C generalized eigenvalues. It is usually preceded by QZHES
27
C and QZIT and may be followed by 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 quasi-triangular matrix. A is a two-
39
C dimensional REAL array, dimensioned A(NM,N).
41
C B contains a real upper triangular matrix. In addition,
42
C location B(N,1) contains the tolerance quantity (EPSB)
43
C computed and saved in QZIT. B is a two-dimensional REAL
44
C array, dimensioned B(NM,N).
46
C MATZ should be set to .TRUE. if the right hand transformations
47
C are to be accumulated for later use in computing
48
C eigenvectors, and to .FALSE. otherwise. MATZ is a LOGICAL
51
C Z contains, if MATZ has been set to .TRUE., the transformation
52
C matrix produced in the reductions by QZHES and QZIT, if
53
C performed, or else the identity matrix. If MATZ has been set
54
C to .FALSE., Z is not referenced. Z is a two-dimensional REAL
55
C array, dimensioned Z(NM,N).
59
C A has been reduced further to a quasi-triangular matrix in
60
C which all nonzero subdiagonal elements correspond to pairs
61
C of complex eigenvalues.
63
C B is still in upper triangular form, although its elements
64
C have been altered. B(N,1) is unaltered.
66
C ALFR and ALFI contain the real and imaginary parts of the
67
C diagonal elements of the triangular matrix that would be
68
C obtained if A were reduced completely to triangular form
69
C by unitary transformations. Non-zero values of ALFI occur
70
C in pairs, the first member positive and the second negative.
71
C ALFR and ALFI are one-dimensional REAL arrays, dimensioned
72
C ALFR(N) and ALFI(N).
74
C BETA contains the diagonal elements of the corresponding B,
75
C normalized to be real and non-negative. The generalized
76
C eigenvalues are then the ratios ((ALFR+I*ALFI)/BETA).
77
C BETA is a one-dimensional REAL array, dimensioned BETA(N).
79
C Z contains the product of the right hand transformations
80
C (for all three steps) if MATZ has been set to .TRUE.
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 890831 Modified array declarations. (WRB)
94
C 890831 REVISION DATE from Version 3.2
95
C 891214 Prologue converted to Version 4.0 format. (BAB)
96
C 920501 Reformatted the REFERENCES section. (WRB)
97
C***END PROLOGUE QZVAL
99
INTEGER I,J,N,EN,NA,NM,NN,ISW
100
REAL A(NM,*),B(NM,*),ALFR(*),ALFI(*),BETA(*),Z(NM,*)
101
REAL C,D,E,R,S,T,AN,A1,A2,BN,CQ,CZ,DI,DR,EI,TI,TR
102
REAL U1,U2,V1,V2,A1I,A11,A12,A2I,A21,A22,B11,B12,B22
103
REAL SQI,SQR,SSI,SSR,SZI,SZR,A11I,A11R,A12I,A12R
107
C***FIRST EXECUTABLE STATEMENT QZVAL
110
C .......... FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES.
111
C FOR EN=N STEP -1 UNTIL 1 DO -- ..........
115
IF (ISW .EQ. 2) GO TO 505
116
IF (EN .EQ. 1) GO TO 410
117
IF (A(EN,NA) .NE. 0.0E0) GO TO 420
118
C .......... 1-BY-1 BLOCK, ONE REAL ROOT ..........
119
410 ALFR(EN) = A(EN,EN)
120
IF (B(EN,EN) .LT. 0.0E0) ALFR(EN) = -ALFR(EN)
121
BETA(EN) = ABS(B(EN,EN))
124
C .......... 2-BY-2 BLOCK ..........
125
420 IF (ABS(B(NA,NA)) .LE. EPSB) GO TO 455
126
IF (ABS(B(EN,EN)) .GT. EPSB) GO TO 430
131
430 AN = ABS(A(NA,NA)) + ABS(A(NA,EN)) + ABS(A(EN,NA))
133
BN = ABS(B(NA,NA)) + ABS(B(NA,EN)) + ABS(B(EN,EN))
143
S = A21 / (B11 * B22)
144
T = (A22 - E * B22) / B22
145
IF (ABS(E) .LE. ABS(EI)) GO TO 431
147
T = (A11 - E * B11) / B11
148
431 C = 0.5E0 * (T - S * B12)
149
D = C * C + S * (A12 - E * B12)
150
IF (D .LT. 0.0E0) GO TO 480
151
C .......... TWO REAL ROOTS.
152
C ZERO BOTH A(EN,NA) AND B(EN,NA) ..........
153
E = E + (C + SIGN(SQRT(D),C))
157
IF (ABS(A11) + ABS(A12) .LT.
158
1 ABS(A21) + ABS(A22)) GO TO 432
164
C .......... CHOOSE AND APPLY REAL Z ..........
165
435 S = ABS(A1) + ABS(A2)
168
R = SIGN(SQRT(U1*U1+U2*U2),U1)
174
T = A(I,EN) + U2 * A(I,NA)
175
A(I,EN) = A(I,EN) + T * V1
176
A(I,NA) = A(I,NA) + T * V2
177
T = B(I,EN) + U2 * B(I,NA)
178
B(I,EN) = B(I,EN) + T * V1
179
B(I,NA) = B(I,NA) + T * V2
182
IF (.NOT. MATZ) GO TO 450
185
T = Z(I,EN) + U2 * Z(I,NA)
186
Z(I,EN) = Z(I,EN) + T * V1
187
Z(I,NA) = Z(I,NA) + T * V2
190
450 IF (BN .EQ. 0.0E0) GO TO 475
191
IF (AN .LT. ABS(E) * BN) GO TO 455
197
C .......... CHOOSE AND APPLY REAL Q ..........
198
460 S = ABS(A1) + ABS(A2)
199
IF (S .EQ. 0.0E0) GO TO 475
202
R = SIGN(SQRT(U1*U1+U2*U2),U1)
208
T = A(NA,J) + U2 * A(EN,J)
209
A(NA,J) = A(NA,J) + T * V1
210
A(EN,J) = A(EN,J) + T * V2
211
T = B(NA,J) + U2 * B(EN,J)
212
B(NA,J) = B(NA,J) + T * V1
213
B(EN,J) = B(EN,J) + T * V2
220
IF (B(NA,NA) .LT. 0.0E0) ALFR(NA) = -ALFR(NA)
221
IF (B(EN,EN) .LT. 0.0E0) ALFR(EN) = -ALFR(EN)
222
BETA(NA) = ABS(B(NA,NA))
223
BETA(EN) = ABS(B(EN,EN))
227
C .......... TWO COMPLEX ROOTS ..........
236
IF (ABS(A11R) + ABS(A11I) + ABS(A12R) + ABS(A12I) .LT.
237
1 ABS(A21) + ABS(A22R) + ABS(A22I)) GO TO 482
247
C .......... CHOOSE COMPLEX Z ..........
248
485 CZ = SQRT(A1*A1+A1I*A1I)
249
IF (CZ .EQ. 0.0E0) GO TO 487
250
SZR = (A1 * A2 + A1I * A2I) / CZ
251
SZI = (A1 * A2I - A1I * A2) / CZ
252
R = SQRT(CZ*CZ+SZR*SZR+SZI*SZI)
259
490 IF (AN .LT. (ABS(E) + EI) * BN) GO TO 492
260
A1 = CZ * B11 + SZR * B12
265
492 A1 = CZ * A11 + SZR * A12
267
A2 = CZ * A21 + SZR * A22
269
C .......... CHOOSE COMPLEX Q ..........
270
495 CQ = SQRT(A1*A1+A1I*A1I)
271
IF (CQ .EQ. 0.0E0) GO TO 497
272
SQR = (A1 * A2 + A1I * A2I) / CQ
273
SQI = (A1 * A2I - A1I * A2) / CQ
274
R = SQRT(CQ*CQ+SQR*SQR+SQI*SQI)
281
C .......... COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT
282
C IF TRANSFORMATIONS WERE APPLIED ..........
283
500 SSR = SQR * SZR + SQI * SZI
284
SSI = SQR * SZI - SQI * SZR
286
TR = CQ * CZ * A11 + CQ * SZR * A12 + SQR * CZ * A21
288
TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22
289
DR = CQ * CZ * B11 + CQ * SZR * B12 + SSR * B22
290
DI = CQ * SZI * B12 + SSI * B22
293
TR = SSR * A11 - SQR * CZ * A12 - CQ * SZR * A21
295
TI = -SSI * A11 - SQI * CZ * A12 + CQ * SZI * A21
296
DR = SSR * B11 - SQR * CZ * B12 + CQ * CZ * B22
297
DI = -SSI * B11 - SQI * CZ * B12
298
503 T = TI * DR - TR * DI
300
IF (T .LT. 0.0E0) J = EN
301
R = SQRT(DR*DR+DI*DI)
303
ALFR(J) = AN * (TR * DR + TI * DI) / R
305
IF (I .EQ. 1) GO TO 502
added
removed
deps/openlibm/slatec/qzval.f
