1 /* qzval.f -- translated by f2c (version 19961017).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
qzval_(integer * nm,integer * n,doublereal * a,doublereal * b,doublereal * alfr,doublereal * alfi,doublereal * beta,logical * matz,doublereal * z__)8 /* Subroutine */ int qzval_(integer *nm, integer *n, doublereal *a,
9 doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta,
10 logical *matz, doublereal *z__)
11 {
12 /* System generated locals */
13 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i__1, i__2;
14 doublereal d__1, d__2, d__3, d__4;
15
16 /* Builtin functions */
17 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
18
19 /* Local variables */
20 doublereal epsb, c__, d__, e=0.0;
21 integer i__, j;
22 doublereal r__, s, t, a1, a2, u1, u2, v1, v2, a11, a12, a21, a22,
23 b11, b12, b22, di, ei;
24 integer na;
25 doublereal an=0.0, bn;
26 integer en;
27 doublereal cq, dr;
28 integer nn;
29 doublereal cz, ti, tr, a1i, a2i, a11i, a12i, a22i, a11r, a12r,
30 a22r, sqi, ssi;
31 integer isw;
32 doublereal sqr, szi, ssr, szr;
33
34
35
36 /* THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM */
37 /* FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */
38 /* SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */
39
40 /* THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM */
41 /* IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM. */
42 /* IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY */
43 /* REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX */
44 /* EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE */
45 /* GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY QZHES */
46 /* AND QZIT AND MAY BE FOLLOWED BY QZVEC. */
47
48 /* ON INPUT */
49
50 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
51 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
52 /* DIMENSION STATEMENT. */
53
54 /* N IS THE ORDER OF THE MATRICES. */
55
56 /* A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. */
57
58 /* B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, */
59 /* LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) */
60 /* COMPUTED AND SAVED IN QZIT. */
61
62 /* MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS
63 */
64 /* ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */
65 /* EIGENVECTORS, AND TO .FALSE. OTHERWISE. */
66
67 /* Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE */
68 /* TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES */
69 /* AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. */
70 /* IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. */
71
72 /* ON OUTPUT */
73
74 /* A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX */
75 /* IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO */
76 /* PAIRS OF COMPLEX EIGENVALUES. */
77
78 /* B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS */
79 /* HAVE BEEN ALTERED. B(N,1) IS UNALTERED. */
80
81 /* ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE */
82 /* DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE */
83 /* OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM */
84 /* BY UNITARY TRANSFORMATIONS. NON-ZERO VALUES OF ALFI OCCUR */
85 /* IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE.
86 */
87
88 /* BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B, */
89 /* NORMALIZED TO BE REAL AND NON-NEGATIVE. THE GENERALIZED */
90 /* EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA). */
91
92 /* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS */
93 /* (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE. */
94
95 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
96 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
97 */
98
99 /* THIS VERSION DATED AUGUST 1983. */
100
101 /* ------------------------------------------------------------------
102 */
103
104 /* Parameter adjustments */
105 z_dim1 = *nm;
106 z_offset = z_dim1 + 1;
107 z__ -= z_offset;
108 --beta;
109 --alfi;
110 --alfr;
111 b_dim1 = *nm;
112 b_offset = b_dim1 + 1;
113 b -= b_offset;
114 a_dim1 = *nm;
115 a_offset = a_dim1 + 1;
116 a -= a_offset;
117
118 /* Function Body */
119 epsb = b[*n + b_dim1];
120 isw = 1;
121 /* .......... FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES. */
122 /* FOR EN=N STEP -1 UNTIL 1 DO -- .......... */
123 i__1 = *n;
124 for (nn = 1; nn <= i__1; ++nn) {
125 en = *n + 1 - nn;
126 na = en - 1;
127 if (isw == 2) {
128 goto L505;
129 }
130 if (en == 1) {
131 goto L410;
132 }
133 if (a[en + na * a_dim1] != 0.) {
134 goto L420;
135 }
136 /* .......... 1-BY-1 BLOCK, ONE REAL ROOT .......... */
137 L410:
138 alfr[en] = a[en + en * a_dim1];
139 if (b[en + en * b_dim1] < 0.) {
140 alfr[en] = -alfr[en];
141 }
142 beta[en] = (d__1 = b[en + en * b_dim1], abs(d__1));
143 alfi[en] = 0.;
144 goto L510;
145 /* .......... 2-BY-2 BLOCK .......... */
146 L420:
147 if ((d__1 = b[na + na * b_dim1], abs(d__1)) <= epsb) {
148 goto L455;
149 }
150 if ((d__1 = b[en + en * b_dim1], abs(d__1)) > epsb) {
151 goto L430;
152 }
153 a1 = a[en + en * a_dim1];
154 a2 = a[en + na * a_dim1];
155 bn = 0.;
156 goto L435;
157 L430:
158 an = (d__1 = a[na + na * a_dim1], abs(d__1)) + (d__2 = a[na + en *
159 a_dim1], abs(d__2)) + (d__3 = a[en + na * a_dim1], abs(d__3))
160 + (d__4 = a[en + en * a_dim1], abs(d__4));
161 bn = (d__1 = b[na + na * b_dim1], abs(d__1)) + (d__2 = b[na + en *
162 b_dim1], abs(d__2)) + (d__3 = b[en + en * b_dim1], abs(d__3));
163 a11 = a[na + na * a_dim1] / an;
164 a12 = a[na + en * a_dim1] / an;
165 a21 = a[en + na * a_dim1] / an;
166 a22 = a[en + en * a_dim1] / an;
167 b11 = b[na + na * b_dim1] / bn;
168 b12 = b[na + en * b_dim1] / bn;
169 b22 = b[en + en * b_dim1] / bn;
170 e = a11 / b11;
171 ei = a22 / b22;
172 s = a21 / (b11 * b22);
173 t = (a22 - e * b22) / b22;
174 if (abs(e) <= abs(ei)) {
175 goto L431;
176 }
177 e = ei;
178 t = (a11 - e * b11) / b11;
179 L431:
180 c__ = (t - s * b12) * .5;
181 d__ = c__ * c__ + s * (a12 - e * b12);
182 if (d__ < 0.) {
183 goto L480;
184 }
185 /* .......... TWO REAL ROOTS. */
186 /* ZERO BOTH A(EN,NA) AND B(EN,NA) .......... */
187 d__1 = sqrt(d__);
188 e += c__ + d_sign(&d__1, &c__);
189 a11 -= e * b11;
190 a12 -= e * b12;
191 a22 -= e * b22;
192 if (abs(a11) + abs(a12) < abs(a21) + abs(a22)) {
193 goto L432;
194 }
195 a1 = a12;
196 a2 = a11;
197 goto L435;
198 L432:
199 a1 = a22;
200 a2 = a21;
201 /* .......... CHOOSE AND APPLY REAL Z .......... */
202 L435:
203 s = abs(a1) + abs(a2);
204 u1 = a1 / s;
205 u2 = a2 / s;
206 d__1 = sqrt(u1 * u1 + u2 * u2);
207 r__ = d_sign(&d__1, &u1);
208 v1 = -(u1 + r__) / r__;
209 v2 = -u2 / r__;
210 u2 = v2 / v1;
211
212 i__2 = en;
213 for (i__ = 1; i__ <= i__2; ++i__) {
214 t = a[i__ + en * a_dim1] + u2 * a[i__ + na * a_dim1];
215 a[i__ + en * a_dim1] += t * v1;
216 a[i__ + na * a_dim1] += t * v2;
217 t = b[i__ + en * b_dim1] + u2 * b[i__ + na * b_dim1];
218 b[i__ + en * b_dim1] += t * v1;
219 b[i__ + na * b_dim1] += t * v2;
220 /* L440: */
221 }
222
223 if (! (*matz)) {
224 goto L450;
225 }
226
227 i__2 = *n;
228 for (i__ = 1; i__ <= i__2; ++i__) {
229 t = z__[i__ + en * z_dim1] + u2 * z__[i__ + na * z_dim1];
230 z__[i__ + en * z_dim1] += t * v1;
231 z__[i__ + na * z_dim1] += t * v2;
232 /* L445: */
233 }
234
235 L450:
236 if (bn == 0.) {
237 goto L475;
238 }
239 if (an < abs(e) * bn) {
240 goto L455;
241 }
242 a1 = b[na + na * b_dim1];
243 a2 = b[en + na * b_dim1];
244 goto L460;
245 L455:
246 a1 = a[na + na * a_dim1];
247 a2 = a[en + na * a_dim1];
248 /* .......... CHOOSE AND APPLY REAL Q .......... */
249 L460:
250 s = abs(a1) + abs(a2);
251 if (s == 0.) {
252 goto L475;
253 }
254 u1 = a1 / s;
255 u2 = a2 / s;
256 d__1 = sqrt(u1 * u1 + u2 * u2);
257 r__ = d_sign(&d__1, &u1);
258 v1 = -(u1 + r__) / r__;
259 v2 = -u2 / r__;
260 u2 = v2 / v1;
261
262 i__2 = *n;
263 for (j = na; j <= i__2; ++j) {
264 t = a[na + j * a_dim1] + u2 * a[en + j * a_dim1];
265 a[na + j * a_dim1] += t * v1;
266 a[en + j * a_dim1] += t * v2;
267 t = b[na + j * b_dim1] + u2 * b[en + j * b_dim1];
268 b[na + j * b_dim1] += t * v1;
269 b[en + j * b_dim1] += t * v2;
270 /* L470: */
271 }
272
273 L475:
274 a[en + na * a_dim1] = 0.;
275 b[en + na * b_dim1] = 0.;
276 alfr[na] = a[na + na * a_dim1];
277 alfr[en] = a[en + en * a_dim1];
278 if (b[na + na * b_dim1] < 0.) {
279 alfr[na] = -alfr[na];
280 }
281 if (b[en + en * b_dim1] < 0.) {
282 alfr[en] = -alfr[en];
283 }
284 beta[na] = (d__1 = b[na + na * b_dim1], abs(d__1));
285 beta[en] = (d__1 = b[en + en * b_dim1], abs(d__1));
286 alfi[en] = 0.;
287 alfi[na] = 0.;
288 goto L505;
289 /* .......... TWO COMPLEX ROOTS .......... */
290 L480:
291 e += c__;
292 ei = sqrt(-d__);
293 a11r = a11 - e * b11;
294 a11i = ei * b11;
295 a12r = a12 - e * b12;
296 a12i = ei * b12;
297 a22r = a22 - e * b22;
298 a22i = ei * b22;
299 if (abs(a11r) + abs(a11i) + abs(a12r) + abs(a12i) < abs(a21) + abs(
300 a22r) + abs(a22i)) {
301 goto L482;
302 }
303 a1 = a12r;
304 a1i = a12i;
305 a2 = -a11r;
306 a2i = -a11i;
307 goto L485;
308 L482:
309 a1 = a22r;
310 a1i = a22i;
311 a2 = -a21;
312 a2i = 0.;
313 /* .......... CHOOSE COMPLEX Z .......... */
314 L485:
315 cz = sqrt(a1 * a1 + a1i * a1i);
316 if (cz == 0.) {
317 goto L487;
318 }
319 szr = (a1 * a2 + a1i * a2i) / cz;
320 szi = (a1 * a2i - a1i * a2) / cz;
321 r__ = sqrt(cz * cz + szr * szr + szi * szi);
322 cz /= r__;
323 szr /= r__;
324 szi /= r__;
325 goto L490;
326 L487:
327 szr = 1.;
328 szi = 0.;
329 L490:
330 if (an < (abs(e) + ei) * bn) {
331 goto L492;
332 }
333 a1 = cz * b11 + szr * b12;
334 a1i = szi * b12;
335 a2 = szr * b22;
336 a2i = szi * b22;
337 goto L495;
338 L492:
339 a1 = cz * a11 + szr * a12;
340 a1i = szi * a12;
341 a2 = cz * a21 + szr * a22;
342 a2i = szi * a22;
343 /* .......... CHOOSE COMPLEX Q .......... */
344 L495:
345 cq = sqrt(a1 * a1 + a1i * a1i);
346 if (cq == 0.) {
347 goto L497;
348 }
349 sqr = (a1 * a2 + a1i * a2i) / cq;
350 sqi = (a1 * a2i - a1i * a2) / cq;
351 r__ = sqrt(cq * cq + sqr * sqr + sqi * sqi);
352 cq /= r__;
353 sqr /= r__;
354 sqi /= r__;
355 goto L500;
356 L497:
357 sqr = 1.;
358 sqi = 0.;
359 /* .......... COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT */
360 /* IF TRANSFORMATIONS WERE APPLIED .......... */
361 L500:
362 ssr = sqr * szr + sqi * szi;
363 ssi = sqr * szi - sqi * szr;
364 i__ = 1;
365 tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 + ssr * a22;
366 ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22;
367 dr = cq * cz * b11 + cq * szr * b12 + ssr * b22;
368 di = cq * szi * b12 + ssi * b22;
369 goto L503;
370 L502:
371 i__ = 2;
372 tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 + cq * cz * a22;
373 ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21;
374 dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22;
375 di = -ssi * b11 - sqi * cz * b12;
376 L503:
377 t = ti * dr - tr * di;
378 j = na;
379 if (t < 0.) {
380 j = en;
381 }
382 r__ = sqrt(dr * dr + di * di);
383 beta[j] = bn * r__;
384 alfr[j] = an * (tr * dr + ti * di) / r__;
385 alfi[j] = an * t / r__;
386 if (i__ == 1) {
387 goto L502;
388 }
389 L505:
390 isw = 3 - isw;
391 L510:
392 ;
393 }
394 b[*n + b_dim1] = epsb;
395
396 return 0;
397 } /* qzval_ */
398
399