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

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

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