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

Parent Directory | Revision Log | Patch

	revision 1183, Wed May 11 12:59:21 2011 UTC	revision 1184, Wed May 11 12:59:41 2011 UTC
#	Line 102	Line 102
102	Diderot_vec3_t sqr, sum;	Diderot_vec3_t sqr, sum;
103	int idx;	int idx;
104
105	VEC_COPY(sum, r0);	sum = r0;
106	if (VEC_DOT(sum, r1) > 0) {	if (dot(sum, r1) > 0)
107	VEC_ADD(sum, r1);	sum += r1;
108	} else {	else
109	VEC_SUB(sum, r1);	sum -= r1;
110	}	if (dot(sum, r2) > 0)
111	if (VEC_DOT(sum, r2) > 0) {	sum += r2;
112	VEC_ADD(sum, r2);	else
113	} else {	sum -= r2;
114	VEC_SUB(sum, r2);	/* find largest component, to get most stable expression for a perpendicular vector */
115	}	sqr = sum*sum;
	/* find largest component, to get most stable expression for a
	perpendicular vector */
	sqr[0] = sum[0]*sum[0];
	sqr[1] = sum[1]*sum[1];
	sqr[2] = sum[2]*sum[2];
116	idx = 0;	idx = 0;
117	if (sqr[0] < sqr[1])	if (sqr[0] < sqr[1])
118	idx = 1;	idx = 1;
119	if (sqr[idx] < sqr[2])	if (sqr[idx] < sqr[2])
120	idx = 2;	idx = 2;
121	/* reta will be perpendicular to sum */	/* rets[0] will be perpendicular to sum */
122	if (0 == idx) {	if (0 == idx) {
123	rets[0] = vec3(sum[1] - sum[2], -sum[0], sum[0]);	rets[0] = vec3(sum[1] - sum[2], -sum[0], sum[0]);
124	} else if (1 == idx) {	} else if (1 == idx) {
#	Line 131	Line 126
126	} else {	} else {
127	rets[0] = vec3(-sum[2], sum[2], sum[0] - sum[1]);	rets[0] = vec3(-sum[2], sum[2], sum[0] - sum[1]);
128	}	}
129	/* and now retb will be perpendicular to both reta and sum */	/* and now rets[1] will be perpendicular to both rets[0] and sum */
130	rets[1] = cross(rets[0], sum);	rets[1] = cross(rets[0], sum);
131	return;	return;
132	}	}
#	Line 318	Line 313
313	evec[1] = vec3(0, 1, 0);	evec[1] = vec3(0, 1, 0);
314	evec[2] = vec3(0, 0, 1);	evec[2] = vec3(0, 0, 1);
315	}	}
316	/* we always make sure its really orthonormal; keeping fixed the	/* we always make sure it's really orthonormal; keeping fixed the
317	* eigenvector associated with the largest-magnitude eigenvalue	* eigenvector associated with the largest-magnitude eigenvalue
318	*/	*/
319	if (ABS(eval[0]) > ABS(eval[2])) {	if (ABS(eval[0]) > ABS(eval[2])) {
320	/* normalize evec+0 but don't move it */	/* normalize evec[0] but don't move it */
321	VEC_NORM(evec+0, len);	evec[0] = normalize3(evec[0]);
322	dot = VEC_DOT(evec+0, evec+3); VEC_SCL_SUB(evec+3, dot, evec+0);	evec[1] -= scale3(dot3(evec[0], evec[1]), evec[0]);
323	VEC_NORM(evec+3, len);	evec[1] = normalize3(evec[1]);
324	dot = VEC_DOT(evec+0, evec+6); VEC_SCL_SUB(evec+6, dot, evec+0);	evec[2] -= scale3(dot3(evec[0], evec[2]), evec[0]);
325	dot = VEC_DOT(evec+3, evec+6); VEC_SCL_SUB(evec+6, dot, evec+3);	evec[2] -= scale3(dot3(evec[1], evec[2]), evec[1]);
326	VEC_NORM(evec+6, len);	evec[2] = normalize3(evec[2]);
327	} else {	}
328	/* normalize evec+6 but don't move it */	else {
329	VEC_NORM(evec+6, len);	/* normalize evec[2] but don't move it */
330	dot = VEC_DOT(evec+6, evec+3); VEC_SCL_SUB(evec+3, dot, evec+6);	evec[2] = normalize3(evec[2]);
331	VEC_NORM(evec+3, len);	evec[1] -= scale3(dot3(evec[2], evec[1]), evec[2]);
332	dot = VEC_DOT(evec+3, evec+0); VEC_SCL_SUB(evec+0, dot, evec+3);	evec[1] = normalize3(evec[1]);
333	dot = VEC_DOT(evec+6, evec+0); VEC_SCL_SUB(evec+0, dot, evec+6);	evec[0] -= scale3(dot3(evec[1], evec[0]), evec[1]);
334	VEC_NORM(evec+0, len);	evec[0] -= scale3(dot3(evec[2], evec[0]), evec[2]);
335		evec[0] = normalize3(evec[0]);
336	}	}
337
338	/* to be nice, make it right-handed */	/* to be nice, make it right-handed */
339		crs = cross(evec[0], evec[1]);
340	VEC_CROSS(crs, evec+0, evec+3);	VEC_CROSS(crs, evec+0, evec+3);
341	if (0 > VEC_DOT(crs, evec+6)) {	if (0 > dot3(crs, evec[2])) {
342	VEC_SCL(evec+6, -1);	evec[2] = -evec[2];
343	}	}
344
345	/* multiply back by eigenvalue L2 norm */	/* multiply back by eigenvalue L2 norm */

---

Colored Diff Long Colored Diff Unidiff Context Diff Side by Side

Legend: Removed from v.1183   changed lines   Added in v.1184

---

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