removed Lisp-specific help
[tsl.git] / lib / eigen.c
blobe45d4f12471c2eb60540370bbe57730753d3a1a5
1 /* eigen.f -- translated by f2c (version of 19 December 1990 16:50:21).
2 and modified. */
4 #include "xmath.h"
6 #define integer int
7 #define real float
8 #define doublereal double
9 #define Min(x,y) ((x) > (y) ? (y) : (x))
10 #define Max(x,y) ((x) > (y) ? (x) : (y))
11 #define Abs(x) ((x) >= 0 ? (x) : -(x))
13 static /* Subroutine */ int tred2(), tql2();
14 static doublereal pythag(), d_sign();
16 /* Table of constant values */
18 static doublereal c_b39 = 1.;
20 /* Subroutine */ int eigen(nm, n, a, w, z, fv1, ierr)
21 integer *nm, *n;
22 doublereal *a, *w, *z, *fv1;
23 integer *ierr;
25 /* System generated locals */
26 integer a_dim1, a_offset, z_dim1, z_offset;
28 /* THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF */
29 /* SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) */
30 /* TO FIND THE EIGENVALUES AND EIGENVECTORS */
31 /* OF A REAL SYMMETRIC MATRIX. */
33 /* ON INPUT */
35 /* NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL */
36 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
37 /* DIMENSION STATEMENT. */
39 /* N IS THE ORDER OF THE MATRIX A. */
41 /* A CONTAINS THE REAL SYMMETRIC MATRIX. */
43 /* ON OUTPUT */
45 /* W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. */
47 /* Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. */
49 /* IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR */
50 /* COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT */
51 /* AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. */
53 /* FV1 IS A TEMPORARY STORAGE ARRAY. */
55 /* ------------------------------------------------------------------
58 /* Parameter adjustments */
59 --fv1;
60 z_dim1 = *nm;
61 z_offset = z_dim1 + 1;
62 z -= z_offset;
63 --w;
64 a_dim1 = *nm;
65 a_offset = a_dim1 + 1;
66 a -= a_offset;
68 /* Function Body */
69 if (*n <= *nm) {
70 goto L10;
72 *ierr = *n * 10;
73 goto L50;
74 /* .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... */
75 L10:
76 tred2(nm, n, &a[a_offset], &w[1], &fv1[1], &z[z_offset]);
77 tql2(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
78 L50:
79 return 0;
80 } /* eigen */
82 static /* Subroutine */ int tred2(nm, n, a, d, e, z)
83 integer *nm, *n;
84 doublereal *a, *d, *e, *z;
86 /* System generated locals */
87 integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
88 doublereal d__1;
90 /* Local variables */
91 static doublereal f, g, h;
92 static integer i, j, k, l;
93 static doublereal scale, hh;
94 static integer ii, jp1;
98 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2, */
99 /* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. */
100 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */
102 /* THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A */
103 /* SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING */
104 /* ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
106 /* ON INPUT */
108 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
109 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
110 /* DIMENSION STATEMENT. */
112 /* N IS THE ORDER OF THE MATRIX. */
114 /* A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE */
115 /* LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */
117 /* ON OUTPUT */
119 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
121 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
122 /* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */
124 /* Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX */
125 /* PRODUCED IN THE REDUCTION. */
127 /* A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED. */
129 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
130 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
133 /* THIS VERSION DATED AUGUST 1983. */
135 /* ------------------------------------------------------------------
138 /* Parameter adjustments */
139 z_dim1 = *nm;
140 z_offset = z_dim1 + 1;
141 z -= z_offset;
142 --e;
143 --d;
144 a_dim1 = *nm;
145 a_offset = a_dim1 + 1;
146 a -= a_offset;
148 /* Function Body */
149 i__1 = *n;
150 for (i = 1; i <= i__1; ++i) {
152 i__2 = *n;
153 for (j = i; j <= i__2; ++j) {
154 /* L80: */
155 z[j + i * z_dim1] = a[j + i * a_dim1];
158 d[i] = a[*n + i * a_dim1];
159 /* L100: */
162 if (*n == 1) {
163 goto L510;
165 /* .......... FOR I=N STEP -1 UNTIL 2 DO -- .......... */
166 i__1 = *n;
167 for (ii = 2; ii <= i__1; ++ii) {
168 i = *n + 2 - ii;
169 l = i - 1;
170 h = 0.;
171 scale = 0.;
172 if (l < 2) {
173 goto L130;
175 /* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */
176 i__2 = l;
177 for (k = 1; k <= i__2; ++k) {
178 /* L120: */
179 scale += (d__1 = d[k], Abs(d__1));
182 if (scale != 0.) {
183 goto L140;
185 L130:
186 e[i] = d[l];
188 i__2 = l;
189 for (j = 1; j <= i__2; ++j) {
190 d[j] = z[l + j * z_dim1];
191 z[i + j * z_dim1] = 0.;
192 z[j + i * z_dim1] = 0.;
193 /* L135: */
196 goto L290;
198 L140:
199 i__2 = l;
200 for (k = 1; k <= i__2; ++k) {
201 d[k] /= scale;
202 h += d[k] * d[k];
203 /* L150: */
206 f = d[l];
207 d__1 = sqrt(h);
208 g = -d_sign(&d__1, &f);
209 e[i] = scale * g;
210 h -= f * g;
211 d[l] = f - g;
212 /* .......... FORM A*U .......... */
213 i__2 = l;
214 for (j = 1; j <= i__2; ++j) {
215 /* L170: */
216 e[j] = 0.;
219 i__2 = l;
220 for (j = 1; j <= i__2; ++j) {
221 f = d[j];
222 z[j + i * z_dim1] = f;
223 g = e[j] + z[j + j * z_dim1] * f;
224 jp1 = j + 1;
225 if (l < jp1) {
226 goto L220;
229 i__3 = l;
230 for (k = jp1; k <= i__3; ++k) {
231 g += z[k + j * z_dim1] * d[k];
232 e[k] += z[k + j * z_dim1] * f;
233 /* L200: */
236 L220:
237 e[j] = g;
238 /* L240: */
240 /* .......... FORM P .......... */
241 f = 0.;
243 i__2 = l;
244 for (j = 1; j <= i__2; ++j) {
245 e[j] /= h;
246 f += e[j] * d[j];
247 /* L245: */
250 hh = f / (h + h);
251 /* .......... FORM Q .......... */
252 i__2 = l;
253 for (j = 1; j <= i__2; ++j) {
254 /* L250: */
255 e[j] -= hh * d[j];
257 /* .......... FORM REDUCED A .......... */
258 i__2 = l;
259 for (j = 1; j <= i__2; ++j) {
260 f = d[j];
261 g = e[j];
263 i__3 = l;
264 for (k = j; k <= i__3; ++k) {
265 /* L260: */
266 z[k + j * z_dim1] = z[k + j * z_dim1] - f * e[k] - g * d[k];
269 d[j] = z[l + j * z_dim1];
270 z[i + j * z_dim1] = 0.;
271 /* L280: */
274 L290:
275 d[i] = h;
276 /* L300: */
278 /* .......... ACCUMULATION OF TRANSFORMATION MATRICES .......... */
279 i__1 = *n;
280 for (i = 2; i <= i__1; ++i) {
281 l = i - 1;
282 z[*n + l * z_dim1] = z[l + l * z_dim1];
283 z[l + l * z_dim1] = 1.;
284 h = d[i];
285 if (h == 0.) {
286 goto L380;
289 i__2 = l;
290 for (k = 1; k <= i__2; ++k) {
291 /* L330: */
292 d[k] = z[k + i * z_dim1] / h;
295 i__2 = l;
296 for (j = 1; j <= i__2; ++j) {
297 g = 0.;
299 i__3 = l;
300 for (k = 1; k <= i__3; ++k) {
301 /* L340: */
302 g += z[k + i * z_dim1] * z[k + j * z_dim1];
305 i__3 = l;
306 for (k = 1; k <= i__3; ++k) {
307 z[k + j * z_dim1] -= g * d[k];
308 /* L360: */
312 L380:
313 i__3 = l;
314 for (k = 1; k <= i__3; ++k) {
315 /* L400: */
316 z[k + i * z_dim1] = 0.;
319 /* L500: */
322 L510:
323 i__1 = *n;
324 for (i = 1; i <= i__1; ++i) {
325 d[i] = z[*n + i * z_dim1];
326 z[*n + i * z_dim1] = 0.;
327 /* L520: */
330 z[*n + *n * z_dim1] = 1.;
331 e[1] = 0.;
332 return 0;
333 } /* tred2 */
335 static /* Subroutine */ int tql2(nm, n, d, e, z, ierr)
336 integer *nm, *n;
337 doublereal *d, *e, *z;
338 integer *ierr;
340 /* System generated locals */
341 integer z_dim1, z_offset, i__1, i__2, i__3;
342 doublereal d__1, d__2;
344 /* Local variables */
345 static doublereal c, f, g, h;
346 static integer i, j, k, l, m;
347 static doublereal p, r, s, c2, c3;
348 static integer l1, l2;
349 static doublereal s2;
350 static integer ii;
351 static doublereal dl1, el1;
352 static integer mml;
353 static doublereal tst1, tst2;
357 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, */
358 /* NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND */
359 /* WILKINSON. */
360 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). */
362 /* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
363 /* OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. */
364 /* THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO */
365 /* BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS */
366 /* FULL MATRIX TO TRIDIAGONAL FORM. */
368 /* ON INPUT */
370 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
371 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
372 /* DIMENSION STATEMENT. */
374 /* N IS THE ORDER OF THE MATRIX. */
376 /* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
378 /* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
379 /* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */
381 /* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */
382 /* REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS */
383 /* OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN */
384 /* THE IDENTITY MATRIX. */
386 /* ON OUTPUT */
388 /* D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN */
389 /* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT */
390 /* UNORDERED FOR INDICES 1,2,...,IERR-1. */
392 /* E HAS BEEN DESTROYED. */
394 /* Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC */
395 /* TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, */
396 /* Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED */
397 /* EIGENVALUES. */
399 /* IERR IS SET TO */
400 /* ZERO FOR NORMAL RETURN, */
401 /* J IF THE J-TH EIGENVALUE HAS NOT BEEN */
402 /* DETERMINED AFTER 30 ITERATIONS. */
404 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */
406 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
407 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
410 /* THIS VERSION DATED AUGUST 1983. */
412 /* ------------------------------------------------------------------
415 /* Parameter adjustments */
416 z_dim1 = *nm;
417 z_offset = z_dim1 + 1;
418 z -= z_offset;
419 --e;
420 --d;
422 /* Function Body */
423 *ierr = 0;
424 if (*n == 1) {
425 goto L1001;
428 i__1 = *n;
429 for (i = 2; i <= i__1; ++i) {
430 /* L100: */
431 e[i - 1] = e[i];
434 f = 0.;
435 tst1 = 0.;
436 e[*n] = 0.;
438 i__1 = *n;
439 for (l = 1; l <= i__1; ++l) {
440 j = 0;
441 h = (d__1 = d[l], Abs(d__1)) + (d__2 = e[l], Abs(d__2));
442 if (tst1 < h) {
443 tst1 = h;
445 /* .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
446 i__2 = *n;
447 for (m = l; m <= i__2; ++m) {
448 tst2 = tst1 + (d__1 = e[m], Abs(d__1));
449 if (tst2 == tst1) {
450 goto L120;
452 /* .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT */
453 /* THROUGH THE BOTTOM OF THE LOOP .......... */
454 /* L110: */
457 L120:
458 if (m == l) {
459 goto L220;
461 L130:
462 if (j == 30) {
463 goto L1000;
465 ++j;
466 /* .......... FORM SHIFT .......... */
467 l1 = l + 1;
468 l2 = l1 + 1;
469 g = d[l];
470 p = (d[l1] - g) / (e[l] * 2.);
471 r = pythag(&p, &c_b39);
472 d[l] = e[l] / (p + d_sign(&r, &p));
473 d[l1] = e[l] * (p + d_sign(&r, &p));
474 dl1 = d[l1];
475 h = g - d[l];
476 if (l2 > *n) {
477 goto L145;
480 i__2 = *n;
481 for (i = l2; i <= i__2; ++i) {
482 /* L140: */
483 d[i] -= h;
486 L145:
487 f += h;
488 /* .......... QL TRANSFORMATION .......... */
489 p = d[m];
490 c = 1.;
491 c2 = c;
492 el1 = e[l1];
493 s = 0.;
494 mml = m - l;
495 /* .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
496 i__2 = mml;
497 for (ii = 1; ii <= i__2; ++ii) {
498 c3 = c2;
499 c2 = c;
500 s2 = s;
501 i = m - ii;
502 g = c * e[i];
503 h = c * p;
504 r = pythag(&p, &e[i]);
505 e[i + 1] = s * r;
506 s = e[i] / r;
507 c = p / r;
508 p = c * d[i] - s * g;
509 d[i + 1] = h + s * (c * g + s * d[i]);
510 /* .......... FORM VECTOR .......... */
511 i__3 = *n;
512 for (k = 1; k <= i__3; ++k) {
513 h = z[k + (i + 1) * z_dim1];
514 z[k + (i + 1) * z_dim1] = s * z[k + i * z_dim1] + c * h;
515 z[k + i * z_dim1] = c * z[k + i * z_dim1] - s * h;
516 /* L180: */
519 /* L200: */
522 p = -s * s2 * c3 * el1 * e[l] / dl1;
523 e[l] = s * p;
524 d[l] = c * p;
525 tst2 = tst1 + (d__1 = e[l], Abs(d__1));
526 if (tst2 > tst1) {
527 goto L130;
529 L220:
530 d[l] += f;
531 /* L240: */
533 /* .......... ORDER EIGENVALUES AND EIGENVECTORS .......... */
534 i__1 = *n;
535 for (ii = 2; ii <= i__1; ++ii) {
536 i = ii - 1;
537 k = i;
538 p = d[i];
540 i__2 = *n;
541 for (j = ii; j <= i__2; ++j) {
542 if (d[j] >= p) {
543 goto L260;
545 k = j;
546 p = d[j];
547 L260:
551 if (k == i) {
552 goto L300;
554 d[k] = d[i];
555 d[i] = p;
557 i__2 = *n;
558 for (j = 1; j <= i__2; ++j) {
559 p = z[j + i * z_dim1];
560 z[j + i * z_dim1] = z[j + k * z_dim1];
561 z[j + k * z_dim1] = p;
562 /* L280: */
565 L300:
569 goto L1001;
570 /* .......... SET ERROR -- NO CONVERGENCE TO AN */
571 /* EIGENVALUE AFTER 30 ITERATIONS .......... */
572 L1000:
573 *ierr = l;
574 L1001:
575 return 0;
576 } /* tql2 */
578 static doublereal pythag(a, b)
579 doublereal *a, *b;
581 /* System generated locals */
582 doublereal ret_val, d__1, d__2, d__3;
584 /* Local variables */
585 static doublereal p, r, s, t, u;
588 /* FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */
591 /* Computing MAX */
592 d__1 = Abs(*a), d__2 = Abs(*b);
593 p = Max(d__1,d__2);
594 if (p == 0.) {
595 goto L20;
597 /* Computing MIN */
598 d__2 = Abs(*a), d__3 = Abs(*b);
599 /* Computing 2nd power */
600 d__1 = Min(d__2,d__3) / p;
601 r = d__1 * d__1;
602 L10:
603 t = r + 4.;
604 if (t == 4.) {
605 goto L20;
607 s = r / t;
608 u = s * 2. + 1.;
609 p = u * p;
610 /* Computing 2nd power */
611 d__1 = s / u;
612 r = d__1 * d__1 * r;
613 goto L10;
614 L20:
615 ret_val = p;
616 return ret_val;
617 } /* pythag */
619 static double d_sign(a,b)
620 doublereal *a, *b;
622 double x;
623 x = (*a >= 0 ? *a : - *a);
624 return( *b >= 0 ? x : -x);