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

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

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