lib/eigen.c

   1 /* eigen.f -- translated by f2c (version of 19 December 1990  16:50:21).
   2    and modified. */
   3
   4 #include "xmath.h"
   5
   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))
  12
  13 static /* Subroutine */ int tred2(), tql2();
  14 static doublereal pythag(), d_sign();
  15
  16 /* Table of constant values */
  17
  18 static doublereal c_b39 = 1.;
  19
  20 /* Subroutine */ int eigen(nm, n, a, w, z, fv1, ierr)
  21 integer *nm, *n;
  22 doublereal *a, *w, *z, *fv1;
  23 integer *ierr;
  24 {
  25     /* System generated locals */
  26     integer a_dim1, a_offset, z_dim1, z_offset;
  27
  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. */
  32
  33 /*     ON INPUT */
  34
  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. */
  38
  39 /*        N  IS THE ORDER OF THE MATRIX  A. */
  40
  41 /*        A  CONTAINS THE REAL SYMMETRIC MATRIX. */
  42
  43 /*     ON OUTPUT */
  44
  45 /*        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER. */
  46
  47 /*        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. */
  48
  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. */
  52
  53 /*        FV1 IS A TEMPORARY STORAGE ARRAY. */
  54
  55 /*     ------------------------------------------------------------------
  56 */
  57
  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;
  67
  68     /* Function Body */
  69     if (*n <= *nm) {
  70         goto L10;
  71     }
  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 */
  81
  82 static /* Subroutine */ int tred2(nm, n, a, d, e, z)
  83 integer *nm, *n;
  84 doublereal *a, *d, *e, *z;
  85 {
  86     /* System generated locals */
  87     integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
  88     doublereal d__1;
  89
  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;
  95
  96
  97
  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). */
 101
 102 /*     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A */
 103 /*     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING */
 104 /*     ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
 105
 106 /*     ON INPUT */
 107
 108 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
 109 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
 110 /*          DIMENSION STATEMENT. */
 111
 112 /*        N IS THE ORDER OF THE MATRIX. */
 113
 114 /*        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE */
 115 /*          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */
 116
 117 /*     ON OUTPUT */
 118
 119 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
 120
 121 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
 122 /*          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO. */
 123
 124 /*        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX */
 125 /*          PRODUCED IN THE REDUCTION. */
 126
 127 /*        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED. */
 128
 129 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
 130 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
 131 */
 132
 133 /*     THIS VERSION DATED AUGUST 1983. */
 134
 135 /*     ------------------------------------------------------------------
 136 */
 137
 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;
 147
 148     /* Function Body */
 149     i__1 = *n;
 150     for (i = 1; i <= i__1; ++i) {
 151
 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];
 156         }
 157
 158         d[i] = a[*n + i * a_dim1];
 159 /* L100: */
 160     }
 161
 162     if (*n == 1) {
 163         goto L510;
 164     }
 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;
 174         }
 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));
 180         }
 181
 182         if (scale != 0.) {
 183             goto L140;
 184         }
 185 L130:
 186         e[i] = d[l];
 187
 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: */
 194         }
 195
 196         goto L290;
 197
 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: */
 204         }
 205
 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.;
 217         }
 218
 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;
 227             }
 228
 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: */
 234             }
 235
 236 L220:
 237             e[j] = g;
 238 /* L240: */
 239         }
 240 /*     .......... FORM P .......... */
 241         f = 0.;
 242
 243         i__2 = l;
 244         for (j = 1; j <= i__2; ++j) {
 245             e[j] /= h;
 246             f += e[j] * d[j];
 247 /* L245: */
 248         }
 249
 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];
 256         }
 257 /*     .......... FORM REDUCED A .......... */
 258         i__2 = l;
 259         for (j = 1; j <= i__2; ++j) {
 260             f = d[j];
 261             g = e[j];
 262
 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];
 267             }
 268
 269             d[j] = z[l + j * z_dim1];
 270             z[i + j * z_dim1] = 0.;
 271 /* L280: */
 272         }
 273
 274 L290:
 275         d[i] = h;
 276 /* L300: */
 277     }
 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;
 287         }
 288
 289         i__2 = l;
 290         for (k = 1; k <= i__2; ++k) {
 291 /* L330: */
 292             d[k] = z[k + i * z_dim1] / h;
 293         }
 294
 295         i__2 = l;
 296         for (j = 1; j <= i__2; ++j) {
 297             g = 0.;
 298
 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];
 303             }
 304
 305             i__3 = l;
 306             for (k = 1; k <= i__3; ++k) {
 307                 z[k + j * z_dim1] -= g * d[k];
 308 /* L360: */
 309             }
 310         }
 311
 312 L380:
 313         i__3 = l;
 314         for (k = 1; k <= i__3; ++k) {
 315 /* L400: */
 316             z[k + i * z_dim1] = 0.;
 317         }
 318
 319 /* L500: */
 320     }
 321
 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: */
 328     }
 329
 330     z[*n + *n * z_dim1] = 1.;
 331     e[1] = 0.;
 332     return 0;
 333 } /* tred2 */
 334
 335 static /* Subroutine */ int tql2(nm, n, d, e, z, ierr)
 336 integer *nm, *n;
 337 doublereal *d, *e, *z;
 338 integer *ierr;
 339 {
 340     /* System generated locals */
 341     integer z_dim1, z_offset, i__1, i__2, i__3;
 342     doublereal d__1, d__2;
 343
 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;
 354
 355
 356
 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). */
 361
 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. */
 367
 368 /*     ON INPUT */
 369
 370 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
 371 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
 372 /*          DIMENSION STATEMENT. */
 373
 374 /*        N IS THE ORDER OF THE MATRIX. */
 375
 376 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
 377
 378 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
 379 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
 380
 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. */
 385
 386 /*      ON OUTPUT */
 387
 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. */
 391
 392 /*        E HAS BEEN DESTROYED. */
 393
 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. */
 398
 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. */
 403
 404 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
 405
 406 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
 407 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
 408 */
 409
 410 /*     THIS VERSION DATED AUGUST 1983. */
 411
 412 /*     ------------------------------------------------------------------
 413 */
 414
 415     /* Parameter adjustments */
 416     z_dim1 = *nm;
 417     z_offset = z_dim1 + 1;
 418     z -= z_offset;
 419     --e;
 420     --d;
 421
 422     /* Function Body */
 423     *ierr = 0;
 424     if (*n == 1) {
 425         goto L1001;
 426     }
 427
 428     i__1 = *n;
 429     for (i = 2; i <= i__1; ++i) {
 430 /* L100: */
 431         e[i - 1] = e[i];
 432     }
 433
 434     f = 0.;
 435     tst1 = 0.;
 436     e[*n] = 0.;
 437
 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;
 444         }
 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;
 451             }
 452 /*     .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT */
 453 /*                THROUGH THE BOTTOM OF THE LOOP .......... */
 454 /* L110: */
 455         }
 456
 457 L120:
 458         if (m == l) {
 459             goto L220;
 460         }
 461 L130:
 462         if (j == 30) {
 463             goto L1000;
 464         }
 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;
 478         }
 479
 480         i__2 = *n;
 481         for (i = l2; i <= i__2; ++i) {
 482 /* L140: */
 483             d[i] -= h;
 484         }
 485
 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: */
 517             }
 518
 519 /* L200: */
 520         }
 521
 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;
 528         }
 529 L220:
 530         d[l] += f;
 531 /* L240: */
 532     }
 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];
 539
 540         i__2 = *n;
 541         for (j = ii; j <= i__2; ++j) {
 542             if (d[j] >= p) {
 543                 goto L260;
 544             }
 545             k = j;
 546             p = d[j];
 547 L260:
 548             ;
 549         }
 550
 551         if (k == i) {
 552             goto L300;
 553         }
 554         d[k] = d[i];
 555         d[i] = p;
 556
 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: */
 563         }
 564
 565 L300:
 566         ;
 567     }
 568
 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 */
 577
 578 static doublereal pythag(a, b)
 579 doublereal *a, *b;
 580 {
 581     /* System generated locals */
 582     doublereal ret_val, d__1, d__2, d__3;
 583
 584     /* Local variables */
 585     static doublereal p, r, s, t, u;
 586
 587
 588 /*     FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */
 589
 590
 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;
 596     }
 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;
 606     }
 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 */
 618
 619 static double d_sign(a,b)
 620      doublereal *a, *b;
 621 {
 622   double x;
 623   x = (*a >= 0 ? *a : - *a);
 624   return( *b >= 0 ? x : -x);
 625 }