1 /* eispack/hqr2.f -- translated by f2c (version 20050501).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /* Table of constant values */
19
20 static doublereal c_b49 = 0.;
21
22 /*< subroutine hqr2(nm,n,low,igh,h,wr,wi,z,ierr) >*/
hqr2_(integer * nm,integer * n,integer * low,integer * igh,doublereal * h__,doublereal * wr,doublereal * wi,doublereal * z__,integer * ierr)23 /* Subroutine */ int hqr2_(integer *nm, integer *n, integer *low, integer *
24 igh, doublereal *h__, doublereal *wr, doublereal *wi, doublereal *z__,
25 integer *ierr)
26 {
27 /* System generated locals */
28 integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3;
29 doublereal d__1, d__2, d__3, d__4;
30
31 /* Builtin functions */
32 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
33
34 /* Local variables */
35 integer i__, j, k, l=0, m=0;
36 doublereal p, q, r__=0, s=0, t, w, x, y;
37 integer na, ii, en, jj;
38 doublereal ra, sa;
39 integer ll, mm, nn;
40 doublereal vi, vr, zz;
41 integer mp2, itn, its, enm2;
42 doublereal tst1, tst2;
43 extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
44 , doublereal *, doublereal *, doublereal *);
45 doublereal norm;
46 logical notlas;
47
48
49 /*< >*/
50 /*< double precision h(nm,n),wr(n),wi(n),z(nm,n) >*/
51 /*< double precision p,q,r,s,t,w,x,y,ra,sa,vi,vr,zz,norm,tst1,tst2 >*/
52 /*< logical notlas >*/
53
54 /* this subroutine is a translation of the algol procedure hqr2, */
55 /* num. math. 16, 181-204(1970) by peters and wilkinson. */
56 /* handbook for auto. comp., vol.ii-linear algebra, 372-395(1971). */
57
58 /* this subroutine finds the eigenvalues and eigenvectors */
59 /* of a real upper hessenberg matrix by the qr method. the */
60 /* eigenvectors of a real general matrix can also be found */
61 /* if elmhes and eltran or orthes and ortran have */
62 /* been used to reduce this general matrix to hessenberg form */
63 /* and to accumulate the similarity transformations. */
64
65 /* on input */
66
67 /* nm must be set to the row dimension of two-dimensional */
68 /* array parameters as declared in the calling program */
69 /* dimension statement. */
70
71 /* n is the order of the matrix. */
72
73 /* low and igh are integers determined by the balancing */
74 /* subroutine balanc. if balanc has not been used, */
75 /* set low=1, igh=n. */
76
77 /* h contains the upper hessenberg matrix. */
78
79 /* z contains the transformation matrix produced by eltran */
80 /* after the reduction by elmhes, or by ortran after the */
81 /* reduction by orthes, if performed. if the eigenvectors */
82 /* of the hessenberg matrix are desired, z must contain the */
83 /* identity matrix. */
84
85 /* on output */
86
87 /* h has been destroyed. */
88
89 /* wr and wi contain the real and imaginary parts, */
90 /* respectively, of the eigenvalues. the eigenvalues */
91 /* are unordered except that complex conjugate pairs */
92 /* of values appear consecutively with the eigenvalue */
93 /* having the positive imaginary part first. if an */
94 /* error exit is made, the eigenvalues should be correct */
95 /* for indices ierr+1,...,n. */
96
97 /* z contains the real and imaginary parts of the eigenvectors. */
98 /* if the i-th eigenvalue is real, the i-th column of z */
99 /* contains its eigenvector. if the i-th eigenvalue is complex */
100 /* with positive imaginary part, the i-th and (i+1)-th */
101 /* columns of z contain the real and imaginary parts of its */
102 /* eigenvector. the eigenvectors are unnormalized. if an */
103 /* error exit is made, none of the eigenvectors has been found. */
104
105 /* ierr is set to */
106 /* zero for normal return, */
107 /* j if the limit of 30*n iterations is exhausted */
108 /* while the j-th eigenvalue is being sought. */
109
110 /* calls cdiv for complex division. */
111
112 /* questions and comments should be directed to burton s. garbow, */
113 /* mathematics and computer science div, argonne national laboratory */
114
115 /* this version dated august 1983. */
116
117 /* ------------------------------------------------------------------ */
118
119 /*< ierr = 0 >*/
120 /* Parameter adjustments */
121 z_dim1 = *nm;
122 z_offset = 1 + z_dim1;
123 z__ -= z_offset;
124 --wi;
125 --wr;
126 h_dim1 = *nm;
127 h_offset = 1 + h_dim1;
128 h__ -= h_offset;
129
130 /* Function Body */
131 *ierr = 0;
132 /*< norm = 0.0d0 >*/
133 norm = 0.;
134 /*< k = 1 >*/
135 k = 1;
136 /* .......... store roots isolated by balanc */
137 /* and compute matrix norm .......... */
138 /*< do 50 i = 1, n >*/
139 i__1 = *n;
140 for (i__ = 1; i__ <= i__1; ++i__) {
141
142 /*< do 40 j = k, n >*/
143 i__2 = *n;
144 for (j = k; j <= i__2; ++j) {
145 /*< 40 norm = norm + dabs(h(i,j)) >*/
146 /* L40: */
147 norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1));
148 }
149
150 /*< k = i >*/
151 k = i__;
152 /*< if (i .ge. low .and. i .le. igh) go to 50 >*/
153 if (i__ >= *low && i__ <= *igh) {
154 goto L50;
155 }
156 /*< wr(i) = h(i,i) >*/
157 wr[i__] = h__[i__ + i__ * h_dim1];
158 /*< wi(i) = 0.0d0 >*/
159 wi[i__] = 0.;
160 /*< 50 continue >*/
161 L50:
162 ;
163 }
164
165 /*< en = igh >*/
166 en = *igh;
167 /*< t = 0.0d0 >*/
168 t = 0.;
169 /*< itn = 30*n >*/
170 itn = *n * 30;
171 /* .......... search for next eigenvalues .......... */
172 /*< 60 if (en .lt. low) go to 340 >*/
173 L60:
174 if (en < *low) {
175 goto L340;
176 }
177 /*< its = 0 >*/
178 its = 0;
179 /*< na = en - 1 >*/
180 na = en - 1;
181 /*< enm2 = na - 1 >*/
182 enm2 = na - 1;
183 /* .......... look for single small sub-diagonal element */
184 /* for l=en step -1 until low do -- .......... */
185 /*< 70 do 80 ll = low, en >*/
186 L70:
187 i__1 = en;
188 for (ll = *low; ll <= i__1; ++ll) {
189 /*< l = en + low - ll >*/
190 l = en + *low - ll;
191 /*< if (l .eq. low) go to 100 >*/
192 if (l == *low) {
193 goto L100;
194 }
195 /*< s = dabs(h(l-1,l-1)) + dabs(h(l,l)) >*/
196 s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l
197 + l * h_dim1], abs(d__2));
198 /*< if (s .eq. 0.0d0) s = norm >*/
199 if (s == 0.) {
200 s = norm;
201 }
202 /*< tst1 = s >*/
203 tst1 = s;
204 /*< tst2 = tst1 + dabs(h(l,l-1)) >*/
205 tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1));
206 /*< if (tst2 .eq. tst1) go to 100 >*/
207 if (tst2 == tst1) {
208 goto L100;
209 }
210 /*< 80 continue >*/
211 /* L80: */
212 }
213 /* .......... form shift .......... */
214 /*< 100 x = h(en,en) >*/
215 L100:
216 x = h__[en + en * h_dim1];
217 /*< if (l .eq. en) go to 270 >*/
218 if (l == en) {
219 goto L270;
220 }
221 /*< y = h(na,na) >*/
222 y = h__[na + na * h_dim1];
223 /*< w = h(en,na) * h(na,en) >*/
224 w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
225 /*< if (l .eq. na) go to 280 >*/
226 if (l == na) {
227 goto L280;
228 }
229 /*< if (itn .eq. 0) go to 1000 >*/
230 if (itn == 0) {
231 goto L1000;
232 }
233 /*< if (its .ne. 10 .and. its .ne. 20) go to 130 >*/
234 if (its != 10 && its != 20) {
235 goto L130;
236 }
237 /* .......... form exceptional shift .......... */
238 /*< t = t + x >*/
239 t += x;
240
241 /*< do 120 i = low, en >*/
242 i__1 = en;
243 for (i__ = *low; i__ <= i__1; ++i__) {
244 /*< 120 h(i,i) = h(i,i) - x >*/
245 /* L120: */
246 h__[i__ + i__ * h_dim1] -= x;
247 }
248
249 /*< s = dabs(h(en,na)) + dabs(h(na,enm2)) >*/
250 s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 *
251 h_dim1], abs(d__2));
252 /*< x = 0.75d0 * s >*/
253 x = s * .75;
254 /*< y = x >*/
255 y = x;
256 /*< w = -0.4375d0 * s * s >*/
257 w = s * -.4375 * s;
258 /*< 130 its = its + 1 >*/
259 L130:
260 ++its;
261 /*< itn = itn - 1 >*/
262 --itn;
263 /* .......... look for two consecutive small */
264 /* sub-diagonal elements. */
265 /* for m=en-2 step -1 until l do -- .......... */
266 /*< do 140 mm = l, enm2 >*/
267 i__1 = enm2;
268 for (mm = l; mm <= i__1; ++mm) {
269 /*< m = enm2 + l - mm >*/
270 m = enm2 + l - mm;
271 /*< zz = h(m,m) >*/
272 zz = h__[m + m * h_dim1];
273 /*< r = x - zz >*/
274 r__ = x - zz;
275 /*< s = y - zz >*/
276 s = y - zz;
277 /*< p = (r * s - w) / h(m+1,m) + h(m,m+1) >*/
278 p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) *
279 h_dim1];
280 /*< q = h(m+1,m+1) - zz - r - s >*/
281 q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
282 /*< r = h(m+2,m+1) >*/
283 r__ = h__[m + 2 + (m + 1) * h_dim1];
284 /*< s = dabs(p) + dabs(q) + dabs(r) >*/
285 s = abs(p) + abs(q) + abs(r__);
286 /*< p = p / s >*/
287 p /= s;
288 /*< q = q / s >*/
289 q /= s;
290 /*< r = r / s >*/
291 r__ /= s;
292 /*< if (m .eq. l) go to 150 >*/
293 if (m == l) {
294 goto L150;
295 }
296 /*< tst1 = dabs(p)*(dabs(h(m-1,m-1)) + dabs(zz) + dabs(h(m+1,m+1))) >*/
297 tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) +
298 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
299 /*< tst2 = tst1 + dabs(h(m,m-1))*(dabs(q) + dabs(r)) >*/
300 tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q)
301 + abs(r__));
302 /*< if (tst2 .eq. tst1) go to 150 >*/
303 if (tst2 == tst1) {
304 goto L150;
305 }
306 /*< 140 continue >*/
307 /* L140: */
308 }
309
310 /*< 150 mp2 = m + 2 >*/
311 L150:
312 mp2 = m + 2;
313
314 /*< do 160 i = mp2, en >*/
315 i__1 = en;
316 for (i__ = mp2; i__ <= i__1; ++i__) {
317 /*< h(i,i-2) = 0.0d0 >*/
318 h__[i__ + (i__ - 2) * h_dim1] = 0.;
319 /*< if (i .eq. mp2) go to 160 >*/
320 if (i__ == mp2) {
321 goto L160;
322 }
323 /*< h(i,i-3) = 0.0d0 >*/
324 h__[i__ + (i__ - 3) * h_dim1] = 0.;
325 /*< 160 continue >*/
326 L160:
327 ;
328 }
329 /* .......... double qr step involving rows l to en and */
330 /* columns m to en .......... */
331 /*< do 260 k = m, na >*/
332 i__1 = na;
333 for (k = m; k <= i__1; ++k) {
334 /*< notlas = k .ne. na >*/
335 notlas = k != na;
336 /*< if (k .eq. m) go to 170 >*/
337 if (k == m) {
338 goto L170;
339 }
340 /*< p = h(k,k-1) >*/
341 p = h__[k + (k - 1) * h_dim1];
342 /*< q = h(k+1,k-1) >*/
343 q = h__[k + 1 + (k - 1) * h_dim1];
344 /*< r = 0.0d0 >*/
345 r__ = 0.;
346 /*< if (notlas) r = h(k+2,k-1) >*/
347 if (notlas) {
348 r__ = h__[k + 2 + (k - 1) * h_dim1];
349 }
350 /*< x = dabs(p) + dabs(q) + dabs(r) >*/
351 x = abs(p) + abs(q) + abs(r__);
352 /*< if (x .eq. 0.0d0) go to 260 >*/
353 if (x == 0.) {
354 goto L260;
355 }
356 /*< p = p / x >*/
357 p /= x;
358 /*< q = q / x >*/
359 q /= x;
360 /*< r = r / x >*/
361 r__ /= x;
362 /*< 170 s = dsign(dsqrt(p*p+q*q+r*r),p) >*/
363 L170:
364 d__1 = sqrt(p * p + q * q + r__ * r__);
365 s = d_sign(&d__1, &p);
366 /*< if (k .eq. m) go to 180 >*/
367 if (k == m) {
368 goto L180;
369 }
370 /*< h(k,k-1) = -s * x >*/
371 h__[k + (k - 1) * h_dim1] = -s * x;
372 /*< go to 190 >*/
373 goto L190;
374 /*< 180 if (l .ne. m) h(k,k-1) = -h(k,k-1) >*/
375 L180:
376 if (l != m) {
377 h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
378 }
379 /*< 190 p = p + s >*/
380 L190:
381 p += s;
382 /*< x = p / s >*/
383 x = p / s;
384 /*< y = q / s >*/
385 y = q / s;
386 /*< zz = r / s >*/
387 zz = r__ / s;
388 /*< q = q / p >*/
389 q /= p;
390 /*< r = r / p >*/
391 r__ /= p;
392 /*< if (notlas) go to 225 >*/
393 if (notlas) {
394 goto L225;
395 }
396 /* .......... row modification .......... */
397 /*< do 200 j = k, n >*/
398 i__2 = *n;
399 for (j = k; j <= i__2; ++j) {
400 /*< p = h(k,j) + q * h(k+1,j) >*/
401 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
402 /*< h(k,j) = h(k,j) - p * x >*/
403 h__[k + j * h_dim1] -= p * x;
404 /*< h(k+1,j) = h(k+1,j) - p * y >*/
405 h__[k + 1 + j * h_dim1] -= p * y;
406 /*< 200 continue >*/
407 /* L200: */
408 }
409
410 /*< j = min0(en,k+3) >*/
411 /* Computing MIN */
412 i__2 = en, i__3 = k + 3;
413 j = min(i__2,i__3);
414 /* .......... column modification .......... */
415 /*< do 210 i = 1, j >*/
416 i__2 = j;
417 for (i__ = 1; i__ <= i__2; ++i__) {
418 /*< p = x * h(i,k) + y * h(i,k+1) >*/
419 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
420 /*< h(i,k) = h(i,k) - p >*/
421 h__[i__ + k * h_dim1] -= p;
422 /*< h(i,k+1) = h(i,k+1) - p * q >*/
423 h__[i__ + (k + 1) * h_dim1] -= p * q;
424 /*< 210 continue >*/
425 /* L210: */
426 }
427 /* .......... accumulate transformations .......... */
428 /*< do 220 i = low, igh >*/
429 i__2 = *igh;
430 for (i__ = *low; i__ <= i__2; ++i__) {
431 /*< p = x * z(i,k) + y * z(i,k+1) >*/
432 p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1];
433 /*< z(i,k) = z(i,k) - p >*/
434 z__[i__ + k * z_dim1] -= p;
435 /*< z(i,k+1) = z(i,k+1) - p * q >*/
436 z__[i__ + (k + 1) * z_dim1] -= p * q;
437 /*< 220 continue >*/
438 /* L220: */
439 }
440 /*< go to 255 >*/
441 goto L255;
442 /*< 225 continue >*/
443 L225:
444 /* .......... row modification .......... */
445 /*< do 230 j = k, n >*/
446 i__2 = *n;
447 for (j = k; j <= i__2; ++j) {
448 /*< p = h(k,j) + q * h(k+1,j) + r * h(k+2,j) >*/
449 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
450 k + 2 + j * h_dim1];
451 /*< h(k,j) = h(k,j) - p * x >*/
452 h__[k + j * h_dim1] -= p * x;
453 /*< h(k+1,j) = h(k+1,j) - p * y >*/
454 h__[k + 1 + j * h_dim1] -= p * y;
455 /*< h(k+2,j) = h(k+2,j) - p * zz >*/
456 h__[k + 2 + j * h_dim1] -= p * zz;
457 /*< 230 continue >*/
458 /* L230: */
459 }
460
461 /*< j = min0(en,k+3) >*/
462 /* Computing MIN */
463 i__2 = en, i__3 = k + 3;
464 j = min(i__2,i__3);
465 /* .......... column modification .......... */
466 /*< do 240 i = 1, j >*/
467 i__2 = j;
468 for (i__ = 1; i__ <= i__2; ++i__) {
469 /*< p = x * h(i,k) + y * h(i,k+1) + zz * h(i,k+2) >*/
470 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] +
471 zz * h__[i__ + (k + 2) * h_dim1];
472 /*< h(i,k) = h(i,k) - p >*/
473 h__[i__ + k * h_dim1] -= p;
474 /*< h(i,k+1) = h(i,k+1) - p * q >*/
475 h__[i__ + (k + 1) * h_dim1] -= p * q;
476 /*< h(i,k+2) = h(i,k+2) - p * r >*/
477 h__[i__ + (k + 2) * h_dim1] -= p * r__;
478 /*< 240 continue >*/
479 /* L240: */
480 }
481 /* .......... accumulate transformations .......... */
482 /*< do 250 i = low, igh >*/
483 i__2 = *igh;
484 for (i__ = *low; i__ <= i__2; ++i__) {
485 /*< p = x * z(i,k) + y * z(i,k+1) + zz * z(i,k+2) >*/
486 p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1] +
487 zz * z__[i__ + (k + 2) * z_dim1];
488 /*< z(i,k) = z(i,k) - p >*/
489 z__[i__ + k * z_dim1] -= p;
490 /*< z(i,k+1) = z(i,k+1) - p * q >*/
491 z__[i__ + (k + 1) * z_dim1] -= p * q;
492 /*< z(i,k+2) = z(i,k+2) - p * r >*/
493 z__[i__ + (k + 2) * z_dim1] -= p * r__;
494 /*< 250 continue >*/
495 /* L250: */
496 }
497 /*< 255 continue >*/
498 L255:
499
500 /*< 260 continue >*/
501 L260:
502 ;
503 }
504
505 /*< go to 70 >*/
506 goto L70;
507 /* .......... one root found .......... */
508 /*< 270 h(en,en) = x + t >*/
509 L270:
510 h__[en + en * h_dim1] = x + t;
511 /*< wr(en) = h(en,en) >*/
512 wr[en] = h__[en + en * h_dim1];
513 /*< wi(en) = 0.0d0 >*/
514 wi[en] = 0.;
515 /*< en = na >*/
516 en = na;
517 /*< go to 60 >*/
518 goto L60;
519 /* .......... two roots found .......... */
520 /*< 280 p = (y - x) / 2.0d0 >*/
521 L280:
522 p = (y - x) / 2.;
523 /*< q = p * p + w >*/
524 q = p * p + w;
525 /*< zz = dsqrt(dabs(q)) >*/
526 zz = sqrt((abs(q)));
527 /*< h(en,en) = x + t >*/
528 h__[en + en * h_dim1] = x + t;
529 /*< x = h(en,en) >*/
530 x = h__[en + en * h_dim1];
531 /*< h(na,na) = y + t >*/
532 h__[na + na * h_dim1] = y + t;
533 /*< if (q .lt. 0.0d0) go to 320 >*/
534 if (q < 0.) {
535 goto L320;
536 }
537 /* .......... real pair .......... */
538 /*< zz = p + dsign(zz,p) >*/
539 zz = p + d_sign(&zz, &p);
540 /*< wr(na) = x + zz >*/
541 wr[na] = x + zz;
542 /*< wr(en) = wr(na) >*/
543 wr[en] = wr[na];
544 /*< if (zz .ne. 0.0d0) wr(en) = x - w / zz >*/
545 if (zz != 0.) {
546 wr[en] = x - w / zz;
547 }
548 /*< wi(na) = 0.0d0 >*/
549 wi[na] = 0.;
550 /*< wi(en) = 0.0d0 >*/
551 wi[en] = 0.;
552 /*< x = h(en,na) >*/
553 x = h__[en + na * h_dim1];
554 /*< s = dabs(x) + dabs(zz) >*/
555 s = abs(x) + abs(zz);
556 /*< p = x / s >*/
557 p = x / s;
558 /*< q = zz / s >*/
559 q = zz / s;
560 /*< r = dsqrt(p*p+q*q) >*/
561 r__ = sqrt(p * p + q * q);
562 /*< p = p / r >*/
563 p /= r__;
564 /*< q = q / r >*/
565 q /= r__;
566 /* .......... row modification .......... */
567 /*< do 290 j = na, n >*/
568 i__1 = *n;
569 for (j = na; j <= i__1; ++j) {
570 /*< zz = h(na,j) >*/
571 zz = h__[na + j * h_dim1];
572 /*< h(na,j) = q * zz + p * h(en,j) >*/
573 h__[na + j * h_dim1] = q * zz + p * h__[en + j * h_dim1];
574 /*< h(en,j) = q * h(en,j) - p * zz >*/
575 h__[en + j * h_dim1] = q * h__[en + j * h_dim1] - p * zz;
576 /*< 290 continue >*/
577 /* L290: */
578 }
579 /* .......... column modification .......... */
580 /*< do 300 i = 1, en >*/
581 i__1 = en;
582 for (i__ = 1; i__ <= i__1; ++i__) {
583 /*< zz = h(i,na) >*/
584 zz = h__[i__ + na * h_dim1];
585 /*< h(i,na) = q * zz + p * h(i,en) >*/
586 h__[i__ + na * h_dim1] = q * zz + p * h__[i__ + en * h_dim1];
587 /*< h(i,en) = q * h(i,en) - p * zz >*/
588 h__[i__ + en * h_dim1] = q * h__[i__ + en * h_dim1] - p * zz;
589 /*< 300 continue >*/
590 /* L300: */
591 }
592 /* .......... accumulate transformations .......... */
593 /*< do 310 i = low, igh >*/
594 i__1 = *igh;
595 for (i__ = *low; i__ <= i__1; ++i__) {
596 /*< zz = z(i,na) >*/
597 zz = z__[i__ + na * z_dim1];
598 /*< z(i,na) = q * zz + p * z(i,en) >*/
599 z__[i__ + na * z_dim1] = q * zz + p * z__[i__ + en * z_dim1];
600 /*< z(i,en) = q * z(i,en) - p * zz >*/
601 z__[i__ + en * z_dim1] = q * z__[i__ + en * z_dim1] - p * zz;
602 /*< 310 continue >*/
603 /* L310: */
604 }
605
606 /*< go to 330 >*/
607 goto L330;
608 /* .......... complex pair .......... */
609 /*< 320 wr(na) = x + p >*/
610 L320:
611 wr[na] = x + p;
612 /*< wr(en) = x + p >*/
613 wr[en] = x + p;
614 /*< wi(na) = zz >*/
615 wi[na] = zz;
616 /*< wi(en) = -zz >*/
617 wi[en] = -zz;
618 /*< 330 en = enm2 >*/
619 L330:
620 en = enm2;
621 /*< go to 60 >*/
622 goto L60;
623 /* .......... all roots found. backsubstitute to find */
624 /* vectors of upper triangular form .......... */
625 /*< 340 if (norm .eq. 0.0d0) go to 1001 >*/
626 L340:
627 if (norm == 0.) {
628 goto L1001;
629 }
630 /* .......... for en=n step -1 until 1 do -- .......... */
631 /*< do 800 nn = 1, n >*/
632 i__1 = *n;
633 for (nn = 1; nn <= i__1; ++nn) {
634 /*< en = n + 1 - nn >*/
635 en = *n + 1 - nn;
636 /*< p = wr(en) >*/
637 p = wr[en];
638 /*< q = wi(en) >*/
639 q = wi[en];
640 /*< na = en - 1 >*/
641 na = en - 1;
642 /*< if (q) 710, 600, 800 >*/
643 if (q < 0.) {
644 goto L710;
645 } else if (q == 0) {
646 goto L600;
647 } else {
648 goto L800;
649 }
650 /* .......... real vector .......... */
651 /*< 600 m = en >*/
652 L600:
653 m = en;
654 /*< h(en,en) = 1.0d0 >*/
655 h__[en + en * h_dim1] = 1.;
656 /*< if (na .eq. 0) go to 800 >*/
657 if (na == 0) {
658 goto L800;
659 }
660 /* .......... for i=en-1 step -1 until 1 do -- .......... */
661 /*< do 700 ii = 1, na >*/
662 i__2 = na;
663 for (ii = 1; ii <= i__2; ++ii) {
664 /*< i = en - ii >*/
665 i__ = en - ii;
666 /*< w = h(i,i) - p >*/
667 w = h__[i__ + i__ * h_dim1] - p;
668 /*< r = 0.0d0 >*/
669 r__ = 0.;
670
671 /*< do 610 j = m, en >*/
672 i__3 = en;
673 for (j = m; j <= i__3; ++j) {
674 /*< 610 r = r + h(i,j) * h(j,en) >*/
675 /* L610: */
676 r__ += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
677 }
678
679 /*< if (wi(i) .ge. 0.0d0) go to 630 >*/
680 if (wi[i__] >= 0.) {
681 goto L630;
682 }
683 /*< zz = w >*/
684 zz = w;
685 /*< s = r >*/
686 s = r__;
687 /*< go to 700 >*/
688 goto L700;
689 /*< 630 m = i >*/
690 L630:
691 m = i__;
692 /*< if (wi(i) .ne. 0.0d0) go to 640 >*/
693 if (wi[i__] != 0.) {
694 goto L640;
695 }
696 /*< t = w >*/
697 t = w;
698 /*< if (t .ne. 0.0d0) go to 635 >*/
699 if (t != 0.) {
700 goto L635;
701 }
702 /*< tst1 = norm >*/
703 tst1 = norm;
704 /*< t = tst1 >*/
705 t = tst1;
706 /*< 632 t = 0.01d0 * t >*/
707 L632:
708 t *= .01;
709 /*< tst2 = norm + t >*/
710 tst2 = norm + t;
711 /*< if (tst2 .gt. tst1) go to 632 >*/
712 if (tst2 > tst1) {
713 goto L632;
714 }
715 /*< 635 h(i,en) = -r / t >*/
716 L635:
717 h__[i__ + en * h_dim1] = -r__ / t;
718 /*< go to 680 >*/
719 goto L680;
720 /* .......... solve real equations .......... */
721 /*< 640 x = h(i,i+1) >*/
722 L640:
723 x = h__[i__ + (i__ + 1) * h_dim1];
724 /*< y = h(i+1,i) >*/
725 y = h__[i__ + 1 + i__ * h_dim1];
726 /*< q = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) >*/
727 q = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__];
728 /*< t = (x * s - zz * r) / q >*/
729 t = (x * s - zz * r__) / q;
730 /*< h(i,en) = t >*/
731 h__[i__ + en * h_dim1] = t;
732 /*< if (dabs(x) .le. dabs(zz)) go to 650 >*/
733 if (abs(x) <= abs(zz)) {
734 goto L650;
735 }
736 /*< h(i+1,en) = (-r - w * t) / x >*/
737 h__[i__ + 1 + en * h_dim1] = (-r__ - w * t) / x;
738 /*< go to 680 >*/
739 goto L680;
740 /*< 650 h(i+1,en) = (-s - y * t) / zz >*/
741 L650:
742 h__[i__ + 1 + en * h_dim1] = (-s - y * t) / zz;
743
744 /* .......... overflow control .......... */
745 /*< 680 t = dabs(h(i,en)) >*/
746 L680:
747 t = (d__1 = h__[i__ + en * h_dim1], abs(d__1));
748 /*< if (t .eq. 0.0d0) go to 700 >*/
749 if (t == 0.) {
750 goto L700;
751 }
752 /*< tst1 = t >*/
753 tst1 = t;
754 /*< tst2 = tst1 + 1.0d0/tst1 >*/
755 tst2 = tst1 + 1. / tst1;
756 /*< if (tst2 .gt. tst1) go to 700 >*/
757 if (tst2 > tst1) {
758 goto L700;
759 }
760 /*< do 690 j = i, en >*/
761 i__3 = en;
762 for (j = i__; j <= i__3; ++j) {
763 /*< h(j,en) = h(j,en)/t >*/
764 h__[j + en * h_dim1] /= t;
765 /*< 690 continue >*/
766 /* L690: */
767 }
768
769 /*< 700 continue >*/
770 L700:
771 ;
772 }
773 /* .......... end real vector .......... */
774 /*< go to 800 >*/
775 goto L800;
776 /* .......... complex vector .......... */
777 /*< 710 m = na >*/
778 L710:
779 m = na;
780 /* .......... last vector component chosen imaginary so that */
781 /* eigenvector matrix is triangular .......... */
782 /*< if (dabs(h(en,na)) .le. dabs(h(na,en))) go to 720 >*/
783 if ((d__1 = h__[en + na * h_dim1], abs(d__1)) <= (d__2 = h__[na + en *
784 h_dim1], abs(d__2))) {
785 goto L720;
786 }
787 /*< h(na,na) = q / h(en,na) >*/
788 h__[na + na * h_dim1] = q / h__[en + na * h_dim1];
789 /*< h(na,en) = -(h(en,en) - p) / h(en,na) >*/
790 h__[na + en * h_dim1] = -(h__[en + en * h_dim1] - p) / h__[en + na *
791 h_dim1];
792 /*< go to 730 >*/
793 goto L730;
794 /*< 720 call cdiv(0.0d0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en)) >*/
795 L720:
796 d__1 = -h__[na + en * h_dim1];
797 d__2 = h__[na + na * h_dim1] - p;
798 cdiv_(&c_b49, &d__1, &d__2, &q, &h__[na + na * h_dim1], &h__[na + en *
799 h_dim1]);
800 /*< 730 h(en,na) = 0.0d0 >*/
801 L730:
802 h__[en + na * h_dim1] = 0.;
803 /*< h(en,en) = 1.0d0 >*/
804 h__[en + en * h_dim1] = 1.;
805 /*< enm2 = na - 1 >*/
806 enm2 = na - 1;
807 /*< if (enm2 .eq. 0) go to 800 >*/
808 if (enm2 == 0) {
809 goto L800;
810 }
811 /* .......... for i=en-2 step -1 until 1 do -- .......... */
812 /*< do 795 ii = 1, enm2 >*/
813 i__2 = enm2;
814 for (ii = 1; ii <= i__2; ++ii) {
815 /*< i = na - ii >*/
816 i__ = na - ii;
817 /*< w = h(i,i) - p >*/
818 w = h__[i__ + i__ * h_dim1] - p;
819 /*< ra = 0.0d0 >*/
820 ra = 0.;
821 /*< sa = 0.0d0 >*/
822 sa = 0.;
823
824 /*< do 760 j = m, en >*/
825 i__3 = en;
826 for (j = m; j <= i__3; ++j) {
827 /*< ra = ra + h(i,j) * h(j,na) >*/
828 ra += h__[i__ + j * h_dim1] * h__[j + na * h_dim1];
829 /*< sa = sa + h(i,j) * h(j,en) >*/
830 sa += h__[i__ + j * h_dim1] * h__[j + en * h_dim1];
831 /*< 760 continue >*/
832 /* L760: */
833 }
834
835 /*< if (wi(i) .ge. 0.0d0) go to 770 >*/
836 if (wi[i__] >= 0.) {
837 goto L770;
838 }
839 /*< zz = w >*/
840 zz = w;
841 /*< r = ra >*/
842 r__ = ra;
843 /*< s = sa >*/
844 s = sa;
845 /*< go to 795 >*/
846 goto L795;
847 /*< 770 m = i >*/
848 L770:
849 m = i__;
850 /*< if (wi(i) .ne. 0.0d0) go to 780 >*/
851 if (wi[i__] != 0.) {
852 goto L780;
853 }
854 /*< call cdiv(-ra,-sa,w,q,h(i,na),h(i,en)) >*/
855 d__1 = -ra;
856 d__2 = -sa;
857 cdiv_(&d__1, &d__2, &w, &q, &h__[i__ + na * h_dim1], &h__[i__ +
858 en * h_dim1]);
859 /*< go to 790 >*/
860 goto L790;
861 /* .......... solve complex equations .......... */
862 /*< 780 x = h(i,i+1) >*/
863 L780:
864 x = h__[i__ + (i__ + 1) * h_dim1];
865 /*< y = h(i+1,i) >*/
866 y = h__[i__ + 1 + i__ * h_dim1];
867 /*< vr = (wr(i) - p) * (wr(i) - p) + wi(i) * wi(i) - q * q >*/
868 vr = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__] - q * q;
869 /*< vi = (wr(i) - p) * 2.0d0 * q >*/
870 vi = (wr[i__] - p) * 2. * q;
871 /*< if (vr .ne. 0.0d0 .or. vi .ne. 0.0d0) go to 784 >*/
872 if (vr != 0. || vi != 0.) {
873 goto L784;
874 }
875 /*< >*/
876 tst1 = norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(zz));
877 /*< vr = tst1 >*/
878 vr = tst1;
879 /*< 783 vr = 0.01d0 * vr >*/
880 L783:
881 vr *= .01;
882 /*< tst2 = tst1 + vr >*/
883 tst2 = tst1 + vr;
884 /*< if (tst2 .gt. tst1) go to 783 >*/
885 if (tst2 > tst1) {
886 goto L783;
887 }
888 /*< >*/
889 L784:
890 d__1 = x * r__ - zz * ra + q * sa;
891 d__2 = x * s - zz * sa - q * ra;
892 cdiv_(&d__1, &d__2, &vr, &vi, &h__[i__ + na * h_dim1], &h__[i__ +
893 en * h_dim1]);
894 /*< if (dabs(x) .le. dabs(zz) + dabs(q)) go to 785 >*/
895 if (abs(x) <= abs(zz) + abs(q)) {
896 goto L785;
897 }
898 /*< h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x >*/
899 h__[i__ + 1 + na * h_dim1] = (-ra - w * h__[i__ + na * h_dim1] +
900 q * h__[i__ + en * h_dim1]) / x;
901 /*< h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x >*/
902 h__[i__ + 1 + en * h_dim1] = (-sa - w * h__[i__ + en * h_dim1] -
903 q * h__[i__ + na * h_dim1]) / x;
904 /*< go to 790 >*/
905 goto L790;
906 /*< >*/
907 L785:
908 d__1 = -r__ - y * h__[i__ + na * h_dim1];
909 d__2 = -s - y * h__[i__ + en * h_dim1];
910 cdiv_(&d__1, &d__2, &zz, &q, &h__[i__ + 1 + na * h_dim1], &h__[
911 i__ + 1 + en * h_dim1]);
912
913 /* .......... overflow control .......... */
914 /*< 790 t = dmax1(dabs(h(i,na)), dabs(h(i,en))) >*/
915 L790:
916 /* Computing MAX */
917 d__3 = (d__1 = h__[i__ + na * h_dim1], abs(d__1)), d__4 = (d__2 =
918 h__[i__ + en * h_dim1], abs(d__2));
919 t = max(d__3,d__4);
920 /*< if (t .eq. 0.0d0) go to 795 >*/
921 if (t == 0.) {
922 goto L795;
923 }
924 /*< tst1 = t >*/
925 tst1 = t;
926 /*< tst2 = tst1 + 1.0d0/tst1 >*/
927 tst2 = tst1 + 1. / tst1;
928 /*< if (tst2 .gt. tst1) go to 795 >*/
929 if (tst2 > tst1) {
930 goto L795;
931 }
932 /*< do 792 j = i, en >*/
933 i__3 = en;
934 for (j = i__; j <= i__3; ++j) {
935 /*< h(j,na) = h(j,na)/t >*/
936 h__[j + na * h_dim1] /= t;
937 /*< h(j,en) = h(j,en)/t >*/
938 h__[j + en * h_dim1] /= t;
939 /*< 792 continue >*/
940 /* L792: */
941 }
942
943 /*< 795 continue >*/
944 L795:
945 ;
946 }
947 /* .......... end complex vector .......... */
948 /*< 800 continue >*/
949 L800:
950 ;
951 }
952 /* .......... end back substitution. */
953 /* vectors of isolated roots .......... */
954 /*< do 840 i = 1, n >*/
955 i__1 = *n;
956 for (i__ = 1; i__ <= i__1; ++i__) {
957 /*< if (i .ge. low .and. i .le. igh) go to 840 >*/
958 if (i__ >= *low && i__ <= *igh) {
959 goto L840;
960 }
961
962 /*< do 820 j = i, n >*/
963 i__2 = *n;
964 for (j = i__; j <= i__2; ++j) {
965 /*< 820 z(i,j) = h(i,j) >*/
966 /* L820: */
967 z__[i__ + j * z_dim1] = h__[i__ + j * h_dim1];
968 }
969
970 /*< 840 continue >*/
971 L840:
972 ;
973 }
974 /* .......... multiply by transformation matrix to give */
975 /* vectors of original full matrix. */
976 /* for j=n step -1 until low do -- .......... */
977 /*< do 880 jj = low, n >*/
978 i__1 = *n;
979 for (jj = *low; jj <= i__1; ++jj) {
980 /*< j = n + low - jj >*/
981 j = *n + *low - jj;
982 /*< m = min0(j,igh) >*/
983 m = min(j,*igh);
984
985 /*< do 880 i = low, igh >*/
986 i__2 = *igh;
987 for (i__ = *low; i__ <= i__2; ++i__) {
988 /*< zz = 0.0d0 >*/
989 zz = 0.;
990
991 /*< do 860 k = low, m >*/
992 i__3 = m;
993 for (k = *low; k <= i__3; ++k) {
994 /*< 860 zz = zz + z(i,k) * h(k,j) >*/
995 /* L860: */
996 zz += z__[i__ + k * z_dim1] * h__[k + j * h_dim1];
997 }
998
999 /*< z(i,j) = zz >*/
1000 z__[i__ + j * z_dim1] = zz;
1001 /*< 880 continue >*/
1002 /* L880: */
1003 }
1004 }
1005
1006 /*< go to 1001 >*/
1007 goto L1001;
1008 /* .......... set error -- all eigenvalues have not */
1009 /* converged after 30*n iterations .......... */
1010 /*< 1000 ierr = en >*/
1011 L1000:
1012 *ierr = en;
1013 /*< 1001 return >*/
1014 L1001:
1015 return 0;
1016 /*< end >*/
1017 } /* hqr2_ */
1018
1019 #ifdef __cplusplus
1020 }
1021 #endif
1022