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

SCM Repository

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

Diff of /branches/pure-cfg/src/compiler/cl-target/fragments/eigen3x3.in

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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[0] = nullspace1 (          evecs[0] = nullspace1 (
287                  VEC3(M00 - eval[0], M01, M02),                  VEC3(M00 - eval->s0, M01, M02),
288                  VEC3(M01, M11 - eval[0], M12),                  VEC3(M01, M11 - eval->s0, M12),
289                  VEC3(M02, M12, M22 - eval[0]));                  VEC3(M02, M12, M22 - eval->s0));
290          evecs[1] = nullspace1 (          evecs[1] = nullspace1 (
291                  VEC3(M00 - eval[1], M01, M02),                  VEC3(M00 - eval->s1, M01, M02),
292                  VEC3(M01, M11 - eval[1], M12),                  VEC3(M01, M11 - eval->s1, M12),
293                  VEC3(M02, M12, M22 - eval[1]));                  VEC3(M02, M12, M22 - eval->s1));
294          evecs[2] = nullspace1 (          evecs[2] = nullspace1 (
295                  VEC3(M00 - eval[2], M01, M02),                  VEC3(M00 - eval->s2, M01, M02),
296                  VEC3(M01, M11 - eval[2], M12),                  VEC3(M01, M11 - eval->s2, M12),
297                  VEC3(M02, M12, M22 - eval[2]));                  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[0]) , two small (eval[1,2]) */            /* one big (eval[0]) , two small (eval[1,2]) */
302              evecs[0] = nullspace1 (              evecs[0] = nullspace1 (
303                      VEC3(M00 - eval[0], M01, M02),                      VEC3(M00 - eval->s0, M01, M02),
304                      VEC3(M01, M11 - eval[0], M12),                      VEC3(M01, M11 - eval->s0, M12),
305                      VEC3(M02, M12, M22 - eval[0]));                      VEC3(M02, M12, M22 - eval->s0));
306              nullspace2 (evecs+1,              nullspace2 (evecs+1,
307                      VEC3(M00 - eval[1], M01, M02),                      VEC3(M00 - eval->s1, M01, M02),
308                      VEC3(M01, M11 - eval[1], M12),                      VEC3(M01, M11 - eval->s1, M12),
309                      VEC3(M02, M12, M22 - eval[1]));                      VEC3(M02, M12, M22 - eval->s1));
310          }          }
311          else {          else {
312            /* two big (eval[0,1]), one small (eval[2]) */            /* two big (eval[0,1]), one small (eval[2]) */
313              nullspace2 (evecs,              nullspace2 (evecs,
314                      VEC3(M00 - eval[0], M01, M02),                      VEC3(M00 - eval->s0, M01, M02),
315                      VEC3(M01, M11 - eval[0], M12),                      VEC3(M01, M11 - eval->s0, M12),
316                      VEC3(M02, M12, M22 - eval[0]));                      VEC3(M02, M12, M22 - eval->s0));
317              evecs[2] = nullspace1 (              evecs[2] = nullspace1 (
318                      VEC3(M00 - eval[2], M01, M02),                      VEC3(M00 - eval->s2, M01, M02),
319                      VEC3(M01, M11 - eval[2], M12),                      VEC3(M01, M11 - eval->s2, M12),
320                      VEC3(M02, M12, M22 - eval[2]));                      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[0] but don't move it */        /* normalize evec[0] but don't move it */
334          evecs[0] = normalize(evecs[0]);          evecs[0] = normalize(evecs[0]);
335          evecs[1] -= dot(evecs[1], evecs[0]) * evecs[0];          evecs[1] -= dot(evecs[1], evecs[0]) * evecs[0];

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

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