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

Linux Cross Reference
Tina5/tina-libs/tina/math/mathNum_eigen.c

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

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU Lesser General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU Lesser General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU Lesser General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  **********
 20  *
 21  * Program :    TINA
 22  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathNum_eigen.c,v $
 23  * Date    :  $Date: 2003/10/06 12:29:47 $
 24  * Version :  $Revision: 1.5 $
 25  * CVS Id  :  $Id: mathNum_eigen.c,v 1.5 2003/10/06 12:29:47 neil Exp $
 26  *
 27  * Author  :  Legacy TINA
 28  *
 29  * Notes :Matrix format function to determine all eigen values and vectors of a
 30  * given symmetric matrix A (generic format) the resulting eigen values and
 31  * coulumn eigen vectors are sorted in ascending eigen value
 32  *
 33  *********
 34 */
 35 /** 
 36  *  @file
 37  *  @brief  Matrix format function to determine all eigen values and vectors of a
 38  *  given symmetric matrix A.   
 39  *
 40  *  The resulting eigenvalues and vectors are sorted in ascending eigen value.
 41  *
 42  *  See Numerical Recipes (3rd ed) Ch 11 pp456-495 
 43  *  for more details on the basic workings of the equivalent tred2 and tqli functions.
 44  * 
 45 */
 46 
 47 #include "mathNum_eigen.h"
 48 
 49 #if HAVE_CONFIG_H
 50 #include <config.h>
 51 #endif
 52 
 53 #include <math.h>
 54 #include <tina/sys/sysDef.h>
 55 #include <tina/sys/sysPro.h>
 56 #include <tina/math/math_MatrDef.h>
 57 #include <tina/math/math_MatrPro.h>
 58 #include <tina/math/math_VecDef.h>
 59 #include <tina/math/math_VecPro.h>
 60 
 61 static void     tred2(double **a, int n, double *d, double *e);
 62 static int      tqli(double *d, double *e, int n, double **z);
 63 
 64 int             matrix_eigen_sym(Matrix * A, Vector * Eval, Matrix * Evec)
 65 {
 66         int             n;
 67         int             success;
 68         Vector         *e;
 69         Vector         *d;
 70         Vector         *eval;
 71         Matrix         *evec;
 72         Matrix         *matrix_cast_fill();
 73         Matrix         *matrix_copy_inplace();
 74         Vector         *vector_alloc();
 75 
 76         int             i, j;
 77 
 78         if (A == NULL || A->n != A->m)
 79                 return (0);
 80 
 81         n = A->n;
 82         if (Evec == NULL || Eval == NULL || Evec->n != n || Evec->m != n || Eval->n != n)
 83                 return (0);
 84 
 85         evec = matrix_cast_fill(A, double_v);
 86         eval = vector_alloc(n, double_v);
 87         d = vector_alloc(n, double_v);
 88         e = vector_alloc(n, double_v);
 89 
 90         tred2(evec->el.double_v, n, (double *) eval->data, (double *) e->data);
 91         success = tqli((double *) eval->data, (double *) e->data, n, evec->el.double_v);
 92 
 93         if (success)
 94                 for (i = 0; i < n - 1; ++i)
 95                 {
 96                         double          min_eval = VECTOR_DOUBLE(eval, i);
 97                         int             minc = i;
 98 
 99                         for (j = i + 1; j < n; ++j)
100                                 if (VECTOR_DOUBLE(eval, j) < min_eval)
101                                 {
102                                         min_eval = VECTOR_DOUBLE(eval, j);
103                                         minc = j;
104                                 }
105                         if (minc != i)
106                         {
107                                 double          temp = VECTOR_DOUBLE(eval, i);
108 
109                                 VECTOR_DOUBLE(eval, i) = VECTOR_DOUBLE(eval, minc);
110                                 VECTOR_DOUBLE(eval, minc) = temp;
111                                 (void) matrix_swap_cols(evec, i, minc);
112                         }
113                 }
114 
115         (void) matrix_copy_inplace(evec, Evec);
116 
117         vector_copy_inplace(Eval, eval);        /** copy goes <--- **/
118 
119         matrix_free((Matrix *) evec);
120         vector_free(e);
121         vector_free(d);
122         vector_free(eval);
123         return (success);
124 }
125 
126 /**
127 eigenvalues and eigenvectors of symmetric  n  by n matrix A:
128 Eval = already allocated vector of n doubles
129 Evec = already allocated vector of n pointers to vectors of n doubles
130 VECTOR_DOUBLE(Eval,0)          = biggest eigenvalue
131 (Vector *) VECTOR_PTR(Evec, 0) = associated eigenvector
132 **/
133 int             matrix_eigen(Matrix * A, Vector * Eval, Vector * Evec)
134 {
135         int             n;
136         int             success;
137         Vector         *e;
138         Vector         *d;
139         Vector         *eval;
140         Matrix         *evec;
141         Matrix         *matrix_cast_fill();
142         Vector         *vector_alloc();
143 
144         int             i, j;
145 
146         if (A == NULL || A->n != A->m)
147                 return (0);
148 
149         n = A->n;
150         if (Evec == NULL || Eval == NULL)
151                 return (0);
152 
153         evec = matrix_cast_fill(A, double_v);
154         eval = vector_alloc(n, double_v);
155         d = vector_alloc(n, double_v);
156         e = vector_alloc(n, double_v);
157 
158         tred2(evec->el.double_v, n, (double *) eval->data, (double *) e->data);
159         success = tqli((double *) eval->data, (double *) e->data, n, evec->el.double_v);
160 
161         if (success)
162                 for (i = 0; i < n - 1; ++i)
163                 {
164                         double          min_eval = VECTOR_DOUBLE(eval, i);
165                         int             minc = i;
166 
167                         for (j = i + 1; j < n; ++j)
168                                 if (VECTOR_DOUBLE(eval, j) < min_eval)
169                                 {
170                                         min_eval = VECTOR_DOUBLE(eval, j);
171                                         minc = j;
172                                 }
173                         if (minc != i)
174                         {
175                                 double          temp = VECTOR_DOUBLE(eval, i);
176 
177                                 VECTOR_DOUBLE(eval, i) = VECTOR_DOUBLE(eval, minc);
178                                 VECTOR_DOUBLE(eval, minc) = temp;
179                                 (void) matrix_swap_cols(evec, i, minc);
180                         }
181                 }
182 
183         /** largest first **/
184         for (j = 0; j < n; j++)
185         {
186                 Vector         *ej = VECTOR_PTR(Evec, j);
187 
188                 VECTOR_DOUBLE(Eval, j) = VECTOR_DOUBLE(eval, n - 1 - j);
189                 for (i = 0; i < n; i++)
190                         VECTOR_DOUBLE(ej, i) = evec->el.double_v[i][n - 1 - j];
191         }
192 
193 
194         matrix_free((Matrix *) evec);
195         vector_free(e);
196         vector_free(d);
197         vector_free(eval);
198         return (success);
199 }
200 
201 static void     tred2(double **a, int n, double *d, double *e)
202 {
203         int             l, k, j, i;
204         double          scale, hh, h, g, f;
205 
206         for (i = 0; i < n; ++i)
207                 (a[i])--;
208         a--;
209         d--;
210         e--;
211 
212         for (i = n; i >= 2; i--)
213         {
214                 l = i - 1;
215                 h = scale = 0.0;
216                 if (l > 1)
217                 {
218                         for (k = 1; k <= l; k++)
219                                 scale += fabs(a[i][k]);
220                         if (scale == 0.0)
221                                 e[i] = a[i][l];
222                         else
223                         {
224                                 for (k = 1; k <= l; k++)
225                                 {
226                                         a[i][k] /= scale;
227                                         h += a[i][k] * a[i][k];
228                                 }
229                                 f = a[i][l];
230                                 g = f > 0 ? -sqrt(h) : sqrt(h);
231                                 e[i] = scale * g;
232                                 h -= f * g;
233                                 a[i][l] = f - g;
234                                 f = 0.0;
235                                 for (j = 1; j <= l; j++)
236                                 {
237                                         a[j][i] = a[i][j] / h;
238                                         g = 0.0;
239                                         for (k = 1; k <= j; k++)
240                                                 g += a[j][k] * a[i][k];
241                                         for (k = j + 1; k <= l; k++)
242                                                 g += a[k][j] * a[i][k];
243                                         e[j] = g / h;
244                                         f += e[j] * a[i][j];
245                                 }
246                                 hh = f / (h + h);
247                                 for (j = 1; j <= l; j++)
248                                 {
249                                         f = a[i][j];
250                                         e[j] = g = e[j] - hh * f;
251                                         for (k = 1; k <= j; k++)
252                                                 a[j][k] -= (f * e[k] + g * a[i][k]);
253                                 }
254                         }
255                 } else
256                         e[i] = a[i][l];
257                 d[i] = h;
258         }
259         d[1] = 0.0;
260         e[1] = 0.0;
261         for (i = 1; i <= n; i++)
262         {
263                 l = i - 1;
264                 if (d[i])
265                 {
266                         for (j = 1; j <= l; j++)
267                         {
268                                 g = 0.0;
269                                 for (k = 1; k <= l; k++)
270                                         g += a[i][k] * a[k][j];
271                                 for (k = 1; k <= l; k++)
272                                         a[k][j] -= g * a[k][i];
273                         }
274                 }
275                 d[i] = a[i][i];
276                 a[i][i] = 1.0;
277                 for (j = 1; j <= l; j++)
278                         a[j][i] = a[i][j] = 0.0;
279         }
280 
281         a++;
282         d++;
283         e++;
284         for (i = 0; i < n; ++i)
285                 (a[i])++;
286 }
287 
288 #undef SIGN
289 
290 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
291 
292 static int      tqli(double *d, double *e, int n, double **z)
293 {
294         int             m, l, iter, i, k;
295         double          s, r, p, g, f, dd, c, b;
296 
297         for (i = 0; i < n; ++i)
298                 (z[i])--;
299         z--;
300         d--;
301         e--;
302 
303         for (i = 2; i <= n; i++)
304                 e[i - 1] = e[i];
305         e[n] = 0.0;
306         for (l = 1; l <= n; l++)
307         {
308                 iter = 0;
309                 do
310                 {
311                         for (m = l; m <= n - 1; m++)
312                         {
313                                 dd = fabs(d[m]) + fabs(d[m + 1]);
314                                 if (fabs(e[m]) + dd == dd)
315                                         break;
316                         }
317                         if (m != l)
318                         {
319                                 if (iter++ == 30)
320                                         return (0);
321                                 g = (d[l + 1] - d[l]) / (2.0 * e[l]);
322                                 r = sqrt((g * g) + 1.0);
323                                 g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
324                                 s = c = 1.0;
325                                 p = 0.0;
326                                 for (i = m - 1; i >= l; i--)
327                                 {
328                                         f = s * e[i];
329                                         b = c * e[i];
330                                         if (fabs(f) >= fabs(g))
331                                         {
332                                                 c = g / f;
333                                                 r = sqrt((c * c) + 1.0);
334                                                 e[i + 1] = f * r;
335                                                 c *= (s = 1.0 / r);
336                                         } else
337                                         {
338                                                 s = f / g;
339                                                 r = sqrt((s * s) + 1.0);
340                                                 e[i + 1] = g * r;
341                                                 s *= (c = 1.0 / r);
342                                         }
343                                         g = d[i + 1] - p;
344                                         r = (d[i] - g) * s + 2.0 * c * b;
345                                         p = s * r;
346                                         d[i + 1] = g + p;
347                                         g = c * r - b;
348                                         for (k = 1; k <= n; k++)
349                                         {
350                                                 f = z[k][i + 1];
351                                                 z[k][i + 1] = s * z[k][i] + c * f;
352                                                 z[k][i] = c * z[k][i] - s * f;
353                                         }
354                                 }
355                                 d[l] = d[l] - p;
356                                 e[l] = g;
357                                 e[m] = 0.0;
358                         }
359                 } while (m != l);
360         }
361 
362         z++;
363         d++;
364         e++;
365         for (i = 0; i < n; ++i)
366                 (z[i])++;
367 
368         return (1);
369 }
370 

~ [ 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.