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

Parent Directory | 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