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 1232 - (view) (download) (as text)

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

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