#include #include #include #include "tisean_cec.h" typedef double doublereal; typedef int integer; #define abs(x) (((x)>=0.0)?(x):-(x)) #define min(x,y) (((x)<=(y))?(x):(y)) #define max(x,y) (((x)>=(y))?(x):(y)) static doublereal c_b10 = 1.; extern void check_alloc(void*); double d_sign(double *a,double *b) { double x; x = (*a >= 0 ? *a : - *a); return ( *b >= 0 ? x : -x); } doublereal pythag(doublereal *a, doublereal *b) { doublereal ret_val, d__1, d__2, d__3; static doublereal p, r__, s, t, u; d__1 = abs(*a), d__2 = abs(*b); p = max(d__1,d__2); if (p == 0.) { goto L20; } d__2 = abs(*a), d__3 = abs(*b); d__1 = min(d__2,d__3) / p; r__ = d__1 * d__1; L10: t = r__ + 4.; if (t == 4.) { goto L20; } s = r__ / t; u = s * 2. + 1.; p = u * p; d__1 = s / u; r__ = d__1 * d__1 * r__; goto L10; L20: ret_val = p; return ret_val; } int tred2(integer *nm, integer *n, doublereal *a, doublereal *d__, doublereal *e, doublereal *z__) { integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3; doublereal d__1; double sqrt(doublereal), d_sign(doublereal *, doublereal *); static doublereal f, g, h__; static integer i__, j, k, l; static doublereal hh; static integer ii, jp1; static doublereal scale; /* this subroutine is a translation of the algol procedure tred2, */ /* num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson. */ /* handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */ /* this subroutine reduces a real symmetric matrix to a */ /* symmetric tridiagonal matrix using and accumulating */ /* orthogonal similarity transformations. */ /* on input */ /* nm must be set to the row dimension of two-dimensional */ /* array parameters as declared in the calling program */ /* dimension statement. */ /* n is the order of the matrix. */ /* a contains the real symmetric input matrix. only the */ /* lower triangle of the matrix need be supplied. */ /* on output */ /* d contains the diagonal elements of the tridiagonal matrix. */ /* e contains the subdiagonal elements of the tridiagonal */ /* matrix in its last n-1 positions. e(1) is set to zero. */ /* z contains the orthogonal transformation matrix */ /* produced in the reduction. */ /* a and z may coincide. if distinct, a is unaltered. */ /* questions and comments should be directed to burton s. garbow, */ /* mathematics and computer science div, argonne national laboratory */ /* this version dated august 1983. */ /* ------------------------------------------------------------------ */ z_dim1 = *nm; z_offset = 1 + z_dim1 * 1; z__ -= z_offset; --e; --d__; a_dim1 = *nm; a_offset = 1 + a_dim1 * 1; a -= a_offset; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *n; for (j = i__; j <= i__2; ++j) { z__[j + i__ * z_dim1] = a[j + i__ * a_dim1]; } d__[i__] = a[*n + i__ * a_dim1]; } if (*n == 1) { goto L510; } i__1 = *n; for (ii = 2; ii <= i__1; ++ii) { i__ = *n + 2 - ii; l = i__ - 1; h__ = 0.; scale = 0.; if (l < 2) { goto L130; } i__2 = l; for (k = 1; k <= i__2; ++k) { scale += (d__1 = d__[k], abs(d__1)); } if (scale != 0.) { goto L140; } L130: e[i__] = d__[l]; i__2 = l; for (j = 1; j <= i__2; ++j) { d__[j] = z__[l + j * z_dim1]; z__[i__ + j * z_dim1] = 0.; z__[j + i__ * z_dim1] = 0.; } goto L290; L140: i__2 = l; for (k = 1; k <= i__2; ++k) { d__[k] /= scale; h__ += d__[k] * d__[k]; } f = d__[l]; d__1 = sqrt(h__); g = -d_sign(&d__1, &f); e[i__] = scale * g; h__ -= f * g; d__[l] = f - g; i__2 = l; for (j = 1; j <= i__2; ++j) { e[j] = 0.; } i__2 = l; for (j = 1; j <= i__2; ++j) { f = d__[j]; z__[j + i__ * z_dim1] = f; g = e[j] + z__[j + j * z_dim1] * f; jp1 = j + 1; if (l < jp1) { goto L220; } i__3 = l; for (k = jp1; k <= i__3; ++k) { g += z__[k + j * z_dim1] * d__[k]; e[k] += z__[k + j * z_dim1] * f; } L220: e[j] = g; } f = 0.; i__2 = l; for (j = 1; j <= i__2; ++j) { e[j] /= h__; f += e[j] * d__[j]; } hh = f / (h__ + h__); i__2 = l; for (j = 1; j <= i__2; ++j) { e[j] -= hh * d__[j]; } i__2 = l; for (j = 1; j <= i__2; ++j) { f = d__[j]; g = e[j]; i__3 = l; for (k = j; k <= i__3; ++k) { z__[k + j * z_dim1] = z__[k + j * z_dim1] - f * e[k] - g * d__[k]; } d__[j] = z__[l + j * z_dim1]; z__[i__ + j * z_dim1] = 0.; } L290: d__[i__] = h__; } i__1 = *n; for (i__ = 2; i__ <= i__1; ++i__) { l = i__ - 1; z__[*n + l * z_dim1] = z__[l + l * z_dim1]; z__[l + l * z_dim1] = 1.; h__ = d__[i__]; if (h__ == 0.) { goto L380; } i__2 = l; for (k = 1; k <= i__2; ++k) { d__[k] = z__[k + i__ * z_dim1] / h__; } i__2 = l; for (j = 1; j <= i__2; ++j) { g = 0.; i__3 = l; for (k = 1; k <= i__3; ++k) { g += z__[k + i__ * z_dim1] * z__[k + j * z_dim1]; } i__3 = l; for (k = 1; k <= i__3; ++k) { z__[k + j * z_dim1] -= g * d__[k]; } } L380: i__3 = l; for (k = 1; k <= i__3; ++k) { z__[k + i__ * z_dim1] = 0.; } } L510: i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { d__[i__] = z__[*n + i__ * z_dim1]; z__[*n + i__ * z_dim1] = 0.; } z__[*n + *n * z_dim1] = 1.; e[1] = 0.; return 0; } int tql2(integer *nm, integer *n, doublereal *d__, doublereal *e, doublereal *z__, integer *ierr) { integer z_dim1, z_offset, i__1, i__2, i__3; doublereal d__1, d__2; double d_sign(doublereal *, doublereal *); static doublereal c__, f, g, h__; static integer i__, j, k, l, m; static doublereal p, r__, s, c2, c3; static integer l1, l2; static doublereal s2; static integer ii; static doublereal dl1, el1; static integer mml; static doublereal tst1, tst2; extern doublereal pythag_(doublereal *, doublereal *); /* this subroutine is a translation of the algol procedure tql2, */ /* num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and */ /* wilkinson. */ /* handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). */ /* this subroutine finds the eigenvalues and eigenvectors */ /* of a symmetric tridiagonal matrix by the ql method. */ /* the eigenvectors of a full symmetric matrix can also */ /* be found if tred2 has been used to reduce this */ /* full matrix to tridiagonal form. */ /* on input */ /* nm must be set to the row dimension of two-dimensional */ /* array parameters as declared in the calling program */ /* dimension statement. */ /* n is the order of the matrix. */ /* d contains the diagonal elements of the input matrix. */ /* e contains the subdiagonal elements of the input matrix */ /* in its last n-1 positions. e(1) is arbitrary. */ /* z contains the transformation matrix produced in the */ /* reduction by tred2, if performed. if the eigenvectors */ /* of the tridiagonal matrix are desired, z must contain */ /* the identity matrix. */ /* on output */ /* d contains the eigenvalues in ascending order. if an */ /* error exit is made, the eigenvalues are correct but */ /* unordered for indices 1,2,...,ierr-1. */ /* e has been destroyed. */ /* z contains orthonormal eigenvectors of the symmetric */ /* tridiagonal (or full) matrix. if an error exit is made, */ /* z contains the eigenvectors associated with the stored */ /* eigenvalues. */ /* ierr is set to */ /* zero for normal return, */ /* j if the j-th eigenvalue has not been */ /* determined after 30 iterations. */ /* calls pythag for dsqrt(a*a + b*b) . */ /* questions and comments should be directed to burton s. garbow, */ /* mathematics and computer science div, argonne national laboratory */ /* this version dated august 1983. */ /* ------------------------------------------------------------------ */ z_dim1 = *nm; z_offset = 1 + z_dim1 * 1; z__ -= z_offset; --e; --d__; *ierr = 0; if (*n == 1) { goto L1001; } i__1 = *n; for (i__ = 2; i__ <= i__1; ++i__) { e[i__ - 1] = e[i__]; } f = 0.; tst1 = 0.; e[*n] = 0.; i__1 = *n; for (l = 1; l <= i__1; ++l) { j = 0; h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2)); if (tst1 < h__) { tst1 = h__; } i__2 = *n; for (m = l; m <= i__2; ++m) { tst2 = tst1 + (d__1 = e[m], abs(d__1)); if (tst2 == tst1) { goto L120; } } L120: if (m == l) { goto L220; } L130: if (j == 30) { goto L1000; } ++j; l1 = l + 1; l2 = l1 + 1; g = d__[l]; p = (d__[l1] - g) / (e[l] * 2.); r__ = pythag(&p, &c_b10); d__[l] = e[l] / (p + d_sign(&r__, &p)); d__[l1] = e[l] * (p + d_sign(&r__, &p)); dl1 = d__[l1]; h__ = g - d__[l]; if (l2 > *n) { goto L145; } i__2 = *n; for (i__ = l2; i__ <= i__2; ++i__) { d__[i__] -= h__; } L145: f += h__; p = d__[m]; c__ = 1.; c2 = c__; el1 = e[l1]; s = 0.; mml = m - l; i__2 = mml; for (ii = 1; ii <= i__2; ++ii) { c3 = c2; c2 = c__; s2 = s; i__ = m - ii; g = c__ * e[i__]; h__ = c__ * p; r__ = pythag(&p, &e[i__]); e[i__ + 1] = s * r__; s = e[i__] / r__; c__ = p / r__; p = c__ * d__[i__] - s * g; d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]); i__3 = *n; for (k = 1; k <= i__3; ++k) { h__ = z__[k + (i__ + 1) * z_dim1]; z__[k + (i__ + 1) * z_dim1] = s * z__[k + i__ * z_dim1] + c__ * h__; z__[k + i__ * z_dim1] = c__ * z__[k + i__ * z_dim1] - s * h__; } } p = -s * s2 * c3 * el1 * e[l] / dl1; e[l] = s * p; d__[l] = c__ * p; tst2 = tst1 + (d__1 = e[l], abs(d__1)); if (tst2 > tst1) { goto L130; } L220: d__[l] += f; } i__1 = *n; for (ii = 2; ii <= i__1; ++ii) { i__ = ii - 1; k = i__; p = d__[i__]; i__2 = *n; for (j = ii; j <= i__2; ++j) { if (d__[j] >= p) { goto L260; } k = j; p = d__[j]; L260: ; } if (k == i__) { goto L300; } d__[k] = d__[i__]; d__[i__] = p; i__2 = *n; for (j = 1; j <= i__2; ++j) { p = z__[j + i__ * z_dim1]; z__[j + i__ * z_dim1] = z__[j + k * z_dim1]; z__[j + k * z_dim1] = p; } L300: ; } goto L1001; L1000: *ierr = l; L1001: return 0; } void eigen(double **mat,unsigned long n,double *eig) { double *trans,*off; int ierr,i,j,nm=(int)n; check_alloc(trans=(double*)malloc(sizeof(double)*nm*nm)); check_alloc(off=(double*)malloc(sizeof(double)*nm)); tred2(&nm,&nm,&mat[0][0],eig,off,trans); tql2(&nm,&nm,eig,off,trans,&ierr); if (ierr != 0) { fprintf(stderr,"Non converging eigenvalues! Exiting\n"); exit(EIG2_TOO_MANY_ITERATIONS); } for (i=0;i