4 #include "tisean_cec.h"
6 typedef double doublereal;
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))
13 static doublereal c_b10 = 1.;
15 extern void check_alloc(void*);
17 double d_sign(double *a,double *b)
20 x = (*a >= 0 ? *a : - *a);
21 return ( *b >= 0 ? x : -x);
24 doublereal pythag(doublereal *a, doublereal *b)
26 doublereal ret_val, d__1, d__2, d__3;
27 static doublereal p, r__, s, t, u;
29 d__1 = abs(*a), d__2 = abs(*b);
34 d__2 = abs(*a), d__3 = abs(*b);
35 d__1 = min(d__2,d__3) / p;
46 r__ = d__1 * d__1 * r__;
54 int tred2(integer *nm, integer *n, doublereal *a,
55 doublereal *d__, doublereal *e, doublereal *z__)
57 integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
60 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
62 static doublereal f, g, h__;
63 static integer i__, j, k, l;
65 static integer ii, jp1;
66 static doublereal scale;
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). */
74 /* this subroutine reduces a real symmetric matrix to a */
75 /* symmetric tridiagonal matrix using and accumulating */
76 /* orthogonal similarity transformations. */
80 /* nm must be set to the row dimension of two-dimensional */
81 /* array parameters as declared in the calling program */
82 /* dimension statement. */
84 /* n is the order of the matrix. */
86 /* a contains the real symmetric input matrix. only the */
87 /* lower triangle of the matrix need be supplied. */
91 /* d contains the diagonal elements of the tridiagonal matrix. */
93 /* e contains the subdiagonal elements of the tridiagonal */
94 /* matrix in its last n-1 positions. e(1) is set to zero. */
96 /* z contains the orthogonal transformation matrix */
97 /* produced in the reduction. */
99 /* a and z may coincide. if distinct, a is unaltered. */
101 /* questions and comments should be directed to burton s. garbow, */
102 /* mathematics and computer science div, argonne national laboratory */
104 /* this version dated august 1983. */
106 /* ------------------------------------------------------------------ */
109 z_offset = 1 + z_dim1 * 1;
114 a_offset = 1 + a_dim1 * 1;
118 for (i__ = 1; i__ <= i__1; ++i__) {
121 for (j = i__; j <= i__2; ++j) {
122 z__[j + i__ * z_dim1] = a[j + i__ * a_dim1];
125 d__[i__] = a[*n + i__ * a_dim1];
132 for (ii = 2; ii <= i__1; ++ii) {
141 for (k = 1; k <= i__2; ++k) {
142 scale += (d__1 = d__[k], abs(d__1));
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.;
162 for (k = 1; k <= i__2; ++k) {
164 h__ += d__[k] * d__[k];
169 g = -d_sign(&d__1, &f);
174 for (j = 1; j <= i__2; ++j) {
179 for (j = 1; j <= i__2; ++j) {
181 z__[j + i__ * z_dim1] = f;
182 g = e[j] + z__[j + j * z_dim1] * f;
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;
200 for (j = 1; j <= i__2; ++j) {
205 hh = f / (h__ + h__);
207 for (j = 1; j <= i__2; ++j) {
211 for (j = 1; j <= i__2; ++j) {
216 for (k = j; k <= i__3; ++k) {
217 z__[k + j * z_dim1] = z__[k + j * z_dim1] - f * e[k] - g *
221 d__[j] = z__[l + j * z_dim1];
222 z__[i__ + j * z_dim1] = 0.;
229 for (i__ = 2; i__ <= i__1; ++i__) {
231 z__[*n + l * z_dim1] = z__[l + l * z_dim1];
232 z__[l + l * z_dim1] = 1.;
239 for (k = 1; k <= i__2; ++k) {
240 d__[k] = z__[k + i__ * z_dim1] / h__;
244 for (j = 1; j <= i__2; ++j) {
248 for (k = 1; k <= i__3; ++k) {
249 g += z__[k + i__ * z_dim1] * z__[k + j * z_dim1];
253 for (k = 1; k <= i__3; ++k) {
254 z__[k + j * z_dim1] -= g * d__[k];
260 for (k = 1; k <= i__3; ++k) {
261 z__[k + i__ * z_dim1] = 0.;
268 for (i__ = 1; i__ <= i__1; ++i__) {
269 d__[i__] = z__[*n + i__ * z_dim1];
270 z__[*n + i__ * z_dim1] = 0.;
273 z__[*n + *n * z_dim1] = 1.;
278 int tql2(integer *nm, integer *n, doublereal *d__,
279 doublereal *e, doublereal *z__, integer *ierr)
281 integer z_dim1, z_offset, i__1, i__2, i__3;
282 doublereal d__1, d__2;
284 double d_sign(doublereal *, doublereal *);
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;
292 static doublereal dl1, el1;
294 static doublereal tst1, tst2;
295 extern doublereal pythag_(doublereal *, doublereal *);
299 /* this subroutine is a translation of the algol procedure tql2, */
300 /* num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and */
302 /* handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). */
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. */
312 /* nm must be set to the row dimension of two-dimensional */
313 /* array parameters as declared in the calling program */
314 /* dimension statement. */
316 /* n is the order of the matrix. */
318 /* d contains the diagonal elements of the input matrix. */
320 /* e contains the subdiagonal elements of the input matrix */
321 /* in its last n-1 positions. e(1) is arbitrary. */
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. */
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. */
334 /* e has been destroyed. */
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 */
342 /* zero for normal return, */
343 /* j if the j-th eigenvalue has not been */
344 /* determined after 30 iterations. */
346 /* calls pythag for dsqrt(a*a + b*b) . */
348 /* questions and comments should be directed to burton s. garbow, */
349 /* mathematics and computer science div, argonne national laboratory */
351 /* this version dated august 1983. */
353 /* ------------------------------------------------------------------ */
356 z_offset = 1 + z_dim1 * 1;
367 for (i__ = 2; i__ <= i__1; ++i__) {
376 for (l = 1; l <= i__1; ++l) {
378 h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2));
383 for (m = l; m <= i__2; ++m) {
384 tst2 = tst1 + (d__1 = e[m], abs(d__1));
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));
413 for (i__ = l2; i__ <= i__2; ++i__) {
426 for (ii = 1; ii <= i__2; ++ii) {
433 r__ = pythag(&p, &e[i__]);
434 e[i__ + 1] = s * r__;
437 p = c__ * d__[i__] - s * g;
438 d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]);
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__
444 z__[k + i__ * z_dim1] = c__ * z__[k + i__ * z_dim1] - s * h__;
449 p = -s * s2 * c3 * el1 * e[l] / dl1;
452 tst2 = tst1 + (d__1 = e[l], abs(d__1));
460 for (ii = 2; ii <= i__1; ++ii) {
466 for (j = ii; j <= i__2; ++j) {
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;
500 void eigen(double **mat,unsigned long n,double *eig)
503 int ierr,i,j,nm=(int)n;
505 check_alloc(trans=(double*)malloc(sizeof(double)*nm*nm));
506 check_alloc(off=(double*)malloc(sizeof(double)*nm));
508 tred2(&nm,&nm,&mat[0][0],eig,off,trans);
509 tql2(&nm,&nm,eig,off,trans,&ierr);
512 fprintf(stderr,"Non converging eigenvalues! Exiting\n");
513 exit(EIG2_TOO_MANY_ITERATIONS);
518 mat[i][j]=trans[i+nm*j];