1 /* hqr2.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
8 /* Table of constant values */
9
10 static doublereal c_b49 = 0.;
11
hqr2_(integer * nm,integer * n,integer * low,integer * igh,doublereal * h__,doublereal * wr,doublereal * wi,doublereal * z__,integer * ierr)12 /* Subroutine */ int hqr2_(integer *nm, integer *n, integer *low, integer *
13 igh, doublereal *h__, doublereal *wr, doublereal *wi, doublereal *z__,
14 integer *ierr)
15 {
16 /* System generated locals */
17 integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3;
18 doublereal d__1, d__2, d__3, d__4;
19
20 /* Builtin functions */
21 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
22
23 /* Local variables */
24 extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
25 , doublereal *, doublereal *, doublereal *);
26 doublereal norm;
27 integer i__, j, k, l=0, m=0;
28 doublereal p, q, r__=0.0, s=0.0, t, w, x, y;
29 integer na, ii, en, jj;
30 doublereal ra, sa;
31 integer ll, mm, nn;
32 doublereal vi, vr, zz;
33 logical notlas;
34 integer mp2, itn, its, enm2;
35 doublereal tst1, tst2;
36
37
38
39 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, */
40 /* NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
41 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
42
43 /* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
44 /* OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD. THE */
45 /* EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND */
46 /* IF ELMHES AND ELTRAN OR ORTHES AND ORTRAN HAVE */
47 /* BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM */
48 /* AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS. */
49
50 /* ON INPUT */
51
52 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
53 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
54 /* DIMENSION STATEMENT. */
55
56 /* N IS THE ORDER OF THE MATRIX. */
57
58 /* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
59 /* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
60 /* SET LOW=1, IGH=N. */
61
62 /* H CONTAINS THE UPPER HESSENBERG MATRIX. */
63
64 /* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ELTRAN */
65 /* AFTER THE REDUCTION BY ELMHES, OR BY ORTRAN AFTER THE */
66 /* REDUCTION BY ORTHES, IF PERFORMED. IF THE EIGENVECTORS */
67 /* OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE */
68 /* IDENTITY MATRIX. */
69
70 /* ON OUTPUT */
71
72 /* H HAS BEEN DESTROYED. */
73
74 /* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
75 /* RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES */
76 /* ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */
77 /* OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */
78 /* HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN */
79 /* ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
80 /* FOR INDICES IERR+1,...,N. */
81
82 /* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. */
83 /* IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z */
84 /* CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX
85 */
86 /* WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH */
87 /* COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS */
88 /* EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN */
89 /* ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND.
90 */
91
92 /* IERR IS SET TO */
93 /* ZERO FOR NORMAL RETURN, */
94 /* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
95 /* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
96
97 /* CALLS CDIV FOR COMPLEX DIVISION. */
98
99 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
100 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
101 */
102
103 /* THIS VERSION DATED AUGUST 1983. */
104
105 /* ------------------------------------------------------------------
106 */
107
108 /* Parameter adjustments */
109 z_dim1 = *nm;
110 z_offset = z_dim1 + 1;
111 z__ -= z_offset;
112 --wi;
113 --wr;
114 h_dim1 = *nm;
115 h_offset = h_dim1 + 1;
116 h__ -= h_offset;
117
118 /* Function Body */
119 *ierr = 0;
120 norm = 0.;
121 k = 1;
122 /* .......... STORE ROOTS ISOLATED BY BALANC */
123 /* AND COMPUTE MATRIX NORM .......... */
124 i__1 = *n;
125 for (i__ = 1; i__ <= i__1; ++i__) {
126
127 i__2 = *n;
128 for (j = k; j <= i__2; ++j) {
129 /* L40: */
130 norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1));
131 }
132
133 k = i__;
134 if (i__ >= *low && i__ <= *igh) {
135 goto L50;
136 }
137 wr[i__] = h__[i__ + i__ * h_dim1];
138 wi[i__] = 0.;
139 L50:
140 ;
141 }
142
143 en = *igh;
144 t = 0.;
145 itn = *n * 30;
146 /* .......... SEARCH FOR NEXT EIGENVALUES .......... */
147 L60:
148 if (en < *low) {
149 goto L340;
150 }
151 its = 0;
152 na = en - 1;
153 enm2 = na - 1;
154 /* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
155 /* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
156 L70:
157 i__1 = en;
158 for (ll = *low; ll <= i__1; ++ll) {
159 l = en + *low - ll;
160 if (l == *low) {
161 goto L100;
162 }
163 s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l
164 + l * h_dim1], abs(d__2));
165 if (s == 0.) {
166 s = norm;
167 }
168 tst1 = s;
169 tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1));
170 if (tst2 == tst1) {
171 goto L100;
172 }
173 /* L80: */
174 }
175 /* .......... FORM SHIFT .......... */
176 L100:
177 x = h__[en + en * h_dim1];
178 if (l == en) {
179 goto L270;
180 }
181 y = h__[na + na * h_dim1];
182 w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
183 if (l == na) {
184 goto L280;
185 }
186 if (itn == 0) {
187 goto L1000;
188 }
189 if (its != 10 && its != 20) {
190 goto L130;
191 }
192 /* .......... FORM EXCEPTIONAL SHIFT .......... */
193 t += x;
194
195 i__1 = en;
196 for (i__ = *low; i__ <= i__1; ++i__) {
197 /* L120: */
198 h__[i__ + i__ * h_dim1] -= x;
199 }
200
201 s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 *
202 h_dim1], abs(d__2));
203 x = s * .75;
204 y = x;
205 w = s * -.4375 * s;
206 L130:
207 ++its;
208 --itn;
209 /* .......... LOOK FOR TWO CONSECUTIVE SMALL */
210 /* SUB-DIAGONAL ELEMENTS. */
211 /* FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */
212 i__1 = enm2;
213 for (mm = l; mm <= i__1; ++mm) {
214 m = enm2 + l - mm;
215 zz = h__[m + m * h_dim1];
216 r__ = x - zz;
217 s = y - zz;
218 p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) *
219 h_dim1];
220 q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
221 r__ = h__[m + 2 + (m + 1) * h_dim1];
222 s = abs(p) + abs(q) + abs(r__);
223 p /= s;
224 q /= s;
225 r__ /= s;
226 if (m == l) {
227 goto L150;
228 }
229 tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) +
230 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
231 tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q)
232 + abs(r__));
233 if (tst2 == tst1) {
234 goto L150;
235 }
236 /* L140: */
237 }
238
239 L150:
240 mp2 = m + 2;
241
242 i__1 = en;
243 for (i__ = mp2; i__ <= i__1; ++i__) {
244 h__[i__ + (i__ - 2) * h_dim1] = 0.;
245 if (i__ == mp2) {
246 goto L160;
247 }
248 h__[i__ + (i__ - 3) * h_dim1] = 0.;
249 L160:
250 ;
251 }
252 /* .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */
253 /* COLUMNS M TO EN .......... */
254 i__1 = na;
255 for (k = m; k <= i__1; ++k) {
256 notlas = k != na;
257 if (k == m) {
258 goto L170;
259 }
260 p = h__[k + (k - 1) * h_dim1];
261 q = h__[k + 1 + (k - 1) * h_dim1];
262 r__ = 0.;
263 if (notlas) {
264 r__ = h__[k + 2 + (k - 1) * h_dim1];
265 }
266 x = abs(p) + abs(q) + abs(r__);
267 if (x == 0.) {
268 goto L260;
269 }
270 p /= x;
271 q /= x;
272 r__ /= x;
273 L170:
274 d__1 = sqrt(p * p + q * q + r__ * r__);
275 s = d_sign(&d__1, &p);
276 if (k == m) {
277 goto L180;
278 }
279 h__[k + (k - 1) * h_dim1] = -s * x;
280 goto L190;
281 L180:
282 if (l != m) {
283 h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
284 }
285 L190:
286 p += s;
287 x = p / s;
288 y = q / s;
289 zz = r__ / s;
290 q /= p;
291 r__ /= p;
292 if (notlas) {
293 goto L225;
294 }
295 /* .......... ROW MODIFICATION .......... */
296 i__2 = *n;
297 for (j = k; j <= i__2; ++j) {
298 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
299 h__[k + j * h_dim1] -= p * x;
300 h__[k + 1 + j * h_dim1] -= p * y;
301 /* L200: */
302 }
303
304 /* Computing MIN */
305 i__2 = en, i__3 = k + 3;
306 j = min(i__2,i__3);
307 /* .......... COLUMN MODIFICATION .......... */
308 i__2 = j;
309 for (i__ = 1; i__ <= i__2; ++i__) {
310 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
311 h__[i__ + k * h_dim1] -= p;
312 h__[i__ + (k + 1) * h_dim1] -= p * q;
313 /* L210: */
314 }
315 /* .......... ACCUMULATE TRANSFORMATIONS .......... */
316 i__2 = *igh;
317 for (i__ = *low; i__ <= i__2; ++i__) {
318 p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1];
319 z__[i__ + k * z_dim1] -= p;
320 z__[i__ + (k + 1) * z_dim1] -= p * q;
321 /* L220: */
322 }
323 goto L255;
324 L225:
325 /* .......... ROW MODIFICATION .......... */
326 i__2 = *n;
327 for (j = k; j <= i__2; ++j) {
328 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
329 k + 2 + j * h_dim1];
330 h__[k + j * h_dim1] -= p * x;
331 h__[k + 1 + j * h_dim1] -= p * y;
332 h__[k + 2 + j * h_dim1] -= p * zz;
333 /* L230: */
334 }
335
336 /* Computing MIN */
337 i__2 = en, i__3 = k + 3;
338 j = min(i__2,i__3);
339 /* .......... COLUMN MODIFICATION .......... */
340 i__2 = j;
341 for (i__ = 1; i__ <= i__2; ++i__) {
342 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] +
343 zz * h__[i__ + (k + 2) * h_dim1];
344 h__[i__ + k * h_dim1] -= p;
345 h__[i__ + (k + 1) * h_dim1] -= p * q;
346 h__[i__ + (k + 2) * h_dim1] -= p * r__;
347 /* L240: */
348 }
349 /* .......... ACCUMULATE TRANSFORMATIONS .......... */
350 i__2 = *igh;
351 for (i__ = *low; i__ <= i__2; ++i__) {
352 p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1] +
353 zz * z__[i__ + (k + 2) * z_dim1];
354 z__[i__ + k * z_dim1] -= p;
355 z__[i__ + (k + 1) * z_dim1] -= p * q;
356 z__[i__ + (k + 2) * z_dim1] -= p * r__;
357 /* L250: */
358 }
359 L255:
360
361 L260:
362 ;
363 }
364
365 goto L70;
366 /* .......... ONE ROOT FOUND .......... */
367 L270:
368 h__[en + en * h_dim1] = x + t;
369 wr[en] = h__[en + en * h_dim1];
370 wi[en] = 0.;
371 en = na;
372 goto L60;
373 /* .......... TWO ROOTS FOUND .......... */
374 L280:
375 p = (y - x) / 2.;
376 q = p * p + w;
377 zz = sqrt((abs(q)));
378 h__[en + en * h_dim1] = x + t;
379 x = h__[en + en * h_dim1];
380 h__[na + na * h_dim1] = y + t;
381 if (q < 0.) {
382 goto L320;
383 }
384 /* .......... REAL PAIR .......... */
385 zz = p + d_sign(&zz, &p);
386 wr[na] = x + zz;
387 wr[en] = wr[na];
388 if (zz != 0.) {
389 wr[en] = x - w / zz;
390 }
391 wi[na] = 0.;
392 wi[en] = 0.;
393 x = h__[en + na * h_dim1];
394 s = abs(x) + abs(zz);
395 p = x / s;
396 q = zz / s;
397 r__ = sqrt(p * p + q * q);
398 p /= r__;
399 q /= r__;
400 /* .......... ROW MODIFICATION .......... */
401 i__1 = *n;
402 for (j = na; j <= i__1; ++j) {
403 zz = h__[na + j * h_dim1];
404 h__[na + j * h_dim1] = q * zz + p * h__[en + j * h_dim1];
405 h__[en + j * h_dim1] = q * h__[en + j * h_dim1] - p * zz;
406 /* L290: */
407 }
408 /* .......... COLUMN MODIFICATION .......... */
409 i__1 = en;
410 for (i__ = 1; i__ <= i__1; ++i__) {
411 zz = h__[i__ + na * h_dim1];
412 h__[i__ + na * h_dim1] = q * zz + p * h__[i__ + en * h_dim1];
413 h__[i__ + en * h_dim1] = q * h__[i__ + en * h_dim1] - p * zz;
414 /* L300: */
415 }
416 /* .......... ACCUMULATE TRANSFORMATIONS .......... */
417 i__1 = *igh;
418 for (i__ = *low; i__ <= i__1; ++i__) {
419 zz = z__[i__ + na * z_dim1];
420 z__[i__ + na * z_dim1] = q * zz + p * z__[i__ + en * z_dim1];
421 z__[i__ + en * z_dim1] = q * z__[i__ + en * z_dim1] - p * zz;
422 /* L310: */
423 }
424
425 goto L330;
426 /* .......... COMPLEX PAIR .......... */
427 L320:
428 wr[na] = x + p;
429 wr[en] = x + p;
430 wi[na] = zz;
431 wi[en] = -zz;
432 L330:
433 en = enm2;
434 goto L60;
435 /* .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND */
436 /* VECTORS OF UPPER TRIANGULAR FORM .......... */
437 L340:
438 if (norm == 0.) {
439 goto L1001;
440 }
441 /* .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... */
442 i__1 = *n;
443 for (nn = 1; nn <= i__1; ++nn) {
444 en = *n + 1 - nn;
445 p = wr[en];
446 q = wi[en];
447 na = en - 1;
448 if (q < 0.) {
449 goto L710;
450 } else if (q == 0) {
451 goto L600;
452 } else {
453 goto L800;
454 }
455 /* .......... REAL VECTOR .......... */
456 L600:
457 m = en;
458 h__[en + en * h_dim1] = 1.;
459 if (na == 0) {
460 goto L800;
461 }
462 /* .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
463 i__2 = na;
464 for (ii = 1; ii <= i__2; ++ii) {
465 i__ = en - ii;
466 w = h__[i__ + i__ * h_dim1] - p;
467 r__ = 0.;
468
469 i__3 = en;
470 for (j = m; j <= i__3; ++j) {
471 /* L610: */
472 r__ += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
473 }
474
475 if (wi[i__] >= 0.) {
476 goto L630;
477 }
478 zz = w;
479 s = r__;
480 goto L700;
481 L630:
482 m = i__;
483 if (wi[i__] != 0.) {
484 goto L640;
485 }
486 t = w;
487 if (t != 0.) {
488 goto L635;
489 }
490 tst1 = norm;
491 t = tst1;
492 L632:
493 t *= .01;
494 tst2 = norm + t;
495 if (tst2 > tst1) {
496 goto L632;
497 }
498 L635:
499 h__[i__ + en * h_dim1] = -r__ / t;
500 goto L680;
501 /* .......... SOLVE REAL EQUATIONS .......... */
502 L640:
503 x = h__[i__ + (i__ + 1) * h_dim1];
504 y = h__[i__ + 1 + i__ * h_dim1];
505 q = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__];
506 t = (x * s - zz * r__) / q;
507 h__[i__ + en * h_dim1] = t;
508 if (abs(x) <= abs(zz)) {
509 goto L650;
510 }
511 h__[i__ + 1 + en * h_dim1] = (-r__ - w * t) / x;
512 goto L680;
513 L650:
514 h__[i__ + 1 + en * h_dim1] = (-s - y * t) / zz;
515
516 /* .......... OVERFLOW CONTROL .......... */
517 L680:
518 t = (d__1 = h__[i__ + en * h_dim1], abs(d__1));
519 if (t == 0.) {
520 goto L700;
521 }
522 tst1 = t;
523 tst2 = tst1 + 1. / tst1;
524 if (tst2 > tst1) {
525 goto L700;
526 }
527 i__3 = en;
528 for (j = i__; j <= i__3; ++j) {
529 h__[j + en * h_dim1] /= t;
530 /* L690: */
531 }
532
533 L700:
534 ;
535 }
536 /* .......... END REAL VECTOR .......... */
537 goto L800;
538 /* .......... COMPLEX VECTOR .......... */
539 L710:
540 m = na;
541 /* .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT */
542 /* EIGENVECTOR MATRIX IS TRIANGULAR .......... */
543 if ((d__1 = h__[en + na * h_dim1], abs(d__1)) <= (d__2 = h__[na + en *
544 h_dim1], abs(d__2))) {
545 goto L720;
546 }
547 h__[na + na * h_dim1] = q / h__[en + na * h_dim1];
548 h__[na + en * h_dim1] = -(h__[en + en * h_dim1] - p) / h__[en + na *
549 h_dim1];
550 goto L730;
551 L720:
552 d__1 = -h__[na + en * h_dim1];
553 d__2 = h__[na + na * h_dim1] - p;
554 cdiv_(&c_b49, &d__1, &d__2, &q, &h__[na + na * h_dim1], &h__[na + en *
555 h_dim1]);
556 L730:
557 h__[en + na * h_dim1] = 0.;
558 h__[en + en * h_dim1] = 1.;
559 enm2 = na - 1;
560 if (enm2 == 0) {
561 goto L800;
562 }
563 /* .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... */
564 i__2 = enm2;
565 for (ii = 1; ii <= i__2; ++ii) {
566 i__ = na - ii;
567 w = h__[i__ + i__ * h_dim1] - p;
568 ra = 0.;
569 sa = 0.;
570
571 i__3 = en;
572 for (j = m; j <= i__3; ++j) {
573 ra += h__[i__ + j * h_dim1] * h__[j + na * h_dim1];
574 sa += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
575 /* L760: */
576 }
577
578 if (wi[i__] >= 0.) {
579 goto L770;
580 }
581 zz = w;
582 r__ = ra;
583 s = sa;
584 goto L795;
585 L770:
586 m = i__;
587 if (wi[i__] != 0.) {
588 goto L780;
589 }
590 d__1 = -ra;
591 d__2 = -sa;
592 cdiv_(&d__1, &d__2, &w, &q, &h__[i__ + na * h_dim1], &h__[i__ +
593 en * h_dim1]);
594 goto L790;
595 /* .......... SOLVE COMPLEX EQUATIONS .......... */
596 L780:
597 x = h__[i__ + (i__ + 1) * h_dim1];
598 y = h__[i__ + 1 + i__ * h_dim1];
599 vr = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__] - q * q;
600 vi = (wr[i__] - p) * 2. * q;
601 if (vr != 0. || vi != 0.) {
602 goto L784;
603 }
604 tst1 = norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(zz));
605 vr = tst1;
606 L783:
607 vr *= .01;
608 tst2 = tst1 + vr;
609 if (tst2 > tst1) {
610 goto L783;
611 }
612 L784:
613 d__1 = x * r__ - zz * ra + q * sa;
614 d__2 = x * s - zz * sa - q * ra;
615 cdiv_(&d__1, &d__2, &vr, &vi, &h__[i__ + na * h_dim1], &h__[i__ +
616 en * h_dim1]);
617 if (abs(x) <= abs(zz) + abs(q)) {
618 goto L785;
619 }
620 h__[i__ + 1 + na * h_dim1] = (-ra - w * h__[i__ + na * h_dim1] +
621 q * h__[i__ + en * h_dim1]) / x;
622 h__[i__ + 1 + en * h_dim1] = (-sa - w * h__[i__ + en * h_dim1] -
623 q * h__[i__ + na * h_dim1]) / x;
624 goto L790;
625 L785:
626 d__1 = -r__ - y * h__[i__ + na * h_dim1];
627 d__2 = -s - y * h__[i__ + en * h_dim1];
628 cdiv_(&d__1, &d__2, &zz, &q, &h__[i__ + 1 + na * h_dim1], &h__[
629 i__ + 1 + en * h_dim1]);
630
631 /* .......... OVERFLOW CONTROL .......... */
632 L790:
633 /* Computing MAX */
634 d__3 = (d__1 = h__[i__ + na * h_dim1], abs(d__1)), d__4 = (d__2 =
635 h__[i__ + en * h_dim1], abs(d__2));
636 t = max(d__3,d__4);
637 if (t == 0.) {
638 goto L795;
639 }
640 tst1 = t;
641 tst2 = tst1 + 1. / tst1;
642 if (tst2 > tst1) {
643 goto L795;
644 }
645 i__3 = en;
646 for (j = i__; j <= i__3; ++j) {
647 h__[j + na * h_dim1] /= t;
648 h__[j + en * h_dim1] /= t;
649 /* L792: */
650 }
651
652 L795:
653 ;
654 }
655 /* .......... END COMPLEX VECTOR .......... */
656 L800:
657 ;
658 }
659 /* .......... END BACK SUBSTITUTION. */
660 /* VECTORS OF ISOLATED ROOTS .......... */
661 i__1 = *n;
662 for (i__ = 1; i__ <= i__1; ++i__) {
663 if (i__ >= *low && i__ <= *igh) {
664 goto L840;
665 }
666
667 i__2 = *n;
668 for (j = i__; j <= i__2; ++j) {
669 /* L820: */
670 z__[i__ + j * z_dim1] = h__[i__ + j * h_dim1];
671 }
672
673 L840:
674 ;
675 }
676 /* .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
677 /* VECTORS OF ORIGINAL FULL MATRIX. */
678 /* FOR J=N STEP -1 UNTIL LOW DO -- .......... */
679 i__1 = *n;
680 for (jj = *low; jj <= i__1; ++jj) {
681 j = *n + *low - jj;
682 m = min(j,*igh);
683
684 i__2 = *igh;
685 for (i__ = *low; i__ <= i__2; ++i__) {
686 zz = 0.;
687
688 i__3 = m;
689 for (k = *low; k <= i__3; ++k) {
690 /* L860: */
691 zz += z__[i__ + k * z_dim1] * h__[k + j * h_dim1];
692 }
693
694 z__[i__ + j * z_dim1] = zz;
695 /* L880: */
696 }
697 }
698
699 goto L1001;
700 /* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
701 /* CONVERGED AFTER 30*N ITERATIONS .......... */
702 L1000:
703 *ierr = en;
704 L1001:
705 return 0;
706 } /* hqr2_ */
707
708