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

SCM Repository

[diderot] Annotation of /branches/pure-cfg/src/lib/c-target/eigen.c
ViewVC logotype

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

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