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

Linux Cross Reference
Tina5/tina-libs/tina/math/mathMatv_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/mathMatv_eigen.c,v $
 23  * Date    :  $Date: 2003/10/06 12:29:47 $
 24  * Version :  $Revision: 1.5 $
 25  * CVS Id  :  $Id: mathMatv_eigen.c,v 1.5 2003/10/06 12:29:47 neil Exp $
 26  *
 27  * Author  :  Legacy TINA
 28  *
 29  * Notes :
 30  *
 31  *********
 32 */
 33 /** 
 34  *  @file
 35  *  @brief  Code to find eigenvectors and eigenvalues of a symmetric n by n matrix.     
 36  *
 37  *  Implementation of methods as described in Numerical Recipes, (see chapter 11 (3rd ed), 
 38  *  esp 11.2 and 11.3) and originally distributed in EISPACK.
 39  * 
 40 */
 41 
 42 #include "mathMatv_eigen.h"
 43 
 44 #if HAVE_CONFIG_H
 45 #include <config.h>
 46 #endif
 47 
 48 #include <math.h>
 49 #include <tina/sys/sysDef.h>
 50 #include <tina/sys/sysPro.h>
 51 #include <tina/math/math_MatvDef.h>
 52 
 53 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
 54 
 55 static void     tred2(int n, double **a, double *d, double *e)
 56 {
 57         int             i, j, k, l;
 58         double          scale, hh, h, g, f;
 59 
 60         /*
 61         start arrays at 0 not 1
 62         */
 63         for (i = 0; i < n; ++i)
 64                 (a[i])--;
 65         a--;
 66         d--;
 67         e--;
 68 
 69         for (i = n; i >= 2; i--)
 70         {
 71                 l = i - 1;
 72                 h = scale = 0.0;
 73                 if (l > 1)
 74                 {
 75                         for (k = 1; k <= l; k++)
 76                                 scale += fabs(a[i][k]);
 77                         if (scale == 0.0)
 78                                 e[i] = a[i][l];
 79                         else
 80                         {
 81                                 for (k = 1; k <= l; k++)
 82                                 {
 83                                         a[i][k] /= scale;
 84                                         h += a[i][k] * a[i][k];
 85                                 }
 86                                 f = a[i][l];
 87                                 g = f > 0 ? -sqrt(h) : sqrt(h);
 88                                 e[i] = scale * g;
 89                                 h -= f * g;
 90                                 a[i][l] = f - g;
 91                                 f = 0.0;
 92                                 for (j = 1; j <= l; j++)
 93                                 {
 94                                         a[j][i] = a[i][j] / h;
 95                                         g = 0.0;
 96                                         for (k = 1; k <= j; k++)
 97                                                 g += a[j][k] * a[i][k];
 98                                         for (k = j + 1; k <= l; k++)
 99                                                 g += a[k][j] * a[i][k];
100                                         e[j] = g / h;
101                                         f += e[j] * a[i][j];
102                                 }
103                                 hh = f / (h + h);
104                                 for (j = 1; j <= l; j++)
105                                 {
106                                         f = a[i][j];
107                                         e[j] = g = e[j] - hh * f;
108                                         for (k = 1; k <= j; k++)
109                                                 a[j][k] -= (f * e[k] + g * a[i][k]);
110                                 }
111                         }
112                 } else
113                         e[i] = a[i][l];
114                 d[i] = h;
115         }
116         d[1] = 0.0;
117         e[1] = 0.0;
118         for (i = 1; i <= n; i++)
119         {
120                 l = i - 1;
121                 if (d[i])
122                 {
123                         for (j = 1; j <= l; j++)
124                         {
125                                 g = 0.0;
126                                 for (k = 1; k <= l; k++)
127                                         g += a[i][k] * a[k][j];
128                                 for (k = 1; k <= l; k++)
129                                         a[k][j] -= g * a[k][i];
130                         }
131                 }
132                 d[i] = a[i][i];
133                 a[i][i] = 1.0;
134                 for (j = 1; j <= l; j++)
135                         a[j][i] = a[i][j] = 0.0;
136         }
137 
138         a++;
139         d++;
140         e++;
141         for (i = 0; i < n; ++i)
142                 (a[i])++;
143 }
144 
145 static int      tqli(int n, double **z, double *d, double *e)
146 {
147         int             m, l, iter, i, k;
148         double          s, r, p, g, f, dd, c, b;
149 
150         /*
151         start arrays at 0 not 1
152         */
153         for (i = 0; i < n; ++i)
154                 (z[i])--;
155         z--;
156         d--;
157         e--;
158 
159         for (i = 2; i <= n; i++)
160                 e[i - 1] = e[i];
161         e[n] = 0.0;
162         for (l = 1; l <= n; l++)
163         {
164                 iter = 0;
165                 do
166                 {
167                         for (m = l; m <= n - 1; m++)
168                         {
169                                 dd = fabs(d[m]) + fabs(d[m + 1]);
170                                 if (fabs(e[m]) + dd == dd)
171                                         break;
172                         }
173                         if (m != l)
174                         {
175                                 if (iter++ == 30)
176                                         return (0);
177                                 g = (d[l + 1] - d[l]) / (2.0 * e[l]);
178                                 r = sqrt((g * g) + 1.0);
179                                 g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
180                                 s = c = 1.0;
181                                 p = 0.0;
182                                 for (i = m - 1; i >= l; i--)
183                                 {
184                                         f = s * e[i];
185                                         b = c * e[i];
186                                         if (fabs(f) >= fabs(g))
187                                         {
188                                                 c = g / f;
189                                                 r = sqrt((c * c) + 1.0);
190                                                 e[i + 1] = f * r;
191                                                 c *= (s = 1.0 / r);
192                                         } else
193                                         {
194                                                 s = f / g;
195                                                 r = sqrt((s * s) + 1.0);
196                                                 e[i + 1] = g * r;
197                                                 s *= (c = 1.0 / r);
198                                         }
199                                         g = d[i + 1] - p;
200                                         r = (d[i] - g) * s + 2.0 * c * b;
201                                         p = s * r;
202                                         d[i + 1] = g + p;
203                                         g = c * r - b;
204                                         for (k = 1; k <= n; k++)
205                                         {
206                                                 f = z[k][i + 1];
207                                                 z[k][i + 1] = s * z[k][i] + c * f;
208                                                 z[k][i] = c * z[k][i] - s * f;
209                                         }
210                                 }
211                                 d[l] = d[l] - p;
212                                 e[l] = g;
213                                 e[m] = 0.0;
214                         }
215                 }
216                 while (m != l);
217         }
218 
219         z++;
220         d++;
221         e++;
222         for (i = 0; i < n; ++i)
223                 (z[i])++;
224 
225         return (1);
226 }
227 
228 /*
229 Find eigenvalues and eigenvectors of symmetric  n  by n matrix a,
230 */
231 static int      eigen(int n, double **a, double **evec, double *eval)
232 {
233         int             i, j, k, success;
234         double         *d = tvector_alloc(0, n, double);
235         double         *e = tvector_alloc(0, n, double);
236 
237         for (i = 0; i < n; i++)
238                 for (j = 0; j < n; j++)
239                         evec[i][j] = a[i][j];
240 
241         tred2(n, evec, eval, e);
242         success = tqli(n, evec, eval, e);
243 
244         if (!success)
245         {
246                 tvector_free(d, 0, double);
247                 tvector_free(e, 0, double);
248                 return (success);
249         }
250         for (j = 0; j < n - 1; j++)
251         {
252                 double          emax = eval[j];
253                 int             jmax = j;
254                 for (k = j + 1; k < n; k++)
255                         if (emax < eval[k])
256                         {
257                                 emax = eval[k];
258                                 jmax = k;
259                         }
260                 if (jmax != j)
261                 {
262                         SWAP(double, eval[j], eval[jmax]);
263                         for (i = 0; i < n; i++)
264                                 SWAP(double, evec[i][j], evec[i][jmax]);
265                 }
266         }
267         tvector_free(d, 0, double);
268         tvector_free(e, 0, double);
269         return (success);
270 }
271 
272 /*
273 eval's in descending order
274 evecs in columns of matrix evec
275 */
276 int             mat_eigen(Mat * a, Mat * evec, Vec * eval)
277 {
278         return (eigen(a->n, a->el, evec->el, eval->el));
279 }
280 
281 
282 #undef SIGN
283 

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