Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] Diff of /branches/pure-cfg/src/compiler/cl-target/fragments/eigen3x3.in
 [diderot] / branches / pure-cfg / src / compiler / cl-target / fragments / eigen3x3.in # Diff of /branches/pure-cfg/src/compiler/cl-target/fragments/eigen3x3.in

revision 1663, Tue Nov 29 09:48:33 2011 UTC revision 1664, Tue Nov 29 10:43:54 2011 UTC
# Line 151  Line 151
151      if (D > EPSILON) {      if (D > EPSILON) {
152        /* three distinct roots- this is the most common case */        /* three distinct roots- this is the most common case */
153          double mm, ss, cc;          double mm, ss, cc;
154          theta = ATAN2(sqrt(D), R)/3.0;          theta = atan2(sqrt(D), R)/3.0;
155          mm = sqrt(Q);          mm = sqrt(Q);
156          ss = SIN(theta);          ss = sin(theta);
157          cc = COS(theta);          cc = cos(theta);
158          eval->s0 = 2*mm*cc;          eval->s0 = 2*mm*cc;
159          eval->s1 = mm*(-cc + M_SQRT3*ss);          eval->s1 = mm*(-cc + M_SQRT3*ss);
160          eval->s2 = mm*(-cc - M_SQRT3*ss);          eval->s2 = mm*(-cc - M_SQRT3*ss);
# Line 250  Line 250
250      if (D > epsilon) {      if (D > epsilon) {
251        /* three distinct roots- this is the most common case */        /* three distinct roots- this is the most common case */
252          double mm, ss, cc;          double mm, ss, cc;
253          theta = ATAN2(sqrt(D), R)/3.0;          theta = atan2(sqrt(D), R)/3.0;
254          mm = sqrt(Q);          mm = sqrt(Q);
255          ss = SIN(theta);          ss = sin(theta);
256          cc = COS(theta);          cc = cos(theta);
257          eval->s0 = 2*mm*cc;          eval->s0 = 2*mm*cc;
258          eval->s1 = mm*(-cc + sqrt(3.0)*ss);          eval->s1 = mm*(-cc + sqrt(3.0)*ss);
259          eval->s2 = mm*(-cc - sqrt(3.0)*ss);          eval->s2 = mm*(-cc - sqrt(3.0)*ss);
# Line 284  Line 284
284      Diderot_vec3_t ev = VEC3(eval->s0, eval->s1, eval->s2);      Diderot_vec3_t ev = VEC3(eval->s0, eval->s1, eval->s2);
285      if (ROOT_THREE == roots) {      if (ROOT_THREE == roots) {
286          evecs = nullspace1 (          evecs = nullspace1 (
287                  VEC3(M00 - eval, M01, M02),                  VEC3(M00 - eval->s0, M01, M02),
288                  VEC3(M01, M11 - eval, M12),                  VEC3(M01, M11 - eval->s0, M12),
289                  VEC3(M02, M12, M22 - eval));                  VEC3(M02, M12, M22 - eval->s0));
290          evecs = nullspace1 (          evecs = nullspace1 (
291                  VEC3(M00 - eval, M01, M02),                  VEC3(M00 - eval->s1, M01, M02),
292                  VEC3(M01, M11 - eval, M12),                  VEC3(M01, M11 - eval->s1, M12),
293                  VEC3(M02, M12, M22 - eval));                  VEC3(M02, M12, M22 - eval->s1));
294          evecs = nullspace1 (          evecs = nullspace1 (
295                  VEC3(M00 - eval, M01, M02),                  VEC3(M00 - eval->s2, M01, M02),
296                  VEC3(M01, M11 - eval, M12),                  VEC3(M01, M11 - eval->s2, M12),
297                  VEC3(M02, M12, M22 - eval));                  VEC3(M02, M12, M22 - eval->s2));
298      }      }
299      else if (ROOT_SINGLE_DOUBLE == roots) {      else if (ROOT_SINGLE_DOUBLE == roots) {
300          if (eval->s1 == eval->s2) {          if (eval->s1 == eval->s2) {
301            /* one big (eval) , two small (eval[1,2]) */            /* one big (eval) , two small (eval[1,2]) */
302              evecs = nullspace1 (              evecs = nullspace1 (
303                      VEC3(M00 - eval, M01, M02),                      VEC3(M00 - eval->s0, M01, M02),
304                      VEC3(M01, M11 - eval, M12),                      VEC3(M01, M11 - eval->s0, M12),
305                      VEC3(M02, M12, M22 - eval));                      VEC3(M02, M12, M22 - eval->s0));
306              nullspace2 (evecs+1,              nullspace2 (evecs+1,
307                      VEC3(M00 - eval, M01, M02),                      VEC3(M00 - eval->s1, M01, M02),
308                      VEC3(M01, M11 - eval, M12),                      VEC3(M01, M11 - eval->s1, M12),
309                      VEC3(M02, M12, M22 - eval));                      VEC3(M02, M12, M22 - eval->s1));
310          }          }
311          else {          else {
312            /* two big (eval[0,1]), one small (eval) */            /* two big (eval[0,1]), one small (eval) */
313              nullspace2 (evecs,              nullspace2 (evecs,
314                      VEC3(M00 - eval, M01, M02),                      VEC3(M00 - eval->s0, M01, M02),
315                      VEC3(M01, M11 - eval, M12),                      VEC3(M01, M11 - eval->s0, M12),
316                      VEC3(M02, M12, M22 - eval));                      VEC3(M02, M12, M22 - eval->s0));
317              evecs = nullspace1 (              evecs = nullspace1 (
318                      VEC3(M00 - eval, M01, M02),                      VEC3(M00 - eval->s2, M01, M02),
319                      VEC3(M01, M11 - eval, M12),                      VEC3(M01, M11 - eval->s2, M12),
320                      VEC3(M02, M12, M22 - eval));                      VEC3(M02, M12, M22 - eval->s2));
321          }          }
322      }      }
323      else {      else {
# Line 329  Line 329
329    /* we always make sure it's really orthonormal; keeping fixed the    /* we always make sure it's really orthonormal; keeping fixed the
330     * eigenvector associated with the largest-magnitude eigenvalue     * eigenvector associated with the largest-magnitude eigenvalue
331     */     */
332      if (ABS(eval->s0) > ABS(eval->s2)) {      if (fabs(eval->s0) > fabs(eval->s2)) {
333        /* normalize evec but don't move it */        /* normalize evec but don't move it */
334          evecs = normalize(evecs);          evecs = normalize(evecs);
335          evecs -= dot(evecs, evecs) * evecs;          evecs -= dot(evecs, evecs) * evecs;

Legend:
 Removed from v.1663 changed lines Added in v.1664