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

Linux Cross Reference
Tina4/src/math/matvec/mat_eigen.c

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

  1 #include <math.h>
  2 #include <tina/sys.h>
  3 #include <tina/sysfuncs.h>
  4 #include <tina/math.h>
  5 #include <tina/mathfuncs.h>
  6 
  7 #define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
  8 
  9 static void tred2(int n, double **a, double *d, double *e)
 10 {
 11     int l, k, j, i;
 12     double scale, hh, h, g, f;
 13 
 14     /*
 15     start arrays at 0 not 1
 16     */
 17     for (i = 0; i < n; ++i)
 18         (a[i])--;
 19     a--;
 20     d--;
 21     e--;
 22 
 23     for (i = n; i >= 2; i--)
 24     {
 25         l = i - 1;
 26         h = scale = 0.0;
 27         if (l > 1)
 28         {
 29             for (k = 1; k <= l; k++)
 30                 scale += fabs(a[i][k]);
 31             if (scale == 0.0)
 32                 e[i] = a[i][l];
 33             else
 34             {
 35                 for (k = 1; k <= l; k++)
 36                 {
 37                     a[i][k] /= scale;
 38                     h += a[i][k] * a[i][k];
 39                 }
 40                 f = a[i][l];
 41                 g = f > 0 ? -sqrt(h) : sqrt(h);
 42                 e[i] = scale * g;
 43                 h -= f * g;
 44                 a[i][l] = f - g;
 45                 f = 0.0;
 46                 for (j = 1; j <= l; j++)
 47                 {
 48                     a[j][i] = a[i][j] / h;
 49                     g = 0.0;
 50                     for (k = 1; k <= j; k++)
 51                         g += a[j][k] * a[i][k];
 52                     for (k = j + 1; k <= l; k++)
 53                         g += a[k][j] * a[i][k];
 54                     e[j] = g / h;
 55                     f += e[j] * a[i][j];
 56                 }
 57                 hh = f / (h + h);
 58                 for (j = 1; j <= l; j++)
 59                 {
 60                     f = a[i][j];
 61                     e[j] = g = e[j] - hh * f;
 62                     for (k = 1; k <= j; k++)
 63                         a[j][k] -= (f * e[k] + g * a[i][k]);
 64                 }
 65             }
 66         }
 67         else
 68             e[i] = a[i][l];
 69         d[i] = h;
 70     }
 71     d[1] = 0.0;
 72     e[1] = 0.0;
 73     for (i = 1; i <= n; i++)
 74     {
 75         l = i - 1;
 76         if (d[i])
 77         {
 78             for (j = 1; j <= l; j++)
 79             {
 80                 g = 0.0;
 81                 for (k = 1; k <= l; k++)
 82                     g += a[i][k] * a[k][j];
 83                 for (k = 1; k <= l; k++)
 84                     a[k][j] -= g * a[k][i];
 85             }
 86         }
 87         d[i] = a[i][i];
 88         a[i][i] = 1.0;
 89         for (j = 1; j <= l; j++)
 90             a[j][i] = a[i][j] = 0.0;
 91     }
 92 
 93     a++;
 94     d++;
 95     e++;
 96     for (i = 0; i < n; ++i)
 97         (a[i])++;
 98 }
 99 
100 static int tqli(int n, double **z, double *d, double *e)
101 {
102     int m, l, iter, i, k;
103     double s, r, p, g, f, dd, c, b;
104 
105     /*
106     start arrays at 0 not 1
107     */
108     for (i = 0; i < n; ++i)
109         (z[i])--;
110     z--;
111     d--;
112     e--;
113 
114     for (i = 2; i <= n; i++)
115         e[i - 1] = e[i];
116     e[n] = 0.0;
117     for (l = 1; l <= n; l++)
118     {
119         iter = 0;
120         do
121         {
122             for (m = l; m <= n - 1; m++)
123             {
124                 dd = fabs(d[m]) + fabs(d[m + 1]);
125                 if (fabs(e[m]) + dd == dd)
126                     break;
127             }
128             if (m != l)
129             {
130                 if (iter++ == 30)
131                     return (0);
132                 g = (d[l + 1] - d[l]) / (2.0 * e[l]);
133                 r = sqrt((g * g) + 1.0);
134                 g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
135                 s = c = 1.0;
136                 p = 0.0;
137                 for (i = m - 1; i >= l; i--)
138                 {
139                     f = s * e[i];
140                     b = c * e[i];
141                     if (fabs(f) >= fabs(g))
142                     {
143                         c = g / f;
144                         r = sqrt((c * c) + 1.0);
145                         e[i + 1] = f * r;
146                         c *= (s = 1.0 / r);
147                     }
148                     else
149                     {
150                         s = f / g;
151                         r = sqrt((s * s) + 1.0);
152                         e[i + 1] = g * r;
153                         s *= (c = 1.0 / r);
154                     }
155                     g = d[i + 1] - p;
156                     r = (d[i] - g) * s + 2.0 * c * b;
157                     p = s * r;
158                     d[i + 1] = g + p;
159                     g = c * r - b;
160                     for (k = 1; k <= n; k++)
161                     {
162                         f = z[k][i + 1];
163                         z[k][i + 1] = s * z[k][i] + c * f;
164                         z[k][i] = c * z[k][i] - s * f;
165                     }
166                 }
167                 d[l] = d[l] - p;
168                 e[l] = g;
169                 e[m] = 0.0;
170             }
171         }
172         while (m != l);
173     }
174 
175     z++;
176     d++;
177     e++;
178     for (i = 0; i < n; ++i)
179         (z[i])++;
180 
181     return (1);
182 }
183 
184 /*
185 Find eigenvalues and eigenvectors of symmetric  n  by n matrix a,
186 */
187 static int eigen(int n, double **a, double **evec, double *eval)
188 {
189     int i, j, k, success;
190     double *d = tvector_alloc(0, n, double);
191     double *e = tvector_alloc(0, n, double);
192 
193     for(i = 0; i < n; i++)
194         for(j = 0; j < n; j++)
195             evec[i][j] = a[i][j];
196 
197     tred2(n, evec, eval, e);
198     success = tqli(n, evec, eval, e);
199 
200     if(!success)
201     {
202         tvector_free(d, 0, double);
203         tvector_free(e, 0, double);
204         return(success);
205     }
206 
207     for (j = 0; j < n - 1; j++)
208     {
209         double emax = eval[j];
210         int jmax = j;
211         for (k = j + 1; k < n; k++)
212             if (emax < eval[k])
213             {
214                 emax = eval[k];
215                 jmax = k;
216             }
217         if (jmax != j)
218         {
219             SWAP(double, eval[j], eval[jmax]);
220             for(i = 0; i < n; i++)
221                 SWAP(double, evec[i][j], evec[i][jmax]);
222         }
223     }
224     tvector_free(d, 0, double);
225     tvector_free(e, 0, double);
226     return (success);
227 }
228 
229 /*
230 eval's in descending order
231 evecs in columns of matrix evec
232 */
233 int mat_eigen(Mat *a, Mat *evec, Vec *eval)
234 {
235     return(eigen(a->n, a->el, evec->el, eval->el));
236 }
237 
238 
239 #undef SIGN
240 

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