~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

# Linux Cross ReferenceTina4/src/math/numeric/eigen.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

```  1 /**@(#)Matrix format function to determine all eigen values and vectors of a
2   @(#)given symmetric matrix A (generic format) the resulting eigen values and
3   @(#)coulumn eigen vectors are sorted in ascending eigen value
4  */
5
6 #include <math.h>
7 #include <tina/sys.h>
8 #include <tina/sysfuncs.h>
9 #include <tina/math.h>
10 #include <tina/mathfuncs.h>
11
12 static void tred2(double **a, int n, double *d, double *e);
13 static int tqli(double *d, double *e, int n, double **z);
14
15 int     matrix_eigen_sym(Matrix * A, Vector * Eval, Matrix * Evec)
16 {
17     int     n;
18     int     success;
19     Vector *e;
20     Vector *d;
21     Vector *eval;
22     Matrix *evec;
23     Matrix *matrix_cast_fill();
24     Matrix *matrix_copy_inplace();
25     Vector *vector_alloc();
26
27     int     i, j;
28
29     if (A == NULL || A->n != A->m)
30         return (0);
31
32     n = A->n;
33     if (Evec == NULL || Eval == NULL || Evec->n != n || Evec->m != n || Eval->n != n)
34         return (0);
35
36     evec = matrix_cast_fill(A, double_v);
37     eval = vector_alloc(n, double_v);
38     d = vector_alloc(n, double_v);
39     e = vector_alloc(n, double_v);
40
41     tred2(evec->el.double_v, n, (double *) eval->data, (double *) e->data);
42     success = tqli((double *) eval->data, (double *) e->data, n, evec->el.double_v);
43
44     if (success)
45         for (i = 0; i < n - 1; ++i)
46         {
47             double  min_eval = VECTOR_DOUBLE(eval, i);
48             int     minc = i;
49
50             for (j = i + 1; j < n; ++j)
51                 if (VECTOR_DOUBLE(eval, j) < min_eval)
52                 {
53                     min_eval = VECTOR_DOUBLE(eval, j);
54                     minc = j;
55                 }
56             if (minc != i)
57             {
58                 double  temp = VECTOR_DOUBLE(eval, i);
59
60                 VECTOR_DOUBLE(eval, i) = VECTOR_DOUBLE(eval, minc);
61                 VECTOR_DOUBLE(eval, minc) = temp;
62                 (void) matrix_swap_cols(evec, i, minc);
63             }
64         }
65
66     (void) matrix_copy_inplace(evec, Evec);
67
68     vector_copy_inplace(Eval, eval);    /** copy goes <--- **/
69
70     matrix_free((Matrix *) evec);
71     vector_free(e);
72     vector_free(d);
73     vector_free(eval);
74     return (success);
75 }
76
77 /**
78 eigenvalues and eigenvectors of symmetric  n  by n matrix A:
79 Eval = already allocated vector of n doubles
80 Evec = already allocated vector of n pointers to vectors of n doubles
81 VECTOR_DOUBLE(Eval,0)          = biggest eigenvalue
82 (Vector *) VECTOR_PTR(Evec, 0) = associated eigenvector
83 **/
84 int     matrix_eigen(Matrix * A, Vector * Eval, Vector * Evec)
85 {
86     int     n;
87     int     success;
88     Vector *e;
89     Vector *d;
90     Vector *eval;
91     Matrix *evec;
92     Matrix *matrix_cast_fill();
93     Vector *vector_alloc();
94
95     int     i, j;
96
97     if (A == NULL || A->n != A->m)
98         return (0);
99
100     n = A->n;
101     if (Evec == NULL || Eval == NULL)
102         return (0);
103
104     evec = matrix_cast_fill(A, double_v);
105     eval = vector_alloc(n, double_v);
106     d = vector_alloc(n, double_v);
107     e = vector_alloc(n, double_v);
108
109     tred2(evec->el.double_v, n, (double *) eval->data, (double *) e->data);
110     success = tqli((double *) eval->data, (double *) e->data, n, evec->el.double_v);
111
112     if (success)
113         for (i = 0; i < n - 1; ++i)
114         {
115             double  min_eval = VECTOR_DOUBLE(eval, i);
116             int     minc = i;
117
118             for (j = i + 1; j < n; ++j)
119                 if (VECTOR_DOUBLE(eval, j) < min_eval)
120                 {
121                     min_eval = VECTOR_DOUBLE(eval, j);
122                     minc = j;
123                 }
124             if (minc != i)
125             {
126                 double  temp = VECTOR_DOUBLE(eval, i);
127
128                 VECTOR_DOUBLE(eval, i) = VECTOR_DOUBLE(eval, minc);
129                 VECTOR_DOUBLE(eval, minc) = temp;
130                 (void) matrix_swap_cols(evec, i, minc);
131             }
132         }
133
134     /** largest first **/
135     for (j = 0; j < n; j++)
136     {
137         Vector *ej = VECTOR_PTR(Evec, j);
138
139         VECTOR_DOUBLE(Eval, j) = VECTOR_DOUBLE(eval, n - 1 - j);
140         for (i = 0; i < n; i++)
141             VECTOR_DOUBLE(ej, i) = evec->el.double_v[i][n - 1 - j];
142     }
143
144
145     matrix_free((Matrix *) evec);
146     vector_free(e);
147     vector_free(d);
148     vector_free(eval);
149     return (success);
150 }
151
152 static void tred2(double **a, int n, double *d, double *e)
153 {
154     int     l, k, j, i;
155     double  scale, hh, h, g, f;
156
157     for (i = 0; i < n; ++i)
158         (a[i])--;
159     a--;
160     d--;
161     e--;
162
163     for (i = n; i >= 2; i--)
164     {
165         l = i - 1;
166         h = scale = 0.0;
167         if (l > 1)
168         {
169             for (k = 1; k <= l; k++)
170                 scale += fabs(a[i][k]);
171             if (scale == 0.0)
172                 e[i] = a[i][l];
173             else
174             {
175                 for (k = 1; k <= l; k++)
176                 {
177                     a[i][k] /= scale;
178                     h += a[i][k] * a[i][k];
179                 }
180                 f = a[i][l];
181                 g = f > 0 ? -sqrt(h) : sqrt(h);
182                 e[i] = scale * g;
183                 h -= f * g;
184                 a[i][l] = f - g;
185                 f = 0.0;
186                 for (j = 1; j <= l; j++)
187                 {
188                     a[j][i] = a[i][j] / h;
189                     g = 0.0;
190                     for (k = 1; k <= j; k++)
191                         g += a[j][k] * a[i][k];
192                     for (k = j + 1; k <= l; k++)
193                         g += a[k][j] * a[i][k];
194                     e[j] = g / h;
195                     f += e[j] * a[i][j];
196                 }
197                 hh = f / (h + h);
198                 for (j = 1; j <= l; j++)
199                 {
200                     f = a[i][j];
201                     e[j] = g = e[j] - hh * f;
202                     for (k = 1; k <= j; k++)
203                         a[j][k] -= (f * e[k] + g * a[i][k]);
204                 }
205             }
206         } else
207             e[i] = a[i][l];
208         d[i] = h;
209     }
210     d[1] = 0.0;
211     e[1] = 0.0;
212     for (i = 1; i <= n; i++)
213     {
214         l = i - 1;
215         if (d[i])
216         {
217             for (j = 1; j <= l; j++)
218             {
219                 g = 0.0;
220                 for (k = 1; k <= l; k++)
221                     g += a[i][k] * a[k][j];
222                 for (k = 1; k <= l; k++)
223                     a[k][j] -= g * a[k][i];
224             }
225         }
226         d[i] = a[i][i];
227         a[i][i] = 1.0;
228         for (j = 1; j <= l; j++)
229             a[j][i] = a[i][j] = 0.0;
230     }
231
232     a++;
233     d++;
234     e++;
235     for (i = 0; i < n; ++i)
236         (a[i])++;
237 }
238
239 #undef SIGN
240
241 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
242
243 static int tqli(double *d, double *e, int n, double **z)
244 {
245     int     m, l, iter, i, k;
246     double  s, r, p, g, f, dd, c, b;
247
248     for (i = 0; i < n; ++i)
249         (z[i])--;
250     z--;
251     d--;
252     e--;
253
254     for (i = 2; i <= n; i++)
255         e[i - 1] = e[i];
256     e[n] = 0.0;
257     for (l = 1; l <= n; l++)
258     {
259         iter = 0;
260         do
261         {
262             for (m = l; m <= n - 1; m++)
263             {
264                 dd = fabs(d[m]) + fabs(d[m + 1]);
265                 if (fabs(e[m]) + dd == dd)
266                     break;
267             }
268             if (m != l)
269             {
270                 if (iter++ == 30)
271                     return (0);
272                 g = (d[l + 1] - d[l]) / (2.0 * e[l]);
273                 r = sqrt((g * g) + 1.0);
274                 g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
275                 s = c = 1.0;
276                 p = 0.0;
277                 for (i = m - 1; i >= l; i--)
278                 {
279                     f = s * e[i];
280                     b = c * e[i];
281                     if (fabs(f) >= fabs(g))
282                     {
283                         c = g / f;
284                         r = sqrt((c * c) + 1.0);
285                         e[i + 1] = f * r;
286                         c *= (s = 1.0 / r);
287                     } else
288                     {
289                         s = f / g;
290                         r = sqrt((s * s) + 1.0);
291                         e[i + 1] = g * r;
292                         s *= (c = 1.0 / r);
293                     }
294                     g = d[i + 1] - p;
295                     r = (d[i] - g) * s + 2.0 * c * b;
296                     p = s * r;
297                     d[i + 1] = g + p;
298                     g = c * r - b;
299                     for (k = 1; k <= n; k++)
300                     {
301                         f = z[k][i + 1];
302                         z[k][i + 1] = s * z[k][i] + c * f;
303                         z[k][i] = c * z[k][i] - s * f;
304                     }
305                 }
306                 d[l] = d[l] - p;
307                 e[l] = g;
308                 e[m] = 0.0;
309             }
310         } while (m != l);
311     }
312
313     z++;
314     d++;
315     e++;
316     for (i = 0; i < n; ++i)
317         (z[i])++;
318
319     return (1);
320 }
321 ```

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~