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

SCM Repository

[diderot] Annotation of /branches/pure-cfg/src/compiler/cl-target/fragments/eigen3x3.in
ViewVC logotype

Annotation of /branches/pure-cfg/src/compiler/cl-target/fragments/eigen3x3.in

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1664 - (view) (download)

1 : jhr 1658 #define ROOT_TRIPLE 2 /* ell_cubic_root_triple */
2 :     #define ROOT_SINGLE_DOUBLE 3 /* ell_cubic_root_single_double */
3 :     #define ROOT_THREE 4 /* ell_cubic_root_three */
4 :    
5 :     // from http://en.wikipedia.org/wiki/Square_root_of_3
6 :     #define M_SQRT3 1.732050807568877293527446341506
7 :    
8 :     /*
9 :     ** All the three given vectors span only a 2D space, and this finds
10 :     ** the normal to that plane. Simply sums up all the pair-wise
11 :     ** cross-products to get a good estimate. Trick is getting the cross
12 :     ** products to line up before summing.
13 :     */
14 :     STATIC_INLINE Diderot_vec3_t nullspace1 (
15 :     const Diderot_vec3_t r0,
16 :     const Diderot_vec3_t r1,
17 :     const Diderot_vec3_t r2)
18 :     {
19 :     Diderot_vec3_t ret, crs;
20 :    
21 :     ret = cross(r0, r1);
22 :     crs = cross(r1, r2);
23 :     /* ret += crs or ret -= crs; whichever makes ret longer */
24 :     if (dot(ret, crs) > 0.0)
25 :     ret += crs;
26 :     else
27 :     ret -= crs;
28 :    
29 :     crs = cross(r0, r2);
30 :     /* ret += crs or ret -= crs; whichever makes ret longer */
31 :     if (dot(ret, crs) > 0.0)
32 :     ret += crs;
33 :     else
34 :     ret -= crs;
35 :    
36 :     return ret;
37 :     }
38 :    
39 :     /*
40 :     ** All vectors are in the same 1D space, we have to find two
41 :     ** mutually vectors perpendicular to that span
42 :     */
43 :     static void nullspace2 (Diderot_vec3_t rets[2],
44 :     const Diderot_vec3_t r0,
45 :     const Diderot_vec3_t r1,
46 :     const Diderot_vec3_t r2)
47 :     {
48 :     Diderot_vec3_t sqr, sum;
49 :     int idx;
50 :    
51 :     sum = r0;
52 :     if (dot(sum, r1) > 0)
53 :     sum += r1;
54 :     else
55 :     sum -= r1;
56 :     if (dot(sum, r2) > 0)
57 :     sum += r2;
58 :     else
59 :     sum -= r2;
60 :     // find largest component, to get most stable expression for a perpendicular vector
61 :     sqr = sum*sum;
62 :     idx = 0;
63 :     if ((sqr.s0 > sqr.s1) && (sqr.s0 > sqr.s2)) {
64 :     // sum[0] is largest component
65 :     rets[0] = VEC3(sum.s1 - sum.s2, -sum.s0, sum.s0);
66 :     }
67 :     else if (sqr.s1 > sqr.s2) {
68 :     // sum[0] is largest component
69 :     rets[0] = VEC3(-sum.s1, sum.s0 - sum.s2, sum.s1);
70 :     }
71 :     else {
72 :     // sum[0] is largest component
73 :     rets[0] = VEC3(-sum.s2, sum.s2, sum.s0 - sum.s1);
74 :     }
75 :    
76 :     rets[1] = cross(rets[0], sum);
77 :     return;
78 :     }
79 :    
80 :     /*
81 :     ** Eigensolver for symmetric 3x3 matrix:
82 :     **
83 :     ** M00 M01 M02
84 :     ** M01 M11 M12
85 :     ** M02 M12 M22
86 :     **
87 :     ** Must be passed pointer to eval vector, and will compute eigenvalues
88 :     ** only if evec[9] is non-NULL. Computed eigenvectors are at evec+0,
89 :     ** evec+3, and evec+6.
90 :     **
91 :     ** Return value indicates something about the eigenvalue solution to
92 :     ** the cubic characteristic equation; see ROOT_ #defines above
93 :     **
94 :     ** HEY: the numerical precision issues here are very subtle, and
95 :     ** merit some more scrutiny. With evals (1.000001, 1, 1), for example,
96 :     ** whether it comes back as a single/double root, vs three distinct roots,
97 :     ** is determines by the comparison between "D" and "epsilon", and the
98 :     ** setting of epsilon seems pretty arbitrary at this point...
99 :     **
100 :     */
101 :     int Diderot_evals3x3 (
102 : lamonts 1662 float3 * eval,
103 : jhr 1658 const Diderot_real_t _M00, const Diderot_real_t _M01, const Diderot_real_t _M02,
104 :     const Diderot_real_t _M11, const Diderot_real_t _M12,
105 :     const Diderot_real_t _M22)
106 :     {
107 :     Diderot_real_t mean, norm, rnorm, Q, R, QQQ, D, theta, M00, M01, M02, M11, M12, M22;
108 :     int roots;
109 :    
110 :     /* copy the given matrix elements */
111 :     M00 = _M00;
112 :     M01 = _M01;
113 :     M02 = _M02;
114 :     M11 = _M11;
115 :     M12 = _M12;
116 :     M22 = _M22;
117 :    
118 :     /*
119 :     ** subtract out the eigenvalue mean (will add back to evals later);
120 :     ** helps with numerical stability
121 :     */
122 :     mean = (M00 + M11 + M22)/3.0;
123 :     M00 -= mean;
124 :     M11 -= mean;
125 :     M22 -= mean;
126 :    
127 :     /*
128 :     ** divide out L2 norm of eigenvalues (will multiply back later);
129 :     ** this too seems to help with stability
130 :     */
131 :     norm = sqrt(M00*M00 + 2*M01*M01 + 2*M02*M02 +
132 :     M11*M11 + 2*M12*M12 +
133 :     M22*M22);
134 :     // QUESTION: what if norm is very close to zero?
135 :     rnorm = norm ? 1.0/norm : 1.0;
136 :     M00 *= rnorm;
137 :     M01 *= rnorm;
138 :     M02 *= rnorm;
139 :     M11 *= rnorm;
140 :     M12 *= rnorm;
141 :     M22 *= rnorm;
142 :    
143 :     /* this code is a mix of prior Teem code and ideas from Eberly's
144 :     * "Eigensystems for 3 x 3 Symmetric Matrices (Revisited)"
145 :     */
146 :     Q = (M01*M01 + M02*M02 + M12*M12 - M00*M11 - M00*M22 - M11*M22)/3.0;
147 :     QQQ = Q*Q*Q;
148 :     R = (M00*M11*M22 + M02*(2*M01*M12 - M02*M11)
149 :     - M00*M12*M12 - M01*M01*M22)/2.0;
150 :     D = QQQ - R*R;
151 :     if (D > EPSILON) {
152 :     /* three distinct roots- this is the most common case */
153 :     double mm, ss, cc;
154 : lamonts 1664 theta = atan2(sqrt(D), R)/3.0;
155 : jhr 1658 mm = sqrt(Q);
156 : lamonts 1664 ss = sin(theta);
157 :     cc = cos(theta);
158 : lamonts 1662 eval->s0 = 2*mm*cc;
159 :     eval->s1 = mm*(-cc + M_SQRT3*ss);
160 :     eval->s2 = mm*(-cc - M_SQRT3*ss);
161 : jhr 1658 roots = ROOT_THREE;
162 :     }
163 :     /* else D is near enough to zero */
164 :     else if ((R < -EPSILON) || (EPSILON < R)) {
165 :     double U;
166 :     /* one double root and one single root */
167 :     U = cbrt(R); /* cube root function */
168 :     if (U > 0) {
169 : lamonts 1662 eval->s0 = 2*U;
170 :     eval->s1 = -U;
171 :     eval->s2 = -U;
172 : jhr 1658 }
173 :     else {
174 : lamonts 1662 eval->s0 = -U;
175 :     eval->s1 = -U;
176 :     eval->s2 = 2*U;
177 : jhr 1658 }
178 :     roots = ROOT_SINGLE_DOUBLE;
179 :     }
180 :     else {
181 :     /* a triple root! */
182 : lamonts 1662 eval->s0 = eval->s1 = eval->s2 = 0.0;
183 : jhr 1658 roots = ROOT_TRIPLE;
184 :     }
185 :    
186 :     /* multiply back by eigenvalue L2 norm */
187 : lamonts 1662 eval->s0 /= rnorm;
188 :     eval->s1 /= rnorm;
189 :     eval->s2 /= rnorm;
190 : jhr 1658
191 :     /* add back in the eigenvalue mean */
192 : lamonts 1662 eval->s0 += mean;
193 :     eval->s1 += mean;
194 :     eval->s2 += mean;
195 : jhr 1658
196 :     return roots;
197 :     }
198 :    
199 :     int Diderot_evecs3x3 (
200 : lamonts 1662 float3 * eval, Diderot_vec3_t evecs[3],
201 : jhr 1658 const Diderot_real_t _M00, const Diderot_real_t _M01, const Diderot_real_t _M02,
202 :     const Diderot_real_t _M11, const Diderot_real_t _M12,
203 :     const Diderot_real_t _M22)
204 :     {
205 :     Diderot_real_t len, dot;
206 :    
207 :     Diderot_real_t mean, norm, rnorm, Q, R, QQQ, D, theta, M00, M01, M02, M11, M12, M22;
208 :     Diderot_real_t epsilon = 1.0E-12;
209 :     int roots;
210 :    
211 :     /* copy the given matrix elements */
212 :     M00 = _M00;
213 :     M01 = _M01;
214 :     M02 = _M02;
215 :     M11 = _M11;
216 :     M12 = _M12;
217 :     M22 = _M22;
218 :    
219 :     /*
220 :     ** subtract out the eigenvalue mean (will add back to evals later);
221 :     ** helps with numerical stability
222 :     */
223 :     mean = (M00 + M11 + M22)/3.0;
224 :     M00 -= mean;
225 :     M11 -= mean;
226 :     M22 -= mean;
227 :    
228 :     /*
229 :     ** divide out L2 norm of eigenvalues (will multiply back later);
230 :     ** this too seems to help with stability
231 :     */
232 :     norm = sqrt(M00*M00 + 2*M01*M01 + 2*M02*M02 +
233 :     M11*M11 + 2*M12*M12 +
234 :     M22*M22);
235 :     rnorm = norm ? 1.0/norm : 1.0;
236 :     M00 *= rnorm;
237 :     M01 *= rnorm;
238 :     M02 *= rnorm;
239 :     M11 *= rnorm;
240 :     M12 *= rnorm;
241 :     M22 *= rnorm;
242 :    
243 :     /* this code is a mix of prior Teem code and ideas from Eberly's
244 :     * "Eigensystems for 3 x 3 Symmetric Matrices (Revisited)"
245 :     */
246 :     Q = (M01*M01 + M02*M02 + M12*M12 - M00*M11 - M00*M22 - M11*M22)/3.0;
247 :     QQQ = Q*Q*Q;
248 :     R = (M00*M11*M22 + M02*(2*M01*M12 - M02*M11) - M00*M12*M12 - M01*M01*M22)/2.0;
249 :     D = QQQ - R*R;
250 :     if (D > epsilon) {
251 :     /* three distinct roots- this is the most common case */
252 :     double mm, ss, cc;
253 : lamonts 1664 theta = atan2(sqrt(D), R)/3.0;
254 : jhr 1658 mm = sqrt(Q);
255 : lamonts 1664 ss = sin(theta);
256 :     cc = cos(theta);
257 : lamonts 1662 eval->s0 = 2*mm*cc;
258 :     eval->s1 = mm*(-cc + sqrt(3.0)*ss);
259 :     eval->s2 = mm*(-cc - sqrt(3.0)*ss);
260 : jhr 1658 roots = ROOT_THREE;
261 :     }
262 :     /* else D is near enough to zero */
263 :     else if (R < -epsilon || epsilon < R) {
264 :     double U;
265 :     /* one double root and one single root */
266 :     U = cbrt(R); /* cube root function */
267 :     if (U > 0) {
268 : lamonts 1662 eval->s0 = 2*U;
269 :     eval->s1 = -U;
270 :     eval->s2 = -U;
271 : jhr 1658 } else {
272 : lamonts 1662 eval->s0 = -U;
273 :     eval->s1 = -U;
274 :     eval->s2 = 2*U;
275 : jhr 1658 }
276 :     roots = ROOT_SINGLE_DOUBLE;
277 :     }
278 :     else {
279 :     /* a triple root! */
280 : lamonts 1662 eval->s0 = eval->s1 = eval->s2 = 0.0;
281 : jhr 1658 roots = ROOT_TRIPLE;
282 :     }
283 :    
284 : lamonts 1662 Diderot_vec3_t ev = VEC3(eval->s0, eval->s1, eval->s2);
285 : jhr 1658 if (ROOT_THREE == roots) {
286 :     evecs[0] = nullspace1 (
287 : lamonts 1664 VEC3(M00 - eval->s0, M01, M02),
288 :     VEC3(M01, M11 - eval->s0, M12),
289 :     VEC3(M02, M12, M22 - eval->s0));
290 : jhr 1658 evecs[1] = nullspace1 (
291 : lamonts 1664 VEC3(M00 - eval->s1, M01, M02),
292 :     VEC3(M01, M11 - eval->s1, M12),
293 :     VEC3(M02, M12, M22 - eval->s1));
294 : jhr 1658 evecs[2] = nullspace1 (
295 : lamonts 1664 VEC3(M00 - eval->s2, M01, M02),
296 :     VEC3(M01, M11 - eval->s2, M12),
297 :     VEC3(M02, M12, M22 - eval->s2));
298 : jhr 1658 }
299 :     else if (ROOT_SINGLE_DOUBLE == roots) {
300 : lamonts 1662 if (eval->s1 == eval->s2) {
301 : jhr 1658 /* one big (eval[0]) , two small (eval[1,2]) */
302 :     evecs[0] = nullspace1 (
303 : lamonts 1664 VEC3(M00 - eval->s0, M01, M02),
304 :     VEC3(M01, M11 - eval->s0, M12),
305 :     VEC3(M02, M12, M22 - eval->s0));
306 : jhr 1658 nullspace2 (evecs+1,
307 : lamonts 1664 VEC3(M00 - eval->s1, M01, M02),
308 :     VEC3(M01, M11 - eval->s1, M12),
309 :     VEC3(M02, M12, M22 - eval->s1));
310 : jhr 1658 }
311 :     else {
312 :     /* two big (eval[0,1]), one small (eval[2]) */
313 :     nullspace2 (evecs,
314 : lamonts 1664 VEC3(M00 - eval->s0, M01, M02),
315 :     VEC3(M01, M11 - eval->s0, M12),
316 :     VEC3(M02, M12, M22 - eval->s0));
317 : jhr 1658 evecs[2] = nullspace1 (
318 : lamonts 1664 VEC3(M00 - eval->s2, M01, M02),
319 :     VEC3(M01, M11 - eval->s2, M12),
320 :     VEC3(M02, M12, M22 - eval->s2));
321 : jhr 1658 }
322 :     }
323 :     else {
324 :     /* ROOT_TRIPLE == roots; use any basis for eigenvectors */
325 :     evecs[0] = VEC3(1, 0, 0);
326 :     evecs[1] = VEC3(0, 1, 0);
327 :     evecs[2] = VEC3(0, 0, 1);
328 :     }
329 :     /* we always make sure it's really orthonormal; keeping fixed the
330 :     * eigenvector associated with the largest-magnitude eigenvalue
331 :     */
332 : lamonts 1664 if (fabs(eval->s0) > fabs(eval->s2)) {
333 : jhr 1658 /* normalize evec[0] but don't move it */
334 :     evecs[0] = normalize(evecs[0]);
335 :     evecs[1] -= dot(evecs[1], evecs[0]) * evecs[0];
336 :     evecs[1] = normalize(evecs[1]);
337 :     evecs[2] = cross(evecs[0], evecs[1]);
338 :     }
339 :     else {
340 :     /* normalize evec[2] but don't move it */
341 :     evecs[2] = normalize(evecs[2]);
342 :     evecs[1] -= dot(evecs[1], evecs[2]) * evecs[2];
343 :     evecs[1] = normalize(evecs[1]);
344 :     evecs[0] = cross(evecs[1], evecs[2]);
345 :     }
346 :     /* note that the right-handedness check has been folded into
347 :     the code above to enforce orthogonality. Indeed, some work
348 :     could be removed by never really bothering to find all three
349 :     eigenvectors; just find two and then use the cross-product.
350 :     The next iteration of the code will do that */
351 :    
352 :     /* multiply back by eigenvalue L2 norm */
353 : lamonts 1662 eval->s0 /= rnorm;
354 :     eval->s1 /= rnorm;
355 :     eval->s2 /= rnorm;
356 : jhr 1658
357 :     /* add back in the eigenvalue mean */
358 : lamonts 1662 eval->s0 += mean;
359 :     eval->s1 += mean;
360 :     eval->s2 += mean;
361 : jhr 1658
362 :     return roots;
363 :     }

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