1 /* hqr.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
hqr_(integer * nm,integer * n,integer * low,integer * igh,doublereal * h__,doublereal * wr,doublereal * wi,integer * ierr)8 /* Subroutine */ int hqr_(integer *nm, integer *n, integer *low, integer *igh,
9 doublereal *h__, doublereal *wr, doublereal *wi, integer *ierr)
10 {
11 /* System generated locals */
12 integer h_dim1, h_offset, i__1, i__2, i__3;
13 doublereal d__1, d__2;
14
15 /* Builtin functions */
16 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
17
18 /* Local variables */
19 doublereal norm;
20 integer i__, j, k, l=0, m=0;
21 doublereal p, q=0.0, r__=0.0, s, t, w, x, y;
22 integer na, en, ll, mm;
23 doublereal zz;
24 logical notlas;
25 integer mp2, itn, its, enm2;
26 doublereal tst1, tst2;
27
28
29
30 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, */
31 /* NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. */
32 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). */
33
34 /* THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL */
35 /* UPPER HESSENBERG MATRIX BY THE QR METHOD. */
36
37 /* ON INPUT */
38
39 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
40 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
41 /* DIMENSION STATEMENT. */
42
43 /* N IS THE ORDER OF THE MATRIX. */
44
45 /* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
46 /* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
47 /* SET LOW=1, IGH=N. */
48
49 /* H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT */
50 /* THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG */
51 /* FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED */
52 /* IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */
53
54 /* ON OUTPUT */
55
56 /* H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED */
57 /* BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND */
58 /* BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. */
59
60 /* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
61 /* RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES */
62 /* ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */
63 /* OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */
64 /* HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN */
65 /* ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
66 /* FOR INDICES IERR+1,...,N. */
67
68 /* IERR IS SET TO */
69 /* ZERO FOR NORMAL RETURN, */
70 /* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */
71 /* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
72
73 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
74 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
75 */
76
77 /* THIS VERSION DATED AUGUST 1983. */
78
79 /* ------------------------------------------------------------------
80 */
81
82 /* Parameter adjustments */
83 --wi;
84 --wr;
85 h_dim1 = *nm;
86 h_offset = h_dim1 + 1;
87 h__ -= h_offset;
88
89 /* Function Body */
90 *ierr = 0;
91 norm = 0.;
92 k = 1;
93 /* .......... STORE ROOTS ISOLATED BY BALANC */
94 /* AND COMPUTE MATRIX NORM .......... */
95 i__1 = *n;
96 for (i__ = 1; i__ <= i__1; ++i__) {
97
98 i__2 = *n;
99 for (j = k; j <= i__2; ++j) {
100 /* L40: */
101 norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1));
102 }
103
104 k = i__;
105 if (i__ >= *low && i__ <= *igh) {
106 goto L50;
107 }
108 wr[i__] = h__[i__ + i__ * h_dim1];
109 wi[i__] = 0.;
110 L50:
111 ;
112 }
113
114 en = *igh;
115 t = 0.;
116 itn = *n * 30;
117 /* .......... SEARCH FOR NEXT EIGENVALUES .......... */
118 L60:
119 if (en < *low) {
120 goto L1001;
121 }
122 its = 0;
123 na = en - 1;
124 enm2 = na - 1;
125 /* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
126 /* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
127 L70:
128 i__1 = en;
129 for (ll = *low; ll <= i__1; ++ll) {
130 l = en + *low - ll;
131 if (l == *low) {
132 goto L100;
133 }
134 s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l
135 + l * h_dim1], abs(d__2));
136 if (s == 0.) {
137 s = norm;
138 }
139 tst1 = s;
140 tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1));
141 if (tst2 == tst1) {
142 goto L100;
143 }
144 /* L80: */
145 }
146 /* .......... FORM SHIFT .......... */
147 L100:
148 x = h__[en + en * h_dim1];
149 if (l == en) {
150 goto L270;
151 }
152 y = h__[na + na * h_dim1];
153 w = h__[en + na * h_dim1] * h__[na + en * h_dim1];
154 if (l == na) {
155 goto L280;
156 }
157 if (itn == 0) {
158 goto L1000;
159 }
160 if (its != 10 && its != 20) {
161 goto L130;
162 }
163 /* .......... FORM EXCEPTIONAL SHIFT .......... */
164 t += x;
165
166 i__1 = en;
167 for (i__ = *low; i__ <= i__1; ++i__) {
168 /* L120: */
169 h__[i__ + i__ * h_dim1] -= x;
170 }
171
172 s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 *
173 h_dim1], abs(d__2));
174 x = s * .75;
175 y = x;
176 w = s * -.4375 * s;
177 L130:
178 ++its;
179 --itn;
180 /* .......... LOOK FOR TWO CONSECUTIVE SMALL */
181 /* SUB-DIAGONAL ELEMENTS. */
182 /* FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */
183 i__1 = enm2;
184 for (mm = l; mm <= i__1; ++mm) {
185 m = enm2 + l - mm;
186 zz = h__[m + m * h_dim1];
187 r__ = x - zz;
188 s = y - zz;
189 p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) *
190 h_dim1];
191 q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s;
192 r__ = h__[m + 2 + (m + 1) * h_dim1];
193 s = abs(p) + abs(q) + abs(r__);
194 p /= s;
195 q /= s;
196 r__ /= s;
197 if (m == l) {
198 goto L150;
199 }
200 tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) +
201 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2)));
202 tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q)
203 + abs(r__));
204 if (tst2 == tst1) {
205 goto L150;
206 }
207 /* L140: */
208 }
209
210 L150:
211 mp2 = m + 2;
212
213 i__1 = en;
214 for (i__ = mp2; i__ <= i__1; ++i__) {
215 h__[i__ + (i__ - 2) * h_dim1] = 0.;
216 if (i__ == mp2) {
217 goto L160;
218 }
219 h__[i__ + (i__ - 3) * h_dim1] = 0.;
220 L160:
221 ;
222 }
223 /* .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */
224 /* COLUMNS M TO EN .......... */
225 i__1 = na;
226 for (k = m; k <= i__1; ++k) {
227 notlas = k != na;
228 if (k == m) {
229 goto L170;
230 }
231 p = h__[k + (k - 1) * h_dim1];
232 q = h__[k + 1 + (k - 1) * h_dim1];
233 r__ = 0.;
234 if (notlas) {
235 r__ = h__[k + 2 + (k - 1) * h_dim1];
236 }
237 x = abs(p) + abs(q) + abs(r__);
238 if (x == 0.) {
239 goto L260;
240 }
241 p /= x;
242 q /= x;
243 r__ /= x;
244 L170:
245 d__1 = sqrt(p * p + q * q + r__ * r__);
246 s = d_sign(&d__1, &p);
247 if (k == m) {
248 goto L180;
249 }
250 h__[k + (k - 1) * h_dim1] = -s * x;
251 goto L190;
252 L180:
253 if (l != m) {
254 h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1];
255 }
256 L190:
257 p += s;
258 x = p / s;
259 y = q / s;
260 zz = r__ / s;
261 q /= p;
262 r__ /= p;
263 if (notlas) {
264 goto L225;
265 }
266 /* .......... ROW MODIFICATION .......... */
267 i__2 = *n;
268 for (j = k; j <= i__2; ++j) {
269 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1];
270 h__[k + j * h_dim1] -= p * x;
271 h__[k + 1 + j * h_dim1] -= p * y;
272 /* L200: */
273 }
274
275 /* Computing MIN */
276 i__2 = en, i__3 = k + 3;
277 j = min(i__2,i__3);
278 /* .......... COLUMN MODIFICATION .......... */
279 i__2 = j;
280 for (i__ = 1; i__ <= i__2; ++i__) {
281 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1];
282 h__[i__ + k * h_dim1] -= p;
283 h__[i__ + (k + 1) * h_dim1] -= p * q;
284 /* L210: */
285 }
286 goto L255;
287 L225:
288 /* .......... ROW MODIFICATION .......... */
289 i__2 = *n;
290 for (j = k; j <= i__2; ++j) {
291 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[
292 k + 2 + j * h_dim1];
293 h__[k + j * h_dim1] -= p * x;
294 h__[k + 1 + j * h_dim1] -= p * y;
295 h__[k + 2 + j * h_dim1] -= p * zz;
296 /* L230: */
297 }
298
299 /* Computing MIN */
300 i__2 = en, i__3 = k + 3;
301 j = min(i__2,i__3);
302 /* .......... COLUMN MODIFICATION .......... */
303 i__2 = j;
304 for (i__ = 1; i__ <= i__2; ++i__) {
305 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] +
306 zz * h__[i__ + (k + 2) * h_dim1];
307 h__[i__ + k * h_dim1] -= p;
308 h__[i__ + (k + 1) * h_dim1] -= p * q;
309 h__[i__ + (k + 2) * h_dim1] -= p * r__;
310 /* L240: */
311 }
312 L255:
313
314 L260:
315 ;
316 }
317
318 goto L70;
319 /* .......... ONE ROOT FOUND .......... */
320 L270:
321 wr[en] = x + t;
322 wi[en] = 0.;
323 en = na;
324 goto L60;
325 /* .......... TWO ROOTS FOUND .......... */
326 L280:
327 p = (y - x) / 2.;
328 q = p * p + w;
329 zz = sqrt((abs(q)));
330 x += t;
331 if (q < 0.) {
332 goto L320;
333 }
334 /* .......... REAL PAIR .......... */
335 zz = p + d_sign(&zz, &p);
336 wr[na] = x + zz;
337 wr[en] = wr[na];
338 if (zz != 0.) {
339 wr[en] = x - w / zz;
340 }
341 wi[na] = 0.;
342 wi[en] = 0.;
343 goto L330;
344 /* .......... COMPLEX PAIR .......... */
345 L320:
346 wr[na] = x + p;
347 wr[en] = x + p;
348 wi[na] = zz;
349 wi[en] = -zz;
350 L330:
351 en = enm2;
352 goto L60;
353 /* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
354 /* CONVERGED AFTER 30*N ITERATIONS .......... */
355 L1000:
356 *ierr = en;
357 L1001:
358 return 0;
359 } /* hqr_ */
360
361