website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_c/routines/eigen.c

   1 #include <math.h>
   2 #include <stdlib.h>
   3 #include <stdio.h>
   4 #include "tisean_cec.h"
   5
   6 typedef double doublereal;
   7 typedef int integer;
   8
   9 #define abs(x) (((x)>=0.0)?(x):-(x))
  10 #define min(x,y) (((x)<=(y))?(x):(y))
  11 #define max(x,y) (((x)>=(y))?(x):(y))
  12
  13 static doublereal c_b10 = 1.;
  14
  15 extern void check_alloc(void*);
  16
  17 double d_sign(double *a,double *b)
  18 {
  19   double x;
  20   x = (*a >= 0 ? *a : - *a);
  21   return ( *b >= 0 ? x : -x);
  22 }
  23
  24 doublereal pythag(doublereal *a, doublereal *b)
  25 {
  26     doublereal ret_val, d__1, d__2, d__3;
  27     static doublereal p, r__, s, t, u;
  28
  29     d__1 = abs(*a), d__2 = abs(*b);
  30     p = max(d__1,d__2);
  31     if (p == 0.) {
  32         goto L20;
  33     }
  34     d__2 = abs(*a), d__3 = abs(*b);
  35     d__1 = min(d__2,d__3) / p;
  36     r__ = d__1 * d__1;
  37 L10:
  38     t = r__ + 4.;
  39     if (t == 4.) {
  40         goto L20;
  41     }
  42     s = r__ / t;
  43     u = s * 2. + 1.;
  44     p = u * p;
  45     d__1 = s / u;
  46     r__ = d__1 * d__1 * r__;
  47     goto L10;
  48 L20:
  49     ret_val = p;
  50     return ret_val;
  51 }
  52
  53
  54 int tred2(integer *nm, integer *n, doublereal *a,
  55         doublereal *d__, doublereal *e, doublereal *z__)
  56 {
  57     integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
  58     doublereal d__1;
  59
  60     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
  61
  62     static doublereal f, g, h__;
  63     static integer i__, j, k, l;
  64     static doublereal hh;
  65     static integer ii, jp1;
  66     static doublereal scale;
  67
  68
  69
  70 /*     this subroutine is a translation of the algol procedure tred2, */
  71 /*     num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson. */
  72 /*     handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */
  73
  74 /*     this subroutine reduces a real symmetric matrix to a */
  75 /*     symmetric tridiagonal matrix using and accumulating */
  76 /*     orthogonal similarity transformations. */
  77
  78 /*     on input */
  79
  80 /*        nm must be set to the row dimension of two-dimensional */
  81 /*          array parameters as declared in the calling program */
  82 /*          dimension statement. */
  83
  84 /*        n is the order of the matrix. */
  85
  86 /*        a contains the real symmetric input matrix.  only the */
  87 /*          lower triangle of the matrix need be supplied. */
  88
  89 /*     on output */
  90
  91 /*        d contains the diagonal elements of the tridiagonal matrix. */
  92
  93 /*        e contains the subdiagonal elements of the tridiagonal */
  94 /*          matrix in its last n-1 positions.  e(1) is set to zero. */
  95
  96 /*        z contains the orthogonal transformation matrix */
  97 /*          produced in the reduction. */
  98
  99 /*        a and z may coincide.  if distinct, a is unaltered. */
 100
 101 /*     questions and comments should be directed to burton s. garbow, */
 102 /*     mathematics and computer science div, argonne national laboratory */
 103
 104 /*     this version dated august 1983. */
 105
 106 /*     ------------------------------------------------------------------ */
 107
 108     z_dim1 = *nm;
 109     z_offset = 1 + z_dim1 * 1;
 110     z__ -= z_offset;
 111     --e;
 112     --d__;
 113     a_dim1 = *nm;
 114     a_offset = 1 + a_dim1 * 1;
 115     a -= a_offset;
 116
 117     i__1 = *n;
 118     for (i__ = 1; i__ <= i__1; ++i__) {
 119
 120         i__2 = *n;
 121         for (j = i__; j <= i__2; ++j) {
 122             z__[j + i__ * z_dim1] = a[j + i__ * a_dim1];
 123         }
 124
 125         d__[i__] = a[*n + i__ * a_dim1];
 126     }
 127
 128     if (*n == 1) {
 129         goto L510;
 130     }
 131     i__1 = *n;
 132     for (ii = 2; ii <= i__1; ++ii) {
 133         i__ = *n + 2 - ii;
 134         l = i__ - 1;
 135         h__ = 0.;
 136         scale = 0.;
 137         if (l < 2) {
 138             goto L130;
 139         }
 140         i__2 = l;
 141         for (k = 1; k <= i__2; ++k) {
 142             scale += (d__1 = d__[k], abs(d__1));
 143         }
 144
 145         if (scale != 0.) {
 146             goto L140;
 147         }
 148 L130:
 149         e[i__] = d__[l];
 150
 151         i__2 = l;
 152         for (j = 1; j <= i__2; ++j) {
 153             d__[j] = z__[l + j * z_dim1];
 154             z__[i__ + j * z_dim1] = 0.;
 155             z__[j + i__ * z_dim1] = 0.;
 156         }
 157
 158         goto L290;
 159
 160 L140:
 161         i__2 = l;
 162         for (k = 1; k <= i__2; ++k) {
 163             d__[k] /= scale;
 164             h__ += d__[k] * d__[k];
 165         }
 166
 167         f = d__[l];
 168         d__1 = sqrt(h__);
 169         g = -d_sign(&d__1, &f);
 170         e[i__] = scale * g;
 171         h__ -= f * g;
 172         d__[l] = f - g;
 173         i__2 = l;
 174         for (j = 1; j <= i__2; ++j) {
 175             e[j] = 0.;
 176         }
 177
 178         i__2 = l;
 179         for (j = 1; j <= i__2; ++j) {
 180             f = d__[j];
 181             z__[j + i__ * z_dim1] = f;
 182             g = e[j] + z__[j + j * z_dim1] * f;
 183             jp1 = j + 1;
 184             if (l < jp1) {
 185                 goto L220;
 186             }
 187
 188             i__3 = l;
 189             for (k = jp1; k <= i__3; ++k) {
 190                 g += z__[k + j * z_dim1] * d__[k];
 191                 e[k] += z__[k + j * z_dim1] * f;
 192             }
 193
 194 L220:
 195             e[j] = g;
 196         }
 197         f = 0.;
 198
 199         i__2 = l;
 200         for (j = 1; j <= i__2; ++j) {
 201             e[j] /= h__;
 202             f += e[j] * d__[j];
 203         }
 204
 205         hh = f / (h__ + h__);
 206         i__2 = l;
 207         for (j = 1; j <= i__2; ++j) {
 208             e[j] -= hh * d__[j];
 209         }
 210         i__2 = l;
 211         for (j = 1; j <= i__2; ++j) {
 212             f = d__[j];
 213             g = e[j];
 214
 215             i__3 = l;
 216             for (k = j; k <= i__3; ++k) {
 217                 z__[k + j * z_dim1] = z__[k + j * z_dim1] - f * e[k] - g *
 218                         d__[k];
 219             }
 220
 221             d__[j] = z__[l + j * z_dim1];
 222             z__[i__ + j * z_dim1] = 0.;
 223         }
 224
 225 L290:
 226         d__[i__] = h__;
 227     }
 228     i__1 = *n;
 229     for (i__ = 2; i__ <= i__1; ++i__) {
 230         l = i__ - 1;
 231         z__[*n + l * z_dim1] = z__[l + l * z_dim1];
 232         z__[l + l * z_dim1] = 1.;
 233         h__ = d__[i__];
 234         if (h__ == 0.) {
 235             goto L380;
 236         }
 237
 238         i__2 = l;
 239         for (k = 1; k <= i__2; ++k) {
 240             d__[k] = z__[k + i__ * z_dim1] / h__;
 241         }
 242
 243         i__2 = l;
 244         for (j = 1; j <= i__2; ++j) {
 245             g = 0.;
 246
 247             i__3 = l;
 248             for (k = 1; k <= i__3; ++k) {
 249                 g += z__[k + i__ * z_dim1] * z__[k + j * z_dim1];
 250             }
 251
 252             i__3 = l;
 253             for (k = 1; k <= i__3; ++k) {
 254                 z__[k + j * z_dim1] -= g * d__[k];
 255             }
 256         }
 257
 258 L380:
 259         i__3 = l;
 260         for (k = 1; k <= i__3; ++k) {
 261             z__[k + i__ * z_dim1] = 0.;
 262         }
 263
 264     }
 265
 266 L510:
 267     i__1 = *n;
 268     for (i__ = 1; i__ <= i__1; ++i__) {
 269         d__[i__] = z__[*n + i__ * z_dim1];
 270         z__[*n + i__ * z_dim1] = 0.;
 271     }
 272
 273     z__[*n + *n * z_dim1] = 1.;
 274     e[1] = 0.;
 275     return 0;
 276 }
 277
 278 int tql2(integer *nm, integer *n, doublereal *d__,
 279         doublereal *e, doublereal *z__, integer *ierr)
 280 {
 281     integer z_dim1, z_offset, i__1, i__2, i__3;
 282     doublereal d__1, d__2;
 283
 284     double d_sign(doublereal *, doublereal *);
 285
 286     static doublereal c__, f, g, h__;
 287     static integer i__, j, k, l, m;
 288     static doublereal p, r__, s, c2, c3;
 289     static integer l1, l2;
 290     static doublereal s2;
 291     static integer ii;
 292     static doublereal dl1, el1;
 293     static integer mml;
 294     static doublereal tst1, tst2;
 295     extern doublereal pythag_(doublereal *, doublereal *);
 296
 297
 298
 299 /*     this subroutine is a translation of the algol procedure tql2, */
 300 /*     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and */
 301 /*     wilkinson. */
 302 /*     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). */
 303
 304 /*     this subroutine finds the eigenvalues and eigenvectors */
 305 /*     of a symmetric tridiagonal matrix by the ql method. */
 306 /*     the eigenvectors of a full symmetric matrix can also */
 307 /*     be found if  tred2  has been used to reduce this */
 308 /*     full matrix to tridiagonal form. */
 309
 310 /*     on input */
 311
 312 /*        nm must be set to the row dimension of two-dimensional */
 313 /*          array parameters as declared in the calling program */
 314 /*          dimension statement. */
 315
 316 /*        n is the order of the matrix. */
 317
 318 /*        d contains the diagonal elements of the input matrix. */
 319
 320 /*        e contains the subdiagonal elements of the input matrix */
 321 /*          in its last n-1 positions.  e(1) is arbitrary. */
 322
 323 /*        z contains the transformation matrix produced in the */
 324 /*          reduction by  tred2, if performed.  if the eigenvectors */
 325 /*          of the tridiagonal matrix are desired, z must contain */
 326 /*          the identity matrix. */
 327
 328 /*      on output */
 329
 330 /*        d contains the eigenvalues in ascending order.  if an */
 331 /*          error exit is made, the eigenvalues are correct but */
 332 /*          unordered for indices 1,2,...,ierr-1. */
 333
 334 /*        e has been destroyed. */
 335
 336 /*        z contains orthonormal eigenvectors of the symmetric */
 337 /*          tridiagonal (or full) matrix.  if an error exit is made, */
 338 /*          z contains the eigenvectors associated with the stored */
 339 /*          eigenvalues. */
 340
 341 /*        ierr is set to */
 342 /*          zero       for normal return, */
 343 /*          j          if the j-th eigenvalue has not been */
 344 /*                     determined after 30 iterations. */
 345
 346 /*     calls pythag for  dsqrt(a*a + b*b) . */
 347
 348 /*     questions and comments should be directed to burton s. garbow, */
 349 /*     mathematics and computer science div, argonne national laboratory */
 350
 351 /*     this version dated august 1983. */
 352
 353 /*     ------------------------------------------------------------------ */
 354
 355     z_dim1 = *nm;
 356     z_offset = 1 + z_dim1 * 1;
 357     z__ -= z_offset;
 358     --e;
 359     --d__;
 360
 361     *ierr = 0;
 362     if (*n == 1) {
 363         goto L1001;
 364     }
 365
 366     i__1 = *n;
 367     for (i__ = 2; i__ <= i__1; ++i__) {
 368         e[i__ - 1] = e[i__];
 369     }
 370
 371     f = 0.;
 372     tst1 = 0.;
 373     e[*n] = 0.;
 374
 375     i__1 = *n;
 376     for (l = 1; l <= i__1; ++l) {
 377         j = 0;
 378         h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2));
 379         if (tst1 < h__) {
 380             tst1 = h__;
 381         }
 382         i__2 = *n;
 383         for (m = l; m <= i__2; ++m) {
 384             tst2 = tst1 + (d__1 = e[m], abs(d__1));
 385             if (tst2 == tst1) {
 386                 goto L120;
 387             }
 388         }
 389
 390 L120:
 391         if (m == l) {
 392             goto L220;
 393         }
 394 L130:
 395         if (j == 30) {
 396             goto L1000;
 397         }
 398         ++j;
 399         l1 = l + 1;
 400         l2 = l1 + 1;
 401         g = d__[l];
 402         p = (d__[l1] - g) / (e[l] * 2.);
 403         r__ = pythag(&p, &c_b10);
 404         d__[l] = e[l] / (p + d_sign(&r__, &p));
 405         d__[l1] = e[l] * (p + d_sign(&r__, &p));
 406         dl1 = d__[l1];
 407         h__ = g - d__[l];
 408         if (l2 > *n) {
 409             goto L145;
 410         }
 411
 412         i__2 = *n;
 413         for (i__ = l2; i__ <= i__2; ++i__) {
 414             d__[i__] -= h__;
 415         }
 416
 417 L145:
 418         f += h__;
 419         p = d__[m];
 420         c__ = 1.;
 421         c2 = c__;
 422         el1 = e[l1];
 423         s = 0.;
 424         mml = m - l;
 425         i__2 = mml;
 426         for (ii = 1; ii <= i__2; ++ii) {
 427             c3 = c2;
 428             c2 = c__;
 429             s2 = s;
 430             i__ = m - ii;
 431             g = c__ * e[i__];
 432             h__ = c__ * p;
 433             r__ = pythag(&p, &e[i__]);
 434             e[i__ + 1] = s * r__;
 435             s = e[i__] / r__;
 436             c__ = p / r__;
 437             p = c__ * d__[i__] - s * g;
 438             d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]);
 439             i__3 = *n;
 440             for (k = 1; k <= i__3; ++k) {
 441                 h__ = z__[k + (i__ + 1) * z_dim1];
 442                 z__[k + (i__ + 1) * z_dim1] = s * z__[k + i__ * z_dim1] + c__
 443                         * h__;
 444                 z__[k + i__ * z_dim1] = c__ * z__[k + i__ * z_dim1] - s * h__;
 445             }
 446
 447         }
 448
 449         p = -s * s2 * c3 * el1 * e[l] / dl1;
 450         e[l] = s * p;
 451         d__[l] = c__ * p;
 452         tst2 = tst1 + (d__1 = e[l], abs(d__1));
 453         if (tst2 > tst1) {
 454             goto L130;
 455         }
 456 L220:
 457         d__[l] += f;
 458     }
 459     i__1 = *n;
 460     for (ii = 2; ii <= i__1; ++ii) {
 461         i__ = ii - 1;
 462         k = i__;
 463         p = d__[i__];
 464
 465         i__2 = *n;
 466         for (j = ii; j <= i__2; ++j) {
 467             if (d__[j] >= p) {
 468                 goto L260;
 469             }
 470             k = j;
 471             p = d__[j];
 472 L260:
 473             ;
 474         }
 475
 476         if (k == i__) {
 477             goto L300;
 478         }
 479         d__[k] = d__[i__];
 480         d__[i__] = p;
 481
 482         i__2 = *n;
 483         for (j = 1; j <= i__2; ++j) {
 484             p = z__[j + i__ * z_dim1];
 485             z__[j + i__ * z_dim1] = z__[j + k * z_dim1];
 486             z__[j + k * z_dim1] = p;
 487         }
 488
 489 L300:
 490         ;
 491     }
 492
 493     goto L1001;
 494 L1000:
 495     *ierr = l;
 496 L1001:
 497     return 0;
 498 }
 499
 500 void eigen(double **mat,unsigned long n,double *eig)
 501 {
 502   double *trans,*off;
 503   int ierr,i,j,nm=(int)n;
 504
 505   check_alloc(trans=(double*)malloc(sizeof(double)*nm*nm));
 506   check_alloc(off=(double*)malloc(sizeof(double)*nm));
 507
 508   tred2(&nm,&nm,&mat[0][0],eig,off,trans);
 509   tql2(&nm,&nm,eig,off,trans,&ierr);
 510
 511   if (ierr != 0) {
 512     fprintf(stderr,"Non converging eigenvalues! Exiting\n");
 513     exit(EIG2_TOO_MANY_ITERATIONS);
 514   }
 515
 516   for (i=0;i<nm;i++)
 517     for (j=0;j<nm;j++)
 518       mat[i][j]=trans[i+nm*j];
 519
 520   free(trans);
 521   free(off);
 522 }