Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Annotation of /trunk/src/lib/c-target/eigen.c
ViewVC logotype

Annotation of /trunk/src/lib/c-target/eigen.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1115 - (view) (download) (as text)

1 : jhr 1115 /*
2 :     Teem: Tools to process and visualize scientific data and images
3 :     Copyright (C) 2011, 2010, 2009, University of Chicago
4 :     Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
5 :     Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
6 :    
7 :     This library is free software; you can redistribute it and/or
8 :     modify it under the terms of the GNU Lesser General Public License
9 :     (LGPL) as published by the Free Software Foundation; either
10 :     version 2.1 of the License, or (at your option) any later version.
11 :     The terms of redistributing and/or modifying this software also
12 :     include exceptions to the LGPL that facilitate static linking.
13 :    
14 :     This library is distributed in the hope that it will be useful,
15 :     but WITHOUT ANY WARRANTY; without even the implied warranty of
16 :     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 :     Lesser General Public License for more details.
18 :    
19 :     You should have received a copy of the GNU Lesser General Public License
20 :     along with this library; if not, write to Free Software Foundation, Inc.,
21 :     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 :     */
23 :    
24 :    
25 :     #include "../ten.h"
26 :    
27 :     char *info = ("tests tenEigensolve_d and new stand-alone function.");
28 :    
29 :     #define ROOT_TRIPLE 2 /* ell_cubic_root_triple */
30 :     #define ROOT_SINGLE_DOUBLE 3 /* ell_cubic_root_single_double */
31 :     #define ROOT_THREE 4 /* ell_cubic_root_three */
32 :    
33 :     #define ABS(a) (((a) > 0.0f ? (a) : -(a)))
34 :     #define VEC_SET(v, a, b, c) \
35 :     ((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))
36 :     #define VEC_DOT(v1, v2) \
37 :     ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
38 :     #define VEC_CROSS(v3, v1, v2) \
39 :     ((v3)[0] = (v1)[1]*(v2)[2] - (v1)[2]*(v2)[1], \
40 :     (v3)[1] = (v1)[2]*(v2)[0] - (v1)[0]*(v2)[2], \
41 :     (v3)[2] = (v1)[0]*(v2)[1] - (v1)[1]*(v2)[0])
42 :     #define VEC_ADD(v1, v2) \
43 :     ((v1)[0] += (v2)[0], \
44 :     (v1)[1] += (v2)[1], \
45 :     (v1)[2] += (v2)[2])
46 :     #define VEC_SUB(v1, v2) \
47 :     ((v1)[0] -= (v2)[0], \
48 :     (v1)[1] -= (v2)[1], \
49 :     (v1)[2] -= (v2)[2])
50 :     #define VEC_SCL(v1, s) \
51 :     ((v1)[0] *= (s), \
52 :     (v1)[1] *= (s), \
53 :     (v1)[2] *= (s))
54 :     #define VEC_LEN(v) (sqrt(VEC_DOT(v,v)))
55 :     #define VEC_NORM(v, len) ((len) = VEC_LEN(v), VEC_SCL(v, 1.0/len))
56 :     #define VEC_SCL_SUB(v1, s, v2) \
57 :     ((v1)[0] -= (s)*(v2)[0], \
58 :     (v1)[1] -= (s)*(v2)[1], \
59 :     (v1)[2] -= (s)*(v2)[2])
60 :     #define VEC_COPY(v1, v2) \
61 :     ((v1)[0] = (v2)[0], \
62 :     (v1)[1] = (v2)[1], \
63 :     (v1)[2] = (v2)[2])
64 :    
65 :     /*
66 :     ** All the three given vectors span only a 2D space, and this finds
67 :     ** the normal to that plane. Simply sums up all the pair-wise
68 :     ** cross-products to get a good estimate. Trick is getting the cross
69 :     ** products to line up before summing.
70 :     */
71 :     void
72 :     nullspace1(double ret[3],
73 :     const double r0[3], const double r1[3], const double r2[3]) {
74 :     double crs[3];
75 :    
76 :     /* ret = r0 x r1 */
77 :     VEC_CROSS(ret, r0, r1);
78 :     /* crs = r1 x r2 */
79 :     VEC_CROSS(crs, r1, r2);
80 :     /* ret += crs or ret -= crs; whichever makes ret longer */
81 :     if (VEC_DOT(ret, crs) > 0) {
82 :     VEC_ADD(ret, crs);
83 :     } else {
84 :     VEC_SUB(ret, crs);
85 :     }
86 :     /* crs = r0 x r2 */
87 :     VEC_CROSS(crs, r0, r2);
88 :     /* ret += crs or ret -= crs; whichever makes ret longer */
89 :     if (VEC_DOT(ret, crs) > 0) {
90 :     VEC_ADD(ret, crs);
91 :     } else {
92 :     VEC_SUB(ret, crs);
93 :     }
94 :    
95 :     return;
96 :     }
97 :    
98 :     /*
99 :     ** All vectors are in the same 1D space, we have to find two
100 :     ** mutually vectors perpendicular to that span
101 :     */
102 :     void
103 :     nullspace2(double reta[3], double retb[3],
104 :     const double r0[3], const double r1[3], const double r2[3]) {
105 :     double sqr[3], sum[3];
106 :     int idx;
107 :    
108 :     VEC_COPY(sum, r0);
109 :     if (VEC_DOT(sum, r1) > 0) {
110 :     VEC_ADD(sum, r1);
111 :     } else {
112 :     VEC_SUB(sum, r1);
113 :     }
114 :     if (VEC_DOT(sum, r2) > 0) {
115 :     VEC_ADD(sum, r2);
116 :     } else {
117 :     VEC_SUB(sum, r2);
118 :     }
119 :     /* find largest component, to get most stable expression for a
120 :     perpendicular vector */
121 :     sqr[0] = sum[0]*sum[0];
122 :     sqr[1] = sum[1]*sum[1];
123 :     sqr[2] = sum[2]*sum[2];
124 :     idx = 0;
125 :     if (sqr[0] < sqr[1])
126 :     idx = 1;
127 :     if (sqr[idx] < sqr[2])
128 :     idx = 2;
129 :     /* reta will be perpendicular to sum */
130 :     if (0 == idx) {
131 :     VEC_SET(reta, sum[1] - sum[2], -sum[0], sum[0]);
132 :     } else if (1 == idx) {
133 :     VEC_SET(reta, -sum[1], sum[0] - sum[2], sum[1]);
134 :     } else {
135 :     VEC_SET(reta, -sum[2], sum[2], sum[0] - sum[1]);
136 :     }
137 :     /* and now retb will be perpendicular to both reta and sum */
138 :     VEC_CROSS(retb, reta, sum);
139 :     return;
140 :     }
141 :    
142 :     /*
143 :     ** Eigensolver for symmetric 3x3 matrix:
144 :     **
145 :     ** M00 M01 M02
146 :     ** M01 M11 M12
147 :     ** M02 M12 M22
148 :     **
149 :     ** Must be passed eval[3] vector, and will compute eigenvalues
150 :     ** only if evec[9] is non-NULL. Computed eigenvectors are at evec+0,
151 :     ** evec+3, and evec+6.
152 :     **
153 :     ** Return value indicates something about the eigenvalue solution to
154 :     ** the cubic characteristic equation; see ROOT_ #defines above
155 :     **
156 :     ** Relies on the ABS and VEC_* macros above, as well as math functions
157 :     ** atan2(), sin(), cos(), sqrt(), and airCbrt(), defined as:
158 :    
159 :     double
160 :     airCbrt(double v) {
161 :     #if defined(_WIN32) || defined(__STRICT_ANSI__)
162 :     return (v < 0.0 ? -pow(-v,1.0/3.0) : pow(v,1.0/3.0));
163 :     #else
164 :     return cbrt(v);
165 :     #endif
166 :     }
167 :    
168 :     **
169 :     ** HEY: the numerical precision issues here are very subtle, and
170 :     ** merit some more scrutiny. With evals (1.000001, 1, 1), for example,
171 :     ** whether it comes back as a single/double root, vs three distinct roots,
172 :     ** is determines by the comparison between "D" and "epsilon", and the
173 :     ** setting of epsilon seems pretty arbitrary at this point...
174 :     **
175 :     */
176 :     int
177 :     evals(double eval[3],
178 :     const double _M00, const double _M01, const double _M02,
179 :     const double _M11, const double _M12,
180 :     const double _M22) {
181 :    
182 :     #include "teigen-evals-A.c"
183 :    
184 :     #include "teigen-evals-B.c"
185 :    
186 :     return roots;
187 :     }
188 :    
189 :     int
190 :     evals_evecs(double eval[3], double evec[9],
191 :     const double _M00, const double _M01, const double _M02,
192 :     const double _M11, const double _M12,
193 :     const double _M22) {
194 :     double r0[3], r1[3], r2[3], crs[3], len, dot;
195 :    
196 :     #include "teigen-evals-A.c"
197 :    
198 :     /* r0, r1, r2 are the vectors we manipulate to
199 :     find the nullspaces of M - lambda*I */
200 :     VEC_SET(r0, 0.0, M01, M02);
201 :     VEC_SET(r1, M01, 0.0, M12);
202 :     VEC_SET(r2, M02, M12, 0.0);
203 :     if (ROOT_THREE == roots) {
204 :     r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
205 :     nullspace1(evec+0, r0, r1, r2);
206 :     r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
207 :     nullspace1(evec+3, r0, r1, r2);
208 :     r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
209 :     nullspace1(evec+6, r0, r1, r2);
210 :     } else if (ROOT_SINGLE_DOUBLE == roots) {
211 :     if (eval[1] == eval[2]) {
212 :     /* one big (eval[0]) , two small (eval[1,2]) */
213 :     r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
214 :     nullspace1(evec+0, r0, r1, r2);
215 :     r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
216 :     nullspace2(evec+3, evec+6, r0, r1, r2);
217 :     }
218 :     else {
219 :     /* two big (eval[0,1]), one small (eval[2]) */
220 :     r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
221 :     nullspace2(evec+0, evec+3, r0, r1, r2);
222 :     r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
223 :     nullspace1(evec+6, r0, r1, r2);
224 :     }
225 :     } else {
226 :     /* ROOT_TRIPLE == roots; use any basis for eigenvectors */
227 :     VEC_SET(evec+0, 1, 0, 0);
228 :     VEC_SET(evec+3, 0, 1, 0);
229 :     VEC_SET(evec+6, 0, 0, 1);
230 :     }
231 :     /* we always make sure its really orthonormal; keeping fixed the
232 :     eigenvector associated with the largest-magnitude eigenvalue */
233 :     if (ABS(eval[0]) > ABS(eval[2])) {
234 :     /* normalize evec+0 but don't move it */
235 :     VEC_NORM(evec+0, len);
236 :     dot = VEC_DOT(evec+0, evec+3); VEC_SCL_SUB(evec+3, dot, evec+0);
237 :     VEC_NORM(evec+3, len);
238 :     dot = VEC_DOT(evec+0, evec+6); VEC_SCL_SUB(evec+6, dot, evec+0);
239 :     dot = VEC_DOT(evec+3, evec+6); VEC_SCL_SUB(evec+6, dot, evec+3);
240 :     VEC_NORM(evec+6, len);
241 :     } else {
242 :     /* normalize evec+6 but don't move it */
243 :     VEC_NORM(evec+6, len);
244 :     dot = VEC_DOT(evec+6, evec+3); VEC_SCL_SUB(evec+3, dot, evec+6);
245 :     VEC_NORM(evec+3, len);
246 :     dot = VEC_DOT(evec+3, evec+0); VEC_SCL_SUB(evec+0, dot, evec+3);
247 :     dot = VEC_DOT(evec+6, evec+0); VEC_SCL_SUB(evec+0, dot, evec+6);
248 :     VEC_NORM(evec+0, len);
249 :     }
250 :    
251 :     /* to be nice, make it right-handed */
252 :     VEC_CROSS(crs, evec+0, evec+3);
253 :     if (0 > VEC_DOT(crs, evec+6)) {
254 :     VEC_SCL(evec+6, -1);
255 :     }
256 :    
257 :     #include "teigen-evals-B.c"
258 :    
259 :     return roots;
260 :     }
261 :    
262 :    
263 :     void
264 :     testeigen(double tt[7], double eval[3], double evec[9]) {
265 :     double mat[9], dot[3], cross[3];
266 :     unsigned int ii;
267 :    
268 :     TEN_T2M(mat, tt);
269 :     printf("evals %g %g %g\n", eval[0], eval[1], eval[2]);
270 :     printf("evec0 (%g) %g %g %g\n",
271 :     ELL_3V_LEN(evec + 0), evec[0], evec[1], evec[2]);
272 :     printf("evec1 (%g) %g %g %g\n",
273 :     ELL_3V_LEN(evec + 3), evec[3], evec[4], evec[5]);
274 :     printf("evec2 (%g) %g %g %g\n",
275 :     ELL_3V_LEN(evec + 6), evec[6], evec[7], evec[8]);
276 :     printf("Mv - lv: (len) X Y Z (should be ~zeros)\n");
277 :     for (ii=0; ii<3; ii++) {
278 :     double uu[3], vv[3], dd[3];
279 :     ELL_3MV_MUL(uu, mat, evec + 3*ii);
280 :     ELL_3V_SCALE(vv, eval[ii], evec + 3*ii);
281 :     ELL_3V_SUB(dd, uu, vv);
282 :     printf("%d: (%g) %g %g %g\n", ii, ELL_3V_LEN(dd), dd[0], dd[1], dd[2]);
283 :     }
284 :     dot[0] = ELL_3V_DOT(evec + 0, evec + 3);
285 :     dot[1] = ELL_3V_DOT(evec + 0, evec + 6);
286 :     dot[2] = ELL_3V_DOT(evec + 3, evec + 6);
287 :     printf("pairwise dots: (%g) %g %g %g\n",
288 :     ELL_3V_LEN(dot), dot[0], dot[1], dot[2]);
289 :     ELL_3V_CROSS(cross, evec+0, evec+3);
290 :     printf("right-handed: %g\n", ELL_3V_DOT(evec+6, cross));
291 :     return;
292 :     }
293 :    
294 :     int
295 :     main(int argc, char *argv[]) {
296 :     char *me;
297 :     hestOpt *hopt=NULL;
298 :     airArray *mop;
299 :    
300 :     double _tt[6], tt[7], ss, pp[3], qq[4], rot[9], mat1[9], mat2[9], tmp,
301 :     evalA[3], evecA[9], evalB[3], evecB[9];
302 :     int roots;
303 :    
304 :     mop = airMopNew();
305 :     me = argv[0];
306 :     hestOptAdd(&hopt, NULL, "m00 m01 m02 m11 m12 m22",
307 :     airTypeDouble, 6, 6, _tt, NULL, "symmtric matrix coeffs");
308 :     hestOptAdd(&hopt, "p", "vec", airTypeDouble, 3, 3, pp, "0 0 0",
309 :     "rotation as P vector");
310 :     hestOptAdd(&hopt, "s", "scl", airTypeDouble, 1, 1, &ss, "1.0",
311 :     "scaling");
312 :     hestParseOrDie(hopt, argc-1, argv+1, NULL,
313 :     me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);
314 :     airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);
315 :     airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);
316 :    
317 :     ELL_6V_COPY(tt + 1, _tt);
318 :     tt[0] = 1.0;
319 :     TEN_T_SCALE(tt, ss, tt);
320 :    
321 :     ELL_4V_SET(qq, 1, pp[0], pp[1], pp[2]);
322 :     ELL_4V_NORM(qq, qq, tmp);
323 :     ell_q_to_3m_d(rot, qq);
324 :     printf("%s: rot\n", me);
325 :     printf(" %g %g %g\n", rot[0], rot[1], rot[2]);
326 :     printf(" %g %g %g\n", rot[3], rot[4], rot[5]);
327 :     printf(" %g %g %g\n", rot[6], rot[7], rot[8]);
328 :    
329 :     TEN_T2M(mat1, tt);
330 :     ell_3m_mul_d(mat2, rot, mat1);
331 :     ELL_3M_TRANSPOSE_IP(rot, tmp);
332 :     ell_3m_mul_d(mat1, mat2, rot);
333 :     TEN_M2T(tt, mat1);
334 :    
335 :     printf("input matrix = \n %g %g %g\n %g %g\n %g\n",
336 :     tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);
337 :    
338 :     printf("================== tenEigensolve_d ==================\n");
339 :     roots = tenEigensolve_d(evalA, evecA, tt);
340 :     printf("%s roots\n", airEnumStr(ell_cubic_root, roots));
341 :     testeigen(tt, evalA, evecA);
342 :    
343 :     printf("================== new eigensolve ==================\n");
344 :     roots = evals(evalB, tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);
345 :     printf("%s roots: %g %g %g\n", airEnumStr(ell_cubic_root, roots),
346 :     evalB[0], evalB[1], evalB[2]);
347 :     roots = evals_evecs(evalB, evecB,
348 :     tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);
349 :     printf("%s roots\n", airEnumStr(ell_cubic_root, roots));
350 :     testeigen(tt, evalB, evecB);
351 :    
352 :     airMopOkay(mop);
353 :     return 0;
354 :     }

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0