Bug Summary

File:src/ell/eigen.c
Location:line 688, column 3
Description:Assigned value is garbage or undefined

Annotated Source Code

1/*
2 Teem: Tools to process and visualize scientific data and images .
3 Copyright (C) 2013, 2012, 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 "ell.h"
26
27/* lop A
28 fprintf(stderr, "_ellAlign3: ----------\n");
29 fprintf(stderr, "_ellAlign3: v0 = %g %g %g\n", (v+0)[0], (v+0)[1], (v+0)[2]);
30 fprintf(stderr, "_ellAlign3: v3 = %g %g %g\n", (v+3)[0], (v+3)[1], (v+3)[2]);
31 fprintf(stderr, "_ellAlign3: v6 = %g %g %g\n", (v+6)[0], (v+6)[1], (v+6)[2]);
32 fprintf(stderr, "_ellAlign3: d = %g %g %g -> %d %d %d\n",
33 d0, d1, d2, Mi, ai, bi);
34 fprintf(stderr, "_ellAlign3: pre dot signs (03, 06, 36): %d %d %d\n",
35 airSgn(ELL_3V_DOT(v+0, v+3)),
36 airSgn(ELL_3V_DOT(v+0, v+6)),
37 airSgn(ELL_3V_DOT(v+3, v+6)));
38 */
39
40/* lop B
41 fprintf(stderr, "_ellAlign3: v0 = %g %g %g\n", (v+0)[0], (v+0)[1], (v+0)[2]);
42 fprintf(stderr, "_ellAlign3: v3 = %g %g %g\n", (v+3)[0], (v+3)[1], (v+3)[2]);
43 fprintf(stderr, "_ellAlign3: v6 = %g %g %g\n", (v+6)[0], (v+6)[1], (v+6)[2]);
44 fprintf(stderr, "_ellAlign3: post dot signs %d %d %d\n",
45 airSgn(ELL_3V_DOT(v+0, v+3)),
46 airSgn(ELL_3V_DOT(v+0, v+6)),
47 airSgn(ELL_3V_DOT(v+3, v+6)));
48 if (airSgn(ELL_3V_DOT(v+0, v+3)) < 0
49 || airSgn(ELL_3V_DOT(v+0, v+6)) < 0
50 || airSgn(ELL_3V_DOT(v+3, v+6)) < 0) {
51 exit(1);
52 }
53 */
54
55/*
56******** ell_quadratic()
57**
58** finds real roots of A*x^2 + B*x + C.
59**
60** records the found roots in the given root array, and returns a
61** value from the ell_quadratic_root* enum:
62**
63** ell_quadratic_root_two:
64** two distinct roots root[0] > root[1]
65** ell_quadratic_root_complex:
66** two complex conjugate roots at root[0] +/- i*root[1]
67** ell_quadratic_root_double:
68** a repeated root root[0] == root[1]
69**
70** HEY simple as this code may seem, it has definitely numerical
71** issues that have not been explored or fixed, such as what if A is
72** near 0. Also correctly handling the transition from double root to
73** complex roots needs to be re-thought, as well as this issue:
74** http://people.csail.mit.edu/bkph/articles/Quadratics.pdf Should
75** also understand http://www.cs.berkeley.edu/~wkahan/Qdrtcs.pdf
76**
77** This does NOT use biff
78*/
79int
80ell_quadratic(double root[2], double A, double B, double C) {
81 /* static const char me[]="ell_quadratic"; */
82 int ret;
83 double disc, rd, tmp, eps=1.0E-12;
84
85 disc = B*B - 4*A*C;
86 if (disc > 0) {
87 rd = sqrt(disc);
88 root[0] = (-B + rd)/(2*A);
89 root[1] = (-B - rd)/(2*A);
90 if (root[0] < root[1]) {
91 ELL_SWAP2(root[0], root[1], tmp)((tmp)=(root[0]),(root[0])=(root[1]),(root[1])=(tmp));
92 }
93 ret = ell_quadratic_root_two;
94 } else if (disc < -eps) {
95 root[0] = -B/(2*A);
96 root[1] = sqrt(-disc)/(2*A);
97 ret = ell_quadratic_root_complex;
98 } else {
99 /* 0 == disc or only *very slightly* negative */
100 root[0] = root[1] = -B/(2*A);
101 ret = ell_quadratic_root_double;
102 }
103 return ret;
104}
105
106int
107ell_2m_eigenvalues_d(double eval[2], const double m[4]) {
108 double A, B, C;
109 int ret;
110
111 A = 1;
112 B = -m[0] - m[3];
113 C = m[0]*m[3] - m[1]*m[2];
114 ret = ell_quadratic(eval, A, B, C);
115 return ret;
116}
117
118void
119ell_2m_1d_nullspace_d(double ans[2], const double _n[4]) {
120 /* static const char me[]="ell_2m_1d_nullspace_d"; */
121 double n[4], dot, len, rowv[2];
122
123 ELL_4V_COPY(n, _n)((n)[0] = (_n)[0], (n)[1] = (_n)[1], (n)[2] = (_n)[2], (n)[3]
= (_n)[3]);
124 dot = ELL_2V_DOT(n + 2*0, n + 2*1)((n + 2*0)[0]*(n + 2*1)[0] + (n + 2*0)[1]*(n + 2*1)[1]);
125 /*
126 fprintf(stderr, "!%s: n = {{%g,%g},{%g,%g}}\n", me,
127 n[0], n[1], n[2], n[3]);
128 fprintf(stderr, "!%s: dot = %g\n", me, dot);
129 */
130 if (dot > 0) {
131 ELL_2V_ADD2(rowv, n + 2*0, n + 2*1)((rowv)[0] = (n + 2*0)[0] + (n + 2*1)[0], (rowv)[1] = (n + 2*
0)[1] + (n + 2*1)[1]);
132 } else {
133 ELL_2V_SUB(rowv, n + 2*0, n + 2*1)((rowv)[0] = (n + 2*0)[0] - (n + 2*1)[0], (rowv)[1] = (n + 2*
0)[1] - (n + 2*1)[1]);
134 }
135 /* fprintf(stderr, "!%s: rowv = %g %g\n", me, rowv[0], rowv[1]); */
136 /* have found good description of what's perpendicular nullspace,
137 so now perpendicularize it */
138 ans[0] = rowv[1];
139 ans[1] = -rowv[0];
140 ELL_2V_NORM(ans, ans, len)(len = (sqrt((((ans))[0]*((ans))[0] + ((ans))[1]*((ans))[1]))
), ((ans)[0] = (1.0/len)*(ans)[0], (ans)[1] = (1.0/len)*(ans)
[1]));
141 /*
142 if (!(AIR_EXISTS(ans[0]) && AIR_EXISTS(ans[1]))) {
143 fprintf(stderr, "!%s: bad! %g %g\n", me, ans[0], ans[1]);
144 }
145 */
146 return;
147}
148
149/*
150******** ell_2m_eigensolve_d
151**
152** Eigensolve 2x2 matrix, which may be asymmetric
153*/
154int
155ell_2m_eigensolve_d(double eval[2], double evec[4], const double m[4]) {
156 /* static const char me[]="ell_2m_eigensolve_d"; */
157 double nul[4], ident[4] = {1,0,0,1};
158 int ret;
159
160 ret = ell_2m_eigenvalues_d(eval, m);
161 /*
162 fprintf(stderr, "!%s: m = {{%.17g,%.17g},{%.17g,%.17g}} -> "
163 "%s evals (%.17g,%.17g)\n", me, m[0], m[1], m[2], m[3],
164 airEnumStr(ell_quadratic_root, ret), eval[0], eval[1]);
165 */
166 switch (ret) {
167 case ell_quadratic_root_two:
168 ELL_4V_SCALE_ADD2(nul, 1.0, m, -eval[0], ident)((nul)[0] = (1.0)*(m)[0] + (-eval[0])*(ident)[0], (nul)[1] = (
1.0)*(m)[1] + (-eval[0])*(ident)[1], (nul)[2] = (1.0)*(m)[2] +
(-eval[0])*(ident)[2], (nul)[3] = (1.0)*(m)[3] + (-eval[0])*
(ident)[3]);
169 ell_2m_1d_nullspace_d(evec + 2*0, nul);
170 /*
171 fprintf(stderr, "!%s: eval=%.17g -> nul {{%.17g,%.17g},{%.17g,%.17g}} "
172 "-> evec %.17g %.17g\n", me, eval[0],
173 nul[0], nul[1], nul[2], nul[3],
174 (evec + 2*0)[0], (evec + 2*0)[1]);
175 */
176 ELL_4V_SCALE_ADD2(nul, 1.0, m, -eval[1], ident)((nul)[0] = (1.0)*(m)[0] + (-eval[1])*(ident)[0], (nul)[1] = (
1.0)*(m)[1] + (-eval[1])*(ident)[1], (nul)[2] = (1.0)*(m)[2] +
(-eval[1])*(ident)[2], (nul)[3] = (1.0)*(m)[3] + (-eval[1])*
(ident)[3]);
177 ell_2m_1d_nullspace_d(evec + 2*1, nul);
178 /*
179 fprintf(stderr, "!%s: eval=%.17g -> nul {{%.17g,%.17g},{%.17g,%.17g}} "
180 "-> evec %.17g %.17g\n", me, eval[1],
181 nul[0], nul[1], nul[2], nul[3],
182 (evec + 2*1)[0], (evec + 2*1)[1]);
183 */
184 break;
185 case ell_quadratic_root_double:
186 /* fprintf(stderr, "!%s: double eval=%.17g\n", me, eval[0]); */
187 ELL_4V_SCALE_ADD2(nul, 1.0, m, -eval[0], ident)((nul)[0] = (1.0)*(m)[0] + (-eval[0])*(ident)[0], (nul)[1] = (
1.0)*(m)[1] + (-eval[0])*(ident)[1], (nul)[2] = (1.0)*(m)[2] +
(-eval[0])*(ident)[2], (nul)[3] = (1.0)*(m)[3] + (-eval[0])*
(ident)[3]);
188 /*
189 fprintf(stderr, "!%s: nul = {{%.17g,%.17g},{%.17g,%.17g}} (len %.17g)\n",
190 me, nul[0], nul[1], nul[2], nul[3], ELL_4V_LEN(nul));
191 */
192 if (ELL_4V_DOT(nul, nul)((nul)[0]*(nul)[0] + (nul)[1]*(nul)[1] + (nul)[2]*(nul)[2] + (
nul)[3]*(nul)[3])) {
193 /* projecting out the nullspace produced non-zero matrix,
194 (possibly from an asymmetric matrix) so there is real
195 orientation to recover */
196 ell_2m_1d_nullspace_d(evec + 2*0, nul);
197 ELL_2V_COPY(evec + 2*1, evec + 2*0)((evec + 2*1)[0] = (evec + 2*0)[0], (evec + 2*1)[1] = (evec +
2*0)[1]);
198 } else {
199 /* so this was isotropic symmetric; invent orientation */
200 ELL_2V_SET(evec + 2*0, 1, 0)((evec + 2*0)[0]=(1), (evec + 2*0)[1]=(0));
201 ELL_2V_SET(evec + 2*1, 0, 1)((evec + 2*1)[0]=(0), (evec + 2*1)[1]=(1));
202 }
203 break;
204 case ell_quadratic_root_complex:
205 /* HEY punting for now */
206 ELL_2V_SET(evec + 2*0, 0.5, 0)((evec + 2*0)[0]=(0.5), (evec + 2*0)[1]=(0));
207 ELL_2V_SET(evec + 2*1, 0, 0.5)((evec + 2*1)[0]=(0), (evec + 2*1)[1]=(0.5));
208 break;
209 default:
210 /* fprintf(stderr, "%s: unexpected solution indicator %d\n", me, ret); */
211 break;
212 }
213 return ret;
214}
215
216
217void
218_ell_align3_d(double v[9]) {
219 double d0, d1, d2;
220 int Mi, ai, bi;
221
222 d0 = ELL_3V_DOT(v+0, v+0)((v+0)[0]*(v+0)[0] + (v+0)[1]*(v+0)[1] + (v+0)[2]*(v+0)[2]);
223 d1 = ELL_3V_DOT(v+3, v+3)((v+3)[0]*(v+3)[0] + (v+3)[1]*(v+3)[1] + (v+3)[2]*(v+3)[2]);
224 d2 = ELL_3V_DOT(v+6, v+6)((v+6)[0]*(v+6)[0] + (v+6)[1]*(v+6)[1] + (v+6)[2]*(v+6)[2]);
225 Mi = ELL_MAX3_IDX(d0, d1, d2)(d0 > d1 ? (d1 > d2 ? 0 : (d0 > d2 ? 0 : 2)) : (d2 >
d1 ? 2 : 1));
226 ai = (Mi + 1) % 3;
227 bi = (Mi + 2) % 3;
228 /* lop A */
229 if (ELL_3V_DOT(v+3*Mi, v+3*ai)((v+3*Mi)[0]*(v+3*ai)[0] + (v+3*Mi)[1]*(v+3*ai)[1] + (v+3*Mi)
[2]*(v+3*ai)[2]) < 0) {
230 ELL_3V_SCALE(v+3*ai, -1, v+3*ai)((v+3*ai)[0] = (-1)*(v+3*ai)[0], (v+3*ai)[1] = (-1)*(v+3*ai)[
1], (v+3*ai)[2] = (-1)*(v+3*ai)[2]);
231 }
232 if (ELL_3V_DOT(v+3*Mi, v+3*bi)((v+3*Mi)[0]*(v+3*bi)[0] + (v+3*Mi)[1]*(v+3*bi)[1] + (v+3*Mi)
[2]*(v+3*bi)[2]) < 0) {
233 ELL_3V_SCALE(v+3*bi, -1, v+3*bi)((v+3*bi)[0] = (-1)*(v+3*bi)[0], (v+3*bi)[1] = (-1)*(v+3*bi)[
1], (v+3*bi)[2] = (-1)*(v+3*bi)[2]);
234 }
235 /* lob B */
236 /* we can't guarantee that dot(v+3*ai,v+3*bi) > 0 . . . */
237}
238
239/*
240** leaves v+3*0 untouched, but makes sure that v+3*0, v+3*1, and v+3*2
241** are mutually orthogonal. Also leaves the magnitudes of all
242** vectors unchanged.
243*/
244void
245_ell_3m_enforce_orthogonality(double v[9]) {
246 double d00, d10, d11, d20, d21, d22, scl, tv[3];
247
248 d00 = ELL_3V_DOT(v+3*0, v+3*0)((v+3*0)[0]*(v+3*0)[0] + (v+3*0)[1]*(v+3*0)[1] + (v+3*0)[2]*(
v+3*0)[2]);
249 d10 = ELL_3V_DOT(v+3*1, v+3*0)((v+3*1)[0]*(v+3*0)[0] + (v+3*1)[1]*(v+3*0)[1] + (v+3*1)[2]*(
v+3*0)[2]);
250 d11 = ELL_3V_DOT(v+3*1, v+3*1)((v+3*1)[0]*(v+3*1)[0] + (v+3*1)[1]*(v+3*1)[1] + (v+3*1)[2]*(
v+3*1)[2]);
251 ELL_3V_SCALE_ADD2(tv, 1, v+3*1, -d10/d00, v+3*0)((tv)[0] = (1)*(v+3*1)[0] + (-d10/d00)*(v+3*0)[0], (tv)[1] = (
1)*(v+3*1)[1] + (-d10/d00)*(v+3*0)[1], (tv)[2] = (1)*(v+3*1)[
2] + (-d10/d00)*(v+3*0)[2]);
252 scl = sqrt(d11/ELL_3V_DOT(tv, tv)((tv)[0]*(tv)[0] + (tv)[1]*(tv)[1] + (tv)[2]*(tv)[2]));
253 ELL_3V_SCALE(v+3*1, scl, tv)((v+3*1)[0] = (scl)*(tv)[0], (v+3*1)[1] = (scl)*(tv)[1], (v+3
*1)[2] = (scl)*(tv)[2]);
254 d20 = ELL_3V_DOT(v+3*2, v+3*0)((v+3*2)[0]*(v+3*0)[0] + (v+3*2)[1]*(v+3*0)[1] + (v+3*2)[2]*(
v+3*0)[2]);
255 d21 = ELL_3V_DOT(v+3*2, v+3*1)((v+3*2)[0]*(v+3*1)[0] + (v+3*2)[1]*(v+3*1)[1] + (v+3*2)[2]*(
v+3*1)[2]);
256 d22 = ELL_3V_DOT(v+3*2, v+3*2)((v+3*2)[0]*(v+3*2)[0] + (v+3*2)[1]*(v+3*2)[1] + (v+3*2)[2]*(
v+3*2)[2]);
257 ELL_3V_SCALE_ADD3(tv, 1, v+3*2, -d20/d00, v+3*0, -d21/d00, v+3*1)((tv)[0] = (1)*(v+3*2)[0] + (-d20/d00)*(v+3*0)[0] + (-d21/d00
)*(v+3*1)[0], (tv)[1] = (1)*(v+3*2)[1] + (-d20/d00)*(v+3*0)[1
] + (-d21/d00)*(v+3*1)[1], (tv)[2] = (1)*(v+3*2)[2] + (-d20/d00
)*(v+3*0)[2] + (-d21/d00)*(v+3*1)[2]);
258 scl = sqrt(d22/ELL_3V_DOT(tv, tv)((tv)[0]*(tv)[0] + (tv)[1]*(tv)[1] + (tv)[2]*(tv)[2]));
259 ELL_3V_SCALE(v+3*2, scl, tv)((v+3*2)[0] = (scl)*(tv)[0], (v+3*2)[1] = (scl)*(tv)[1], (v+3
*2)[2] = (scl)*(tv)[2]);
260 return;
261}
262
263/*
264** makes sure that v+3*2 has a positive dot product with
265** cross product of v+3*0 and v+3*1
266*/
267void
268_ell_3m_make_right_handed_d(double v[9]) {
269 double x[3];
270
271 ELL_3V_CROSS(x, v+3*0, v+3*1)((x)[0] = (v+3*0)[1]*(v+3*1)[2] - (v+3*0)[2]*(v+3*1)[1], (x)[
1] = (v+3*0)[2]*(v+3*1)[0] - (v+3*0)[0]*(v+3*1)[2], (x)[2] = (
v+3*0)[0]*(v+3*1)[1] - (v+3*0)[1]*(v+3*1)[0]);
272 if (0 > ELL_3V_DOT(x, v+3*2)((x)[0]*(v+3*2)[0] + (x)[1]*(v+3*2)[1] + (x)[2]*(v+3*2)[2])) {
273 ELL_3V_SCALE(v+3*2, -1, v+3*2)((v+3*2)[0] = (-1)*(v+3*2)[0], (v+3*2)[1] = (-1)*(v+3*2)[1], (
v+3*2)[2] = (-1)*(v+3*2)[2]);
274 }
275}
276
277/* lop A
278 fprintf(stderr, "=== pre ===\n");
279 fprintf(stderr, "crosses: %g %g %g\n", (t+0)[0], (t+0)[1], (t+0)[2]);
280 fprintf(stderr, " %g %g %g\n", (t+3)[0], (t+3)[1], (t+3)[2]);
281 fprintf(stderr, " %g %g %g\n", (t+6)[0], (t+6)[1], (t+6)[2]);
282 fprintf(stderr, "cross dots: %g %g %g\n",
283 ELL_3V_DOT(t+0, t+3), ELL_3V_DOT(t+0, t+6), ELL_3V_DOT(t+3, t+6));
284*/
285
286/* lop B
287 fprintf(stderr, "=== post ===\n");
288 fprintf(stderr, "crosses: %g %g %g\n", (t+0)[0], (t+0)[1], (t+0)[2]);
289 fprintf(stderr, " %g %g %g\n", (t+3)[0], (t+3)[1], (t+3)[2]);
290 fprintf(stderr, " %g %g %g\n", (t+6)[0], (t+6)[1], (t+6)[2]);
291 fprintf(stderr, "cross dots: %g %g %g\n",
292 ELL_3V_DOT(t+0, t+3), ELL_3V_DOT(t+0, t+6), ELL_3V_DOT(t+3, t+6));
293*/
294
295/*
296******** ell_3m_1d_nullspace_d()
297**
298** the given matrix is assumed to have a nullspace of dimension one.
299** A normalized vector which spans the nullspace is put into ans.
300**
301** The given nullspace matrix is NOT modified.
302**
303** This does NOT use biff
304*/
305void
306ell_3m_1d_nullspace_d(double ans[3], const double _n[9]) {
307 double t[9], n[9], norm;
308
309 ELL_3M_TRANSPOSE(n, _n)((n)[0] = (_n)[0], (n)[1] = (_n)[3], (n)[2] = (_n)[6], (n)[3]
= (_n)[1], (n)[4] = (_n)[4], (n)[5] = (_n)[7], (n)[6] = (_n)
[2], (n)[7] = (_n)[5], (n)[8] = (_n)[8]);
310 /* find the three cross-products of pairs of column vectors of n */
311 ELL_3V_CROSS(t+0, n+0, n+3)((t+0)[0] = (n+0)[1]*(n+3)[2] - (n+0)[2]*(n+3)[1], (t+0)[1] =
(n+0)[2]*(n+3)[0] - (n+0)[0]*(n+3)[2], (t+0)[2] = (n+0)[0]*(
n+3)[1] - (n+0)[1]*(n+3)[0]);
312 ELL_3V_CROSS(t+3, n+0, n+6)((t+3)[0] = (n+0)[1]*(n+6)[2] - (n+0)[2]*(n+6)[1], (t+3)[1] =
(n+0)[2]*(n+6)[0] - (n+0)[0]*(n+6)[2], (t+3)[2] = (n+0)[0]*(
n+6)[1] - (n+0)[1]*(n+6)[0]);
313 ELL_3V_CROSS(t+6, n+3, n+6)((t+6)[0] = (n+3)[1]*(n+6)[2] - (n+3)[2]*(n+6)[1], (t+6)[1] =
(n+3)[2]*(n+6)[0] - (n+3)[0]*(n+6)[2], (t+6)[2] = (n+3)[0]*(
n+6)[1] - (n+3)[1]*(n+6)[0]);
314 /* lop A */
315 _ell_align3_d(t);
316 /* lop B */
317 /* add them up (longer, hence more accurate, should dominate) */
318 ELL_3V_ADD3(ans, t+0, t+3, t+6)((ans)[0] = (t+0)[0] + (t+3)[0] + (t+6)[0], (ans)[1] = (t+0)[
1] + (t+3)[1] + (t+6)[1], (ans)[2] = (t+0)[2] + (t+3)[2] + (t
+6)[2]);
319
320 /* normalize */
321 ELL_3V_NORM(ans, ans, norm)(norm = (sqrt((((ans))[0]*((ans))[0] + ((ans))[1]*((ans))[1] +
((ans))[2]*((ans))[2]))), ((ans)[0] = (1.0/norm)*(ans)[0], (
ans)[1] = (1.0/norm)*(ans)[1], (ans)[2] = (1.0/norm)*(ans)[2]
));
322
323 return;
324}
325
326/*
327******** ell_3m_2d_nullspace_d()
328**
329** the given matrix is assumed to have a nullspace of dimension two.
330**
331** The given nullspace matrix is NOT modified
332**
333** This does NOT use biff
334*/
335void
336ell_3m_2d_nullspace_d(double ans0[3], double ans1[3], const double _n[9]) {
337 double n[9], tmp[3], norm;
338
339 ELL_3M_TRANSPOSE(n, _n)((n)[0] = (_n)[0], (n)[1] = (_n)[3], (n)[2] = (_n)[6], (n)[3]
= (_n)[1], (n)[4] = (_n)[4], (n)[5] = (_n)[7], (n)[6] = (_n)
[2], (n)[7] = (_n)[5], (n)[8] = (_n)[8]);
340 _ell_align3_d(n);
341 ELL_3V_ADD3(tmp, n+0, n+3, n+6)((tmp)[0] = (n+0)[0] + (n+3)[0] + (n+6)[0], (tmp)[1] = (n+0)[
1] + (n+3)[1] + (n+6)[1], (tmp)[2] = (n+0)[2] + (n+3)[2] + (n
+6)[2]);
342 ELL_3V_NORM(tmp, tmp, norm)(norm = (sqrt((((tmp))[0]*((tmp))[0] + ((tmp))[1]*((tmp))[1] +
((tmp))[2]*((tmp))[2]))), ((tmp)[0] = (1.0/norm)*(tmp)[0], (
tmp)[1] = (1.0/norm)*(tmp)[1], (tmp)[2] = (1.0/norm)*(tmp)[2]
));
343
344 /* any two vectors which are perpendicular to the (supposedly 1D)
345 span of the column vectors span the nullspace */
346 ell_3v_perp_d(ans0, tmp);
347 ELL_3V_NORM(ans0, ans0, norm)(norm = (sqrt((((ans0))[0]*((ans0))[0] + ((ans0))[1]*((ans0))
[1] + ((ans0))[2]*((ans0))[2]))), ((ans0)[0] = (1.0/norm)*(ans0
)[0], (ans0)[1] = (1.0/norm)*(ans0)[1], (ans0)[2] = (1.0/norm
)*(ans0)[2]));
348 ELL_3V_CROSS(ans1, tmp, ans0)((ans1)[0] = (tmp)[1]*(ans0)[2] - (tmp)[2]*(ans0)[1], (ans1)[
1] = (tmp)[2]*(ans0)[0] - (tmp)[0]*(ans0)[2], (ans1)[2] = (tmp
)[0]*(ans0)[1] - (tmp)[1]*(ans0)[0]);
349
350 return;
351}
352
353/*
354******** ell_3m_eigenvalues_d()
355**
356** finds eigenvalues of given matrix.
357**
358** returns information about the roots according to ellCubeRoot enum,
359** see header for ellCubic for details.
360**
361** given matrix is NOT modified
362**
363** This does NOT use biff
364**
365** Doing the frobenius normalization proved successfull in avoiding the
366** the creating of NaN eigenvalues when the coefficients of the matrix
367** were really large (> 50000). Also, when the matrix norm was really
368** small, the comparison to "epsilon" in ell_cubic mistook three separate
369** roots for a single and a double, with this matrix in particular:
370** 1.7421892 0.0137642 0.0152975
371** 0.0137642 1.7565432 -0.0062296
372** 0.0152975 -0.0062296 1.7700019
373** (actually, this is prior to tenEigensolve's isotropic removal)
374**
375** HEY: tenEigensolve_d and tenEigensolve_f start by removing the
376** isotropic part of the tensor. It may be that that smarts should
377** be migrated here, but GLK is uncertain how it would change the
378** handling of non-symmetric matrices.
379*/
380int
381ell_3m_eigenvalues_d(double _eval[3], const double _m[9], const int newton) {
382 double A, B, C, scale, frob, m[9], eval[3];
383 int roots;
384
385 frob = ELL_3M_FROB(_m)(sqrt((((_m)+0)[0]*((_m)+0)[0] + ((_m)+0)[1]*((_m)+0)[1] + ((
_m)+0)[2]*((_m)+0)[2]) + (((_m)+3)[0]*((_m)+3)[0] + ((_m)+3)[
1]*((_m)+3)[1] + ((_m)+3)[2]*((_m)+3)[2]) + (((_m)+6)[0]*((_m
)+6)[0] + ((_m)+6)[1]*((_m)+6)[1] + ((_m)+6)[2]*((_m)+6)[2]))
);
386 scale = frob ? 1.0/frob : 1.0;
387 ELL_3M_SCALE(m, scale, _m)((((m)+0)[0] = ((scale))*((_m)+0)[0], ((m)+0)[1] = ((scale))*
((_m)+0)[1], ((m)+0)[2] = ((scale))*((_m)+0)[2]), (((m)+3)[0]
= ((scale))*((_m)+3)[0], ((m)+3)[1] = ((scale))*((_m)+3)[1],
((m)+3)[2] = ((scale))*((_m)+3)[2]), (((m)+6)[0] = ((scale))
*((_m)+6)[0], ((m)+6)[1] = ((scale))*((_m)+6)[1], ((m)+6)[2] =
((scale))*((_m)+6)[2]));
388 /*
389 printf("!%s: m = %g %g %g; %g %g %g; %g %g %g\n", "ell_3m_eigenvalues_d",
390 m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8]);
391 */
392 /*
393 ** from gordon with mathematica; these are the coefficients of the
394 ** cubic polynomial in x: det(x*I - M). The full cubic is
395 ** x^3 + A*x^2 + B*x + C.
396 */
397 A = -m[0] - m[4] - m[8];
398 B = m[0]*m[4] - m[3]*m[1]
399 + m[0]*m[8] - m[6]*m[2]
400 + m[4]*m[8] - m[7]*m[5];
401 C = (m[6]*m[4] - m[3]*m[7])*m[2]
402 + (m[0]*m[7] - m[6]*m[1])*m[5]
403 + (m[3]*m[1] - m[0]*m[4])*m[8];
404 /*
405 printf("!%s: A B C = %g %g %g\n", "ell_3m_eigenvalues_d", A, B, C);
406 */
407 roots = ell_cubic(eval, A, B, C, newton);
408 /* no longer need to sort here */
409 ELL_3V_SCALE(_eval, 1.0/scale, eval)((_eval)[0] = (1.0/scale)*(eval)[0], (_eval)[1] = (1.0/scale)
*(eval)[1], (_eval)[2] = (1.0/scale)*(eval)[2]);
410 return roots;
411}
412
413void
414_ell_3m_evecs_d(double evec[9], double eval[3], int roots,
415 const double m[9]) {
416 double n[9], e0=0, e1=0.0, e2=0.0, t /* , tmpv[3] */ ;
417
418 ELL_3V_GET(e0, e1, e2, eval)((e0) = (eval)[0], (e1) = (eval)[1], (e2) = (eval)[2]);
419 /* if (ell_debug) {
420 printf("ell_3m_evecs_d: numroots = %d\n", numroots);
421 } */
422
423 /* we form m - lambda*I by doing a memcpy from m, and then
424 (repeatedly) over-writing the diagonal elements */
425 ELL_3M_COPY(n, m)((((n)+0)[0] = ((m)+0)[0], ((n)+0)[1] = ((m)+0)[1], ((n)+0)[2
] = ((m)+0)[2]), (((n)+3)[0] = ((m)+3)[0], ((n)+3)[1] = ((m)+
3)[1], ((n)+3)[2] = ((m)+3)[2]), (((n)+6)[0] = ((m)+6)[0], ((
n)+6)[1] = ((m)+6)[1], ((n)+6)[2] = ((m)+6)[2]));
426 switch (roots) {
427 case ell_cubic_root_three:
428 /* if (ell_debug) {
429 printf("ell_3m_evecs_d: evals: %20.15f %20.15f %20.15f\n",
430 eval[0], eval[1], eval[2]);
431 } */
432 ELL_3M_DIAG_SET(n, m[0]-e0, m[4]-e0, m[8]-e0)((n)[0] = (m[0]-e0), (n)[4] = (m[4]-e0), (n)[8] = (m[8]-e0));
433 ell_3m_1d_nullspace_d(evec+0, n);
434 ELL_3M_DIAG_SET(n, m[0]-e1, m[4]-e1, m[8]-e1)((n)[0] = (m[0]-e1), (n)[4] = (m[4]-e1), (n)[8] = (m[8]-e1));
435 ell_3m_1d_nullspace_d(evec+3, n);
436 ELL_3M_DIAG_SET(n, m[0]-e2, m[4]-e2, m[8]-e2)((n)[0] = (m[0]-e2), (n)[4] = (m[4]-e2), (n)[8] = (m[8]-e2));
437 ell_3m_1d_nullspace_d(evec+6, n);
438 _ell_3m_enforce_orthogonality(evec);
439 _ell_3m_make_right_handed_d(evec);
440 ELL_3V_SET(eval, e0, e1, e2)((eval)[0] = (e0), (eval)[1] = (e1), (eval)[2] = (e2));
441 break;
442 case ell_cubic_root_single_double:
443 ELL_SORT3(e0, e1, e2, t)if (e0 > e1) { if (e1 < e2) { if (e0 > e2) { ((t)=(e1
),(e1)=(e2),(e2)=(t)); } else { ((t)=(e0),(e0)=(e2),(e2)=(e1)
,(e1)=(t)); } } } else { if (e1 > e2) { if (e0 > e2) { (
(t)=(e0),(e0)=(e1),(e1)=(t)); } else { ((t)=(e0),(e0)=(e1),(e1
)=(e2),(e2)=(t)); } } else { ((t)=(e0),(e0)=(e2),(e2)=(t)); }
};
444 if (e0 > e1) {
445 /* one big (e0) , two small (e1, e2) : more like a cigar */
446 ELL_3M_DIAG_SET(n, m[0]-e0, m[4]-e0, m[8]-e0)((n)[0] = (m[0]-e0), (n)[4] = (m[4]-e0), (n)[8] = (m[8]-e0));
447 ell_3m_1d_nullspace_d(evec+0, n);
448 ELL_3M_DIAG_SET(n, m[0]-e1, m[4]-e1, m[8]-e1)((n)[0] = (m[0]-e1), (n)[4] = (m[4]-e1), (n)[8] = (m[8]-e1));
449 ell_3m_2d_nullspace_d(evec+3, evec+6, n);
450 }
451 else {
452 /* two big (e0, e1), one small (e2): more like a pancake */
453 ELL_3M_DIAG_SET(n, m[0]-e0, m[4]-e0, m[8]-e0)((n)[0] = (m[0]-e0), (n)[4] = (m[4]-e0), (n)[8] = (m[8]-e0));
454 ell_3m_2d_nullspace_d(evec+0, evec+3, n);
455 ELL_3M_DIAG_SET(n, m[0]-e2, m[4]-e2, m[8]-e2)((n)[0] = (m[0]-e2), (n)[4] = (m[4]-e2), (n)[8] = (m[8]-e2));
456 ell_3m_1d_nullspace_d(evec+6, n);
457 }
458 _ell_3m_enforce_orthogonality(evec);
459 _ell_3m_make_right_handed_d(evec);
460 ELL_3V_SET(eval, e0, e1, e2)((eval)[0] = (e0), (eval)[1] = (e1), (eval)[2] = (e2));
461 break;
462 case ell_cubic_root_triple:
463 /* one triple root; use any basis as the eigenvectors */
464 ELL_3V_SET(evec+0, 1, 0, 0)((evec+0)[0] = (1), (evec+0)[1] = (0), (evec+0)[2] = (0));
465 ELL_3V_SET(evec+3, 0, 1, 0)((evec+3)[0] = (0), (evec+3)[1] = (1), (evec+3)[2] = (0));
466 ELL_3V_SET(evec+6, 0, 0, 1)((evec+6)[0] = (0), (evec+6)[1] = (0), (evec+6)[2] = (1));
467 ELL_3V_SET(eval, e0, e1, e2)((eval)[0] = (e0), (eval)[1] = (e1), (eval)[2] = (e2));
468 break;
469 case ell_cubic_root_single:
470 /* only one real root */
471 ELL_3M_DIAG_SET(n, m[0]-e0, m[4]-e0, m[8]-e0)((n)[0] = (m[0]-e0), (n)[4] = (m[4]-e0), (n)[8] = (m[8]-e0));
472 ell_3m_1d_nullspace_d(evec+0, n);
473 ELL_3V_SET(evec+3, AIR_NAN, AIR_NAN, AIR_NAN)((evec+3)[0] = ((airFloatQNaN.f)), (evec+3)[1] = ((airFloatQNaN
.f)), (evec+3)[2] = ((airFloatQNaN.f)));
474 ELL_3V_SET(evec+6, AIR_NAN, AIR_NAN, AIR_NAN)((evec+6)[0] = ((airFloatQNaN.f)), (evec+6)[1] = ((airFloatQNaN
.f)), (evec+6)[2] = ((airFloatQNaN.f)));
475 ELL_3V_SET(eval, e0, AIR_NAN, AIR_NAN)((eval)[0] = (e0), (eval)[1] = ((airFloatQNaN.f)), (eval)[2] =
((airFloatQNaN.f)));
476 break;
477 }
478 /* if (ell_debug) {
479 printf("ell_3m_evecs_d (numroots = %d): evecs: \n", numroots);
480 ELL_3MV_MUL(tmpv, m, evec[0]);
481 printf(" (%g:%g): %20.15f %20.15f %20.15f\n",
482 eval[0], ELL_3V_DOT(evec[0], tmpv),
483 evec[0][0], evec[0][1], evec[0][2]);
484 ELL_3MV_MUL(tmpv, m, evec[1]);
485 printf(" (%g:%g): %20.15f %20.15f %20.15f\n",
486 eval[1], ELL_3V_DOT(evec[1], tmpv),
487 evec[1][0], evec[1][1], evec[1][2]);
488 ELL_3MV_MUL(tmpv, m, evec[2]);
489 printf(" (%g:%g): %20.15f %20.15f %20.15f\n",
490 eval[2], ELL_3V_DOT(evec[2], tmpv),
491 evec[2][0], evec[2][1], evec[2][2]);
492 } */
493 return;
494}
495
496/*
497******** ell_3m_eigensolve_d()
498**
499** finds eigenvalues and eigenvectors of given matrix m
500**
501** returns information about the roots according to ellCubeRoot enum,
502** see header for ellCubic for details. When eval[i] is set, evec+3*i
503** is set to a corresponding eigenvector. The eigenvectors are
504** (evec+0)[], (evec+3)[], and (evec+6)[]
505**
506** NOTE: Even in the post-Teem-1.7 switch from column-major to
507** row-major- its still the case that the eigenvectors are at
508** evec+0, evec+3, evec+6: this means that they USED to be the
509** "columns" of the matrix, and NOW they're the rows.
510**
511** The eigenvalues (and associated eigenvectors) are sorted in
512** descending order.
513**
514** This does NOT use biff
515*/
516int
517ell_3m_eigensolve_d(double eval[3], double evec[9],
518 const double m[9], const int newton) {
519 int roots;
520
521 /* if (ell_debug) {
522 printf("ell_3m_eigensolve_d: input matrix:\n");
523 printf("{{%20.15f,\t%20.15f,\t%20.15f},\n", m[0], m[1], m[2]);
524 printf(" {%20.15f,\t%20.15f,\t%20.15f},\n", m[3], m[4], m[5]);
525 printf(" {%20.15f,\t%20.15f,\t%20.15f}};\n",m[6], m[7], m[8]);
526 } */
527
528 roots = ell_3m_eigenvalues_d(eval, m, newton);
529 _ell_3m_evecs_d(evec, eval, roots, m);
530
531 return roots;
532}
533
534/* ____________________________ 3m2sub ____________________________ */
535/*
536******** ell_3m2sub_eigenvalues_d
537**
538** for doing eigensolve of the upper-left 2x2 matrix sub-matrix of a
539** 3x3 matrix. The other entries are assumed to be zero. A 0 root is
540** put last (in eval[2]), possibly in defiance of the usual eigenvalue
541** ordering.
542*/
543int
544ell_3m2sub_eigenvalues_d(double eval[3], const double _m[9]) {
545 double A, B, m[4], D, Dsq, eps=1.0E-11;
546 int roots;
547 /* static const char me[]="ell_3m2sub_eigenvalues_d"; */
548
549 m[0] = _m[0];
550 m[1] = _m[1];
551 m[2] = _m[3];
552 m[3] = _m[4];
553
554 /* cubic characteristic equation is L^3 + A*L^2 + B*L = 0 */
555 A = -m[0] - m[3];
556 B = m[0]*m[3] - m[1]*m[2];
557 Dsq = A*A - 4*B;
558 /*
559 fprintf(stderr, "!%s: m = {{%f,%f},{%f,%f}} -> A=%f B=%f Dsq=%.17f %s 0 (%.17f)\n", me,
560 m[0], m[1], m[2], m[3], A, B, Dsq,
561 (Dsq > 0 ? ">" : (Dsq < 0 ? "<" : "==")), eps);
562 fprintf(stderr, "!%s: Dsq = \n", me);
563 airFPFprintf_d(stderr, Dsq);
564 fprintf(stderr, "!%s: eps = \n", me);
565 airFPFprintf_d(stderr, eps);
566 ell_3m_print_d(stderr, _m);
567 */
568 if (Dsq > eps) {
569 D = sqrt(Dsq);
570 eval[0] = (-A + D)/2;
571 eval[1] = (-A - D)/2;
572 eval[2] = 0;
573 roots = ell_cubic_root_three;
574 } else if (Dsq < -eps) {
575 /* no quadratic roots; only the implied zero */
576 ELL_3V_SET(eval, AIR_NAN, AIR_NAN, 0)((eval)[0] = ((airFloatQNaN.f)), (eval)[1] = ((airFloatQNaN.f
)), (eval)[2] = (0));
577 roots = ell_cubic_root_single;
578 } else {
579 /* a quadratic double root */
580 ELL_3V_SET(eval, -A/2, -A/2, 0)((eval)[0] = (-A/2), (eval)[1] = (-A/2), (eval)[2] = (0));
581 roots = ell_cubic_root_single_double;
582 }
583 /*
584 fprintf(stderr, "!%s: Dsq=%f, roots=%d (%f %f %f)\n",
585 me, Dsq, roots, eval[0], eval[1], eval[2]);
586 */
587 return roots;
588}
589
590void
591_ell_22v_enforce_orthogonality(double uu[2], double _vv[2]) {
592 double dot, vv[2], len;
593
594 dot = ELL_2V_DOT(uu, _vv)((uu)[0]*(_vv)[0] + (uu)[1]*(_vv)[1]);
595 ELL_2V_SCALE_ADD2(vv, 1, _vv, -dot, uu)((vv)[0] = (1)*(_vv)[0] + (-dot)*(uu)[0], (vv)[1] = (1)*(_vv)
[1] + (-dot)*(uu)[1]);
596 ELL_2V_NORM(_vv, vv, len)(len = (sqrt((((vv))[0]*((vv))[0] + ((vv))[1]*((vv))[1]))), (
(_vv)[0] = (1.0/len)*(vv)[0], (_vv)[1] = (1.0/len)*(vv)[1]));
597 return;
598}
599
600/*
601** NOTE: assumes that eval and roots have come from
602** ell_3m2sub_eigenvalues_d(m)
603*/
604void
605_ell_3m2sub_evecs_d(double evec[9], double eval[3], int roots,
606 const double m[9]) {
607 double n[4];
608 static const char me[]="_ell_3m2sub_evecs_d";
609
610 if (ell_cubic_root_three == roots) {
611 /* set off-diagonal entries once */
612 n[1] = m[1];
613 n[2] = m[3];
614 /* find first evec */
615 n[0] = m[0] - eval[0];
616 n[3] = m[4] - eval[0];
617 ell_2m_1d_nullspace_d(evec + 3*0, n);
618 (evec + 3*0)[2] = 0;
619 /* find second evec */
620 n[0] = m[0] - eval[1];
621 n[3] = m[4] - eval[1];
622 ell_2m_1d_nullspace_d(evec + 3*1, n);
623 (evec + 3*1)[2] = 0;
624 _ell_22v_enforce_orthogonality(evec + 3*0, evec + 3*1);
625 /* make right-handed */
626 ELL_3V_CROSS(evec + 3*2, evec + 3*0, evec + 3*1)((evec + 3*2)[0] = (evec + 3*0)[1]*(evec + 3*1)[2] - (evec + 3
*0)[2]*(evec + 3*1)[1], (evec + 3*2)[1] = (evec + 3*0)[2]*(evec
+ 3*1)[0] - (evec + 3*0)[0]*(evec + 3*1)[2], (evec + 3*2)[2]
= (evec + 3*0)[0]*(evec + 3*1)[1] - (evec + 3*0)[1]*(evec + 3
*1)[0]);
627 } else if (ell_cubic_root_single_double == roots) {
628 /* can pick any 2D basis */
629 ELL_3V_SET(evec + 3*0, 1, 0, 0)((evec + 3*0)[0] = (1), (evec + 3*0)[1] = (0), (evec + 3*0)[2
] = (0));
630 ELL_3V_SET(evec + 3*1, 0, 1, 0)((evec + 3*1)[0] = (0), (evec + 3*1)[1] = (1), (evec + 3*1)[2
] = (0));
631 ELL_3V_SET(evec + 3*2, 0, 0, 1)((evec + 3*2)[0] = (0), (evec + 3*2)[1] = (0), (evec + 3*2)[2
] = (1));
632 } else {
633 /* ell_cubic_root_single == roots, if assumptions are met */
634 ELL_3V_SET(evec + 3*0, AIR_NAN, AIR_NAN, 0)((evec + 3*0)[0] = ((airFloatQNaN.f)), (evec + 3*0)[1] = ((airFloatQNaN
.f)), (evec + 3*0)[2] = (0));
635 ELL_3V_SET(evec + 3*1, AIR_NAN, AIR_NAN, 0)((evec + 3*1)[0] = ((airFloatQNaN.f)), (evec + 3*1)[1] = ((airFloatQNaN
.f)), (evec + 3*1)[2] = (0));
636 ELL_3V_SET(evec + 3*2, 0, 0, 1)((evec + 3*2)[0] = (0), (evec + 3*2)[1] = (0), (evec + 3*2)[2
] = (1));
637 }
638 if (!ELL_3M_EXISTS(evec)(((((int)(!((((evec) + 0)[0]) - (((evec) + 0)[0]))))) &&
(((int)(!((((evec) + 0)[1]) - (((evec) + 0)[1]))))) &&
(((int)(!((((evec) + 0)[2]) - (((evec) + 0)[2])))))) &&
((((int)(!((((evec) + 3)[0]) - (((evec) + 3)[0]))))) &&
(((int)(!((((evec) + 3)[1]) - (((evec) + 3)[1]))))) &&
(((int)(!((((evec) + 3)[2]) - (((evec) + 3)[2])))))) &&
((((int)(!((((evec) + 6)[0]) - (((evec) + 6)[0]))))) &&
(((int)(!((((evec) + 6)[1]) - (((evec) + 6)[1]))))) &&
(((int)(!((((evec) + 6)[2]) - (((evec) + 6)[2])))))))) {
639 fprintf(stderr__stderrp, "%s: given m = \n", me);
640 ell_3m_print_d(stderr__stderrp, m);
641 fprintf(stderr__stderrp, "%s: got roots = %s (%d) and evecs = \n", me,
642 airEnumStr(ell_cubic_root, roots), roots);
643 ell_3m_print_d(stderr__stderrp, evec);
644 }
645 return;
646}
647
648int
649ell_3m2sub_eigensolve_d(double eval[3], double evec[9],
650 const double m[9]) {
651 int roots;
652
653 roots = ell_3m2sub_eigenvalues_d(eval, m);
654 _ell_3m2sub_evecs_d(evec, eval, roots, m);
655
656 return roots;
657}
658
659/* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 3m2sub ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
660
661/*
662******** ell_3m_svd_d
663**
664** singular value decomposition:
665** mat = uu * diag(sval) * vv
666**
667** singular values are square roots of eigenvalues of mat * mat^T
668** columns of uu are eigenvectors of mat * mat^T
669** rows of vv are eigenvectors of mat^T * mat
670**
671** returns info about singular values according to ellCubeRoot enum
672**
673** HEY: I think this does the wrong thing when given a symmetric
674** matrix with negative eigenvalues . . .
675*/
676int
677ell_3m_svd_d(double uu[9], double sval[3], double vv[9],
678 const double mat[9], const int newton) {
679 double trn[9], msqr[9], eval[3], evec[9];
680 int roots;
681
682 ELL_3M_TRANSPOSE(trn, mat)((trn)[0] = (mat)[0], (trn)[1] = (mat)[3], (trn)[2] = (mat)[6
], (trn)[3] = (mat)[1], (trn)[4] = (mat)[4], (trn)[5] = (mat)
[7], (trn)[6] = (mat)[2], (trn)[7] = (mat)[5], (trn)[8] = (mat
)[8]);
683 ELL_3M_MUL(msqr, mat, trn)((msqr)[0] = (mat)[0]*(trn)[0] + (mat)[1]*(trn)[3] + (mat)[2]
*(trn)[6], (msqr)[1] = (mat)[0]*(trn)[1] + (mat)[1]*(trn)[4] +
(mat)[2]*(trn)[7], (msqr)[2] = (mat)[0]*(trn)[2] + (mat)[1]*
(trn)[5] + (mat)[2]*(trn)[8], (msqr)[3] = (mat)[3]*(trn)[0] +
(mat)[4]*(trn)[3] + (mat)[5]*(trn)[6], (msqr)[4] = (mat)[3]*
(trn)[1] + (mat)[4]*(trn)[4] + (mat)[5]*(trn)[7], (msqr)[5] =
(mat)[3]*(trn)[2] + (mat)[4]*(trn)[5] + (mat)[5]*(trn)[8], (
msqr)[6] = (mat)[6]*(trn)[0] + (mat)[7]*(trn)[3] + (mat)[8]*(
trn)[6], (msqr)[7] = (mat)[6]*(trn)[1] + (mat)[7]*(trn)[4] + (
mat)[8]*(trn)[7], (msqr)[8] = (mat)[6]*(trn)[2] + (mat)[7]*(trn
)[5] + (mat)[8]*(trn)[8]);
684 roots = ell_3m_eigensolve_d(eval, evec, msqr, newton);
685 sval[0] = sqrt(eval[0]);
686 sval[1] = sqrt(eval[1]);
687 sval[2] = sqrt(eval[2]);
688 ELL_3M_TRANSPOSE(uu, evec)((uu)[0] = (evec)[0], (uu)[1] = (evec)[3], (uu)[2] = (evec)[6
], (uu)[3] = (evec)[1], (uu)[4] = (evec)[4], (uu)[5] = (evec)
[7], (uu)[6] = (evec)[2], (uu)[7] = (evec)[5], (uu)[8] = (evec
)[8]);
Within the expansion of the macro 'ELL_3M_TRANSPOSE':
a
Assigned value is garbage or undefined
689 ELL_3M_MUL(msqr, trn, mat)((msqr)[0] = (trn)[0]*(mat)[0] + (trn)[1]*(mat)[3] + (trn)[2]
*(mat)[6], (msqr)[1] = (trn)[0]*(mat)[1] + (trn)[1]*(mat)[4] +
(trn)[2]*(mat)[7], (msqr)[2] = (trn)[0]*(mat)[2] + (trn)[1]*
(mat)[5] + (trn)[2]*(mat)[8], (msqr)[3] = (trn)[3]*(mat)[0] +
(trn)[4]*(mat)[3] + (trn)[5]*(mat)[6], (msqr)[4] = (trn)[3]*
(mat)[1] + (trn)[4]*(mat)[4] + (trn)[5]*(mat)[7], (msqr)[5] =
(trn)[3]*(mat)[2] + (trn)[4]*(mat)[5] + (trn)[5]*(mat)[8], (
msqr)[6] = (trn)[6]*(mat)[0] + (trn)[7]*(mat)[3] + (trn)[8]*(
mat)[6], (msqr)[7] = (trn)[6]*(mat)[1] + (trn)[7]*(mat)[4] + (
trn)[8]*(mat)[7], (msqr)[8] = (trn)[6]*(mat)[2] + (trn)[7]*(mat
)[5] + (trn)[8]*(mat)[8]);
690 _ell_3m_evecs_d(vv, eval, roots, msqr);
691
692 return roots;
693}
694
695/*
696** NOTE: profiling showed that about a quarter of the execution time of
697** ell_6ms_eigensolve_d() is spent here; so reconsider its need and
698** implementation . . . (fabs vs. AIR_ABS() made no difference)
699*/
700static void
701_maxI_sum_find(unsigned int maxI[2], double *sumon, double *sumoff,
702 double mat[6][6]) {
703 double maxm, tmp;
704 unsigned int rrI, ccI;
705
706 /* we hope that all these loops are unrolled by the optimizer */
707 *sumon = *sumoff = 0.0;
708 for (rrI=0; rrI<6; rrI++) {
709 *sumon += AIR_ABS(mat[rrI][rrI])((mat[rrI][rrI]) > 0.0f ? (mat[rrI][rrI]) : -(mat[rrI][rrI
]));
710 }
711 maxm = -1;
712 maxI[0] = maxI[1] = 0;
713 for (rrI=0; rrI<5; rrI++) {
714 for (ccI=rrI+1; ccI<6; ccI++) {
715 tmp = AIR_ABS(mat[rrI][ccI])((mat[rrI][ccI]) > 0.0f ? (mat[rrI][ccI]) : -(mat[rrI][ccI
]));
716 *sumoff += tmp;
717 if (tmp > maxm) {
718 maxm = tmp;
719 maxI[0] = rrI;
720 maxI[1] = ccI;
721 }
722 }
723 }
724
725 /*
726 if (1) {
727 double nrm, trc;
728 nrm = trc = 0;
729 for (rrI=0; rrI<6; rrI++) {
730 trc += mat[rrI][rrI];
731 nrm += mat[rrI][rrI]*mat[rrI][rrI];
732 }
733 for (rrI=0; rrI<5; rrI++) {
734 for (ccI=rrI+1; ccI<6; ccI++) {
735 nrm += 2*mat[rrI][ccI]*mat[rrI][ccI];
736 }
737 }
738 fprintf(stderr, "---------------- invars = %g %g\n", trc, nrm);
739 }
740 */
741
742 return;
743}
744
745static int
746_compar(const void *A_void, const void *B_void) {
747 const double *A, *B;
748 A = AIR_CAST(const double *, A_void)((const double *)(A_void));
749 B = AIR_CAST(const double *, B_void)((const double *)(B_void));
750 return (A[0] < B[0] ? 1 : (A[0] > B[0] ? -1 : 0));
751}
752
753/*
754******* ell_6ms_eigensolve_d
755**
756** uses Jacobi iterations to find eigensystem of 6x6 symmetric matrix,
757** given in sym[21], to within convergence threshold eps. Puts
758** eigenvalues, in descending order, in eval[6], and corresponding
759** eigenvectors in _evec+0, _evec+6, . . ., _evec+30. NOTE: you can
760** pass a NULL _evec if eigenvectors aren't needed.
761**
762** does NOT use biff
763*/
764int
765ell_6ms_eigensolve_d(double eval[6], double _evec[36],
766 const double sym[21], const double eps) {
767 /* char me[]="ell_6ms_eigensolve_d"; */
768 double mat[2][6][6], evec[2][6][6], sumon, sumoff, evtmp[12];
769 unsigned int cur, rrI, ccI, maxI[2], iter;
770
771 if (!( eval && sym && eps >= 0 )) {
772 return 1;
773 }
774 /* copy symmetric matrix sym[] into upper tris of mat[0][][] & mat[1][][] */
775 mat[0][0][0] = sym[ 0];
776 mat[0][0][1] = sym[ 1];
777 mat[0][0][2] = sym[ 2];
778 mat[0][0][3] = sym[ 3];
779 mat[0][0][4] = sym[ 4];
780 mat[0][0][5] = sym[ 5];
781 mat[0][1][1] = sym[ 6];
782 mat[0][1][2] = sym[ 7];
783 mat[0][1][3] = sym[ 8];
784 mat[0][1][4] = sym[ 9];
785 mat[0][1][5] = sym[10];
786 mat[0][2][2] = sym[11];
787 mat[0][2][3] = sym[12];
788 mat[0][2][4] = sym[13];
789 mat[0][2][5] = sym[14];
790 mat[0][3][3] = sym[15];
791 mat[0][3][4] = sym[16];
792 mat[0][3][5] = sym[17];
793 mat[0][4][4] = sym[18];
794 mat[0][4][5] = sym[19];
795 mat[0][5][5] = sym[20];
796 if (_evec) {
797 /* initialize evec[0]; */
798 for (rrI=0; rrI<6; rrI++) {
799 for (ccI=0; ccI<6; ccI++) {
800 evec[0][ccI][rrI] = (rrI == ccI);
801 }
802 }
803 }
804 /*
805 fprintf(stderr, "!%s(INIT): m = [", me);
806 for (rrI=0; rrI<6; rrI++) {
807 for (ccI=0; ccI<6; ccI++) {
808 fprintf(stderr, "%f%s",
809 (rrI <= ccI ? mat[0][rrI][ccI] : mat[0][ccI][rrI]),
810 ccI<5 ? "," : (rrI<5 ? ";" : "]"));
811 }
812 fprintf(stderr, "\n");
813 }
814 */
815 maxI[0] = maxI[1] = UINT_MAX(2147483647 *2U +1U); /* quiet warnings about using maxI unset */
816 _maxI_sum_find(maxI, &sumon, &sumoff, mat[0]);
817 cur = 1; /* fake out anticipating first line of loop */
818 iter = 0;
819 while (sumoff/sumon > eps) {
820 double th, tt, cc, ss;
821 const unsigned int P = maxI[0];
822 const unsigned int Q = maxI[1];
823 //make sure that P and Q are within the bounds for mat[2][6][6]
824 if(P >=6 || Q >= 6){
825 break;
826 }
827 /*
828 fprintf(stderr, "!%s(%u): sumoff/sumon = %g/%g = %g > %g\n", me, iter,
829 sumoff, sumon, sumoff/sumon, eps);
830 */
831 cur = 1 - cur;
832
833 th = (mat[cur][Q][Q] - mat[cur][P][P])/(2*mat[cur][P][Q]);
834 tt = (th > 0 ? +1 : -1)/(AIR_ABS(th)((th) > 0.0f ? (th) : -(th)) + sqrt(th*th + 1));
835 cc = 1/sqrt(tt*tt + 1);
836 ss = cc*tt;
837 /*
838 fprintf(stderr, "!%s(%u): maxI = (P,Q) = (%u,%u) --> ss=%f, cc=%f\n",
839 me, iter, P, Q, ss, cc);
840 fprintf(stderr, " r = [");
841 for (rrI=0; rrI<6; rrI++) {
842 for (ccI=0; ccI<6; ccI++) {
843 fprintf(stderr, "%g%s",
844 (rrI == ccI
845 ? (rrI == P || rrI == Q ? cc : 1.0)
846 : (rrI == P && ccI == Q
847 ? ss
848 : (rrI == Q && ccI == P
849 ? -ss
850 : 0))),
851 ccI<5 ? "," : (rrI<5 ? ";" : "]"));
852 }
853 fprintf(stderr, "\n");
854 }
855 */
856 /* initialize by copying whole matrix */
857 for (rrI=0; rrI<6; rrI++) {
858 for (ccI=rrI; ccI<6; ccI++) {
859 mat[1-cur][rrI][ccI] = mat[cur][rrI][ccI];
860 }
861 }
862 /* perform Jacobi rotation */
863 for (rrI=0; rrI<P; rrI++) {
864 mat[1-cur][rrI][P] = cc*mat[cur][rrI][P] - ss*mat[cur][rrI][Q];
865 }
866 for (ccI=P+1; ccI<6; ccI++) {
867 mat[1-cur][P][ccI] = cc*mat[cur][P][ccI] - ss*(Q <= ccI
868 ? mat[cur][Q][ccI]
869 : mat[cur][ccI][Q]);
870 }
871 for (rrI=0; rrI<Q; rrI++) {
872 mat[1-cur][rrI][Q] = ss*(rrI <= P
873 ? mat[cur][rrI][P]
874 : mat[cur][P][rrI]) + cc*mat[cur][rrI][Q];
875 }
876 for (ccI=Q+1; ccI<6; ccI++) {
877 mat[1-cur][Q][ccI] = ss*mat[cur][P][ccI] + cc*mat[cur][Q][ccI];
878 }
879 /* set special entries */
880 mat[1-cur][P][P] = mat[cur][P][P] - tt*mat[cur][P][Q];
881 mat[1-cur][Q][Q] = mat[cur][Q][Q] + tt*mat[cur][P][Q];
882 mat[1-cur][P][Q] = 0.0;
883 if (_evec) {
884 /* NOTE: the eigenvectors use transpose of indexing of mat */
885 /* start by copying all */
886 for (rrI=0; rrI<6; rrI++) {
887 for (ccI=0; ccI<6; ccI++) {
888 evec[1-cur][ccI][rrI] = evec[cur][ccI][rrI];
889 }
890 }
891 for (rrI=0; rrI<6; rrI++) {
892 evec[1-cur][P][rrI] = cc*evec[cur][P][rrI] - ss*evec[cur][Q][rrI];
893 evec[1-cur][Q][rrI] = ss*evec[cur][P][rrI] + cc*evec[cur][Q][rrI];
894 }
895 }
896
897 _maxI_sum_find(maxI, &sumon, &sumoff, mat[1-cur]);
898
899 /*
900 fprintf(stderr, "!%s(%u): m = [", me, iter);
901 for (rrI=0; rrI<6; rrI++) {
902 for (ccI=0; ccI<6; ccI++) {
903 fprintf(stderr, "%f%s",
904 (rrI <= ccI ? mat[1-cur][rrI][ccI] : mat[1-cur][ccI][rrI]),
905 ccI<5 ? "," : (rrI<5 ? ";" : "]"));
906 }
907 fprintf(stderr, "\n");
908 }
909 */
910 iter++;
911 }
912 /* 1-cur is index of final solution */
913
914 /* sort evals */
915 for (ccI=0; ccI<6; ccI++) {
916 evtmp[0 + 2*ccI] = mat[1-cur][ccI][ccI];
917 evtmp[1 + 2*ccI] = ccI;
918 }
919 qsort(evtmp, 6, 2*sizeof(double), _compar);
920
921 /* copy out solution */
922 for (ccI=0; ccI<6; ccI++) {
923 eval[ccI] = evtmp[0 + 2*ccI];
924 if (_evec) {
925 unsigned eeI;
926 for (rrI=0; rrI<6; rrI++) {
927 eeI = AIR_CAST(unsigned int, evtmp[1 + 2*ccI])((unsigned int)(evtmp[1 + 2*ccI]));
928 _evec[rrI + 6*ccI] = evec[1-cur][eeI][rrI];
929 }
930 }
931 }
932
933 return 0;
934}