Annotation of /branches/pure-cfg/src/lib/c-target/eigen.c

Parent Directory | Revision Log

---

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

1 :	jhr	1172	/*! \file eigen.c
2 :			*
3 :			* \author Gordon Kindlmann
4 :			*/
5 :
6 :	jhr	1106	/*
7 :	jhr	1172	* COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8 :			* All rights reserved.
9 :			*/
10 :
11 :			/*
12 :	jhr	1106	Teem: Tools to process and visualize scientific data and images
13 :			Copyright (C) 2011, 2010, 2009, University of Chicago
14 :			Copyright (C) 2008, 2007, 2006, 2005 Gordon Kindlmann
15 :			Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998 University of Utah
16 :
17 :			This library is free software; you can redistribute it and/or
18 :			modify it under the terms of the GNU Lesser General Public License
19 :			(LGPL) as published by the Free Software Foundation; either
20 :			version 2.1 of the License, or (at your option) any later version.
21 :			The terms of redistributing and/or modifying this software also
22 :			include exceptions to the LGPL that facilitate static linking.
23 :
24 :			This library is distributed in the hope that it will be useful,
25 :			but WITHOUT ANY WARRANTY; without even the implied warranty of
26 :			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
27 :			Lesser General Public License for more details.
28 :
29 :			You should have received a copy of the GNU Lesser General Public License
30 :			along with this library; if not, write to Free Software Foundation, Inc.,
31 :			51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
32 :			*/
33 :
34 :	jhr	1172	#ifdef PORTED
35 :	jhr	1106
36 :	jhr	1172	#include <Diderot/diderot.h>
37 :	jhr	1106
38 :			#define ROOT_TRIPLE 2 /* ell_cubic_root_triple */
39 :			#define ROOT_SINGLE_DOUBLE 3 /* ell_cubic_root_single_double */
40 :			#define ROOT_THREE 4 /* ell_cubic_root_three */
41 :
42 :			#define ABS(a) (((a) > 0.0f ? (a) : -(a)))
43 :			#define VEC_SET(v, a, b, c) \
44 :			((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))
45 :			#define VEC_DOT(v1, v2) \
46 :			((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
47 :			#define VEC_CROSS(v3, v1, v2) \
48 :			((v3)[0] = (v1)[1]*(v2)[2] - (v1)[2]*(v2)[1], \
49 :			(v3)[1] = (v1)[2]*(v2)[0] - (v1)[0]*(v2)[2], \
50 :			(v3)[2] = (v1)[0]*(v2)[1] - (v1)[1]*(v2)[0])
51 :			#define VEC_ADD(v1, v2) \
52 :			((v1)[0] += (v2)[0], \
53 :			(v1)[1] += (v2)[1], \
54 :			(v1)[2] += (v2)[2])
55 :			#define VEC_SUB(v1, v2) \
56 :			((v1)[0] -= (v2)[0], \
57 :			(v1)[1] -= (v2)[1], \
58 :			(v1)[2] -= (v2)[2])
59 :			#define VEC_SCL(v1, s) \
60 :			((v1)[0] *= (s), \
61 :			(v1)[1] *= (s), \
62 :			(v1)[2] *= (s))
63 :			#define VEC_LEN(v) (sqrt(VEC_DOT(v,v)))
64 :			#define VEC_NORM(v, len) ((len) = VEC_LEN(v), VEC_SCL(v, 1.0/len))
65 :			#define VEC_SCL_SUB(v1, s, v2) \
66 :			((v1)[0] -= (s)*(v2)[0], \
67 :			(v1)[1] -= (s)*(v2)[1], \
68 :			(v1)[2] -= (s)*(v2)[2])
69 :			#define VEC_COPY(v1, v2) \
70 :			((v1)[0] = (v2)[0], \
71 :			(v1)[1] = (v2)[1], \
72 :			(v1)[2] = (v2)[2])
73 :
74 :			/*
75 :			** All the three given vectors span only a 2D space, and this finds
76 :			** the normal to that plane. Simply sums up all the pair-wise
77 :			** cross-products to get a good estimate. Trick is getting the cross
78 :			** products to line up before summing.
79 :			*/
80 :			void
81 :			nullspace1(double ret[3],
82 :			const double r0[3], const double r1[3], const double r2[3]) {
83 :			double crs[3];
84 :
85 :			/* ret = r0 x r1 */
86 :			VEC_CROSS(ret, r0, r1);
87 :			/* crs = r1 x r2 */
88 :			VEC_CROSS(crs, r1, r2);
89 :			/* ret += crs or ret -= crs; whichever makes ret longer */
90 :			if (VEC_DOT(ret, crs) > 0) {
91 :			VEC_ADD(ret, crs);
92 :			} else {
93 :			VEC_SUB(ret, crs);
94 :			}
95 :			/* crs = r0 x r2 */
96 :			VEC_CROSS(crs, r0, r2);
97 :			/* ret += crs or ret -= crs; whichever makes ret longer */
98 :			if (VEC_DOT(ret, crs) > 0) {
99 :			VEC_ADD(ret, crs);
100 :			} else {
101 :			VEC_SUB(ret, crs);
102 :			}
103 :
104 :			return;
105 :			}
106 :
107 :			/*
108 :			** All vectors are in the same 1D space, we have to find two
109 :			** mutually vectors perpendicular to that span
110 :			*/
111 :			void
112 :			nullspace2(double reta[3], double retb[3],
113 :			const double r0[3], const double r1[3], const double r2[3]) {
114 :			double sqr[3], sum[3];
115 :			int idx;
116 :
117 :			VEC_COPY(sum, r0);
118 :			if (VEC_DOT(sum, r1) > 0) {
119 :			VEC_ADD(sum, r1);
120 :			} else {
121 :			VEC_SUB(sum, r1);
122 :			}
123 :			if (VEC_DOT(sum, r2) > 0) {
124 :			VEC_ADD(sum, r2);
125 :			} else {
126 :			VEC_SUB(sum, r2);
127 :			}
128 :			/* find largest component, to get most stable expression for a
129 :			perpendicular vector */
130 :			sqr[0] = sum[0]*sum[0];
131 :			sqr[1] = sum[1]*sum[1];
132 :			sqr[2] = sum[2]*sum[2];
133 :			idx = 0;
134 :			if (sqr[0] < sqr[1])
135 :			idx = 1;
136 :			if (sqr[idx] < sqr[2])
137 :			idx = 2;
138 :			/* reta will be perpendicular to sum */
139 :			if (0 == idx) {
140 :			VEC_SET(reta, sum[1] - sum[2], -sum[0], sum[0]);
141 :			} else if (1 == idx) {
142 :			VEC_SET(reta, -sum[1], sum[0] - sum[2], sum[1]);
143 :			} else {
144 :			VEC_SET(reta, -sum[2], sum[2], sum[0] - sum[1]);
145 :			}
146 :			/* and now retb will be perpendicular to both reta and sum */
147 :			VEC_CROSS(retb, reta, sum);
148 :			return;
149 :			}
150 :
151 :			/*
152 :			** Eigensolver for symmetric 3x3 matrix:
153 :			**
154 :			** M00 M01 M02
155 :			** M01 M11 M12
156 :			** M02 M12 M22
157 :			**
158 :			** Must be passed eval[3] vector, and will compute eigenvalues
159 :			** only if evec[9] is non-NULL. Computed eigenvectors are at evec+0,
160 :			** evec+3, and evec+6.
161 :			**
162 :			** Return value indicates something about the eigenvalue solution to
163 :			** the cubic characteristic equation; see ROOT_ #defines above
164 :			**
165 :			** Relies on the ABS and VEC_* macros above, as well as math functions
166 :			** atan2(), sin(), cos(), sqrt(), and airCbrt(), defined as:
167 :
168 :			double
169 :			airCbrt(double v) {
170 :			#if defined(_WIN32) || defined(__STRICT_ANSI__)
171 :			return (v < 0.0 ? -pow(-v,1.0/3.0) : pow(v,1.0/3.0));
172 :			#else
173 :			return cbrt(v);
174 :			#endif
175 :			}
176 :
177 :			**
178 :			** HEY: the numerical precision issues here are very subtle, and
179 :			** merit some more scrutiny. With evals (1.000001, 1, 1), for example,
180 :			** whether it comes back as a single/double root, vs three distinct roots,
181 :			** is determines by the comparison between "D" and "epsilon", and the
182 :			** setting of epsilon seems pretty arbitrary at this point...
183 :			**
184 :			*/
185 :			int
186 :			evals(double eval[3],
187 :			const double _M00, const double _M01, const double _M02,
188 :			const double _M11, const double _M12,
189 :			const double _M22) {
190 :
191 :			#include "teigen-evals-A.c"
192 :
193 :			#include "teigen-evals-B.c"
194 :
195 :			return roots;
196 :			}
197 :
198 :			int
199 :			evals_evecs(double eval[3], double evec[9],
200 :			const double _M00, const double _M01, const double _M02,
201 :			const double _M11, const double _M12,
202 :			const double _M22) {
203 :			double r0[3], r1[3], r2[3], crs[3], len, dot;
204 :
205 :			#include "teigen-evals-A.c"
206 :
207 :			/* r0, r1, r2 are the vectors we manipulate to
208 :			find the nullspaces of M - lambda*I */
209 :			VEC_SET(r0, 0.0, M01, M02);
210 :			VEC_SET(r1, M01, 0.0, M12);
211 :			VEC_SET(r2, M02, M12, 0.0);
212 :			if (ROOT_THREE == roots) {
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 :			nullspace1(evec+3, r0, r1, r2);
217 :			r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
218 :			nullspace1(evec+6, r0, r1, r2);
219 :			} else if (ROOT_SINGLE_DOUBLE == roots) {
220 :			if (eval[1] == eval[2]) {
221 :			/* one big (eval[0]) , two small (eval[1,2]) */
222 :			r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
223 :			nullspace1(evec+0, r0, r1, r2);
224 :			r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
225 :			nullspace2(evec+3, evec+6, r0, r1, r2);
226 :			}
227 :			else {
228 :			/* two big (eval[0,1]), one small (eval[2]) */
229 :			r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
230 :			nullspace2(evec+0, evec+3, r0, r1, r2);
231 :			r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
232 :			nullspace1(evec+6, r0, r1, r2);
233 :			}
234 :			} else {
235 :			/* ROOT_TRIPLE == roots; use any basis for eigenvectors */
236 :			VEC_SET(evec+0, 1, 0, 0);
237 :			VEC_SET(evec+3, 0, 1, 0);
238 :			VEC_SET(evec+6, 0, 0, 1);
239 :			}
240 :			/* we always make sure its really orthonormal; keeping fixed the
241 :			eigenvector associated with the largest-magnitude eigenvalue */
242 :			if (ABS(eval[0]) > ABS(eval[2])) {
243 :			/* normalize evec+0 but don't move it */
244 :			VEC_NORM(evec+0, len);
245 :			dot = VEC_DOT(evec+0, evec+3); VEC_SCL_SUB(evec+3, dot, evec+0);
246 :			VEC_NORM(evec+3, len);
247 :			dot = VEC_DOT(evec+0, evec+6); VEC_SCL_SUB(evec+6, dot, evec+0);
248 :			dot = VEC_DOT(evec+3, evec+6); VEC_SCL_SUB(evec+6, dot, evec+3);
249 :			VEC_NORM(evec+6, len);
250 :			} else {
251 :			/* normalize evec+6 but don't move it */
252 :			VEC_NORM(evec+6, len);
253 :			dot = VEC_DOT(evec+6, evec+3); VEC_SCL_SUB(evec+3, dot, evec+6);
254 :			VEC_NORM(evec+3, len);
255 :			dot = VEC_DOT(evec+3, evec+0); VEC_SCL_SUB(evec+0, dot, evec+3);
256 :			dot = VEC_DOT(evec+6, evec+0); VEC_SCL_SUB(evec+0, dot, evec+6);
257 :			VEC_NORM(evec+0, len);
258 :			}
259 :
260 :			/* to be nice, make it right-handed */
261 :			VEC_CROSS(crs, evec+0, evec+3);
262 :			if (0 > VEC_DOT(crs, evec+6)) {
263 :			VEC_SCL(evec+6, -1);
264 :			}
265 :
266 :			#include "teigen-evals-B.c"
267 :
268 :			return roots;
269 :			}
270 :
271 :	jhr	1172	#endif /* PORTED */

---

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