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

Linux Cross Reference
Tina4/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 ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.