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

SCM Repository

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

Diff of /trunk/src/compiler/cl-target/fragments/eigen3x3.in

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

revision 1874, Wed May 16 13:30:16 2012 UTC revision 1875, Wed May 16 23:07:38 2012 UTC
# Line 155  Line 155 
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);
161          roots = ROOT_THREE;          roots = ROOT_THREE;
162      }      }
163    /* else D is near enough to zero */    /* else D is near enough to zero */
# Line 166  Line 166 
166        /* one double root and one single root */        /* one double root and one single root */
167          U = cbrt(R); /* cube root function */          U = cbrt(R); /* cube root function */
168          if (U > 0) {          if (U > 0) {
169              eval->s0 = 2*U;              (*eval)s0 = 2*U;
170              eval->s1 = -U;              (*eval)s1 = -U;
171              eval->s2 = -U;              (*eval)s2 = -U;
172          }          }
173          else {          else {
174              eval->s0 = -U;              (*eval)s0 = -U;
175              eval->s1 = -U;              (*eval)s1 = -U;
176              eval->s2 = 2*U;              (*eval)s2 = 2*U;
177          }          }
178          roots = ROOT_SINGLE_DOUBLE;          roots = ROOT_SINGLE_DOUBLE;
179      }      }
180      else {      else {
181        /* a triple root! */        /* a triple root! */
182          eval->s0 = eval->s1 = eval->s2 = 0.0;          (*eval)s0 = (*eval)s1 = (*eval)s2 = 0.0;
183          roots = ROOT_TRIPLE;          roots = ROOT_TRIPLE;
184      }      }
185    
186    /* multiply back by eigenvalue L2 norm */    /* multiply back by eigenvalue L2 norm */
187      eval->s0 /= rnorm;      (*eval)s0 /= rnorm;
188      eval->s1 /= rnorm;      (*eval)s1 /= rnorm;
189      eval->s2 /= rnorm;      (*eval)s2 /= rnorm;
190    
191    /* add back in the eigenvalue mean */    /* add back in the eigenvalue mean */
192      eval->s0 += mean;      (*eval)s0 += mean;
193      eval->s1 += mean;      (*eval)s1 += mean;
194      eval->s2 += mean;      (*eval)s2 += mean;
195    
196      return roots;      return roots;
197  }  }
# Line 254  Line 254 
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);
260          roots = ROOT_THREE;          roots = ROOT_THREE;
261      }      }
262    /* else D is near enough to zero */    /* else D is near enough to zero */
# Line 265  Line 265 
265        /* one double root and one single root */        /* one double root and one single root */
266          U = cbrt(R); /* cube root function */          U = cbrt(R); /* cube root function */
267          if (U > 0) {          if (U > 0) {
268              eval->s0 = 2*U;              (*eval)s0 = 2*U;
269              eval->s1 = -U;              (*eval)s1 = -U;
270              eval->s2 = -U;              (*eval)s2 = -U;
271          } else {          } else {
272              eval->s0 = -U;              (*eval)s0 = -U;
273              eval->s1 = -U;              (*eval)s1 = -U;
274              eval->s2 = 2*U;              (*eval)s2 = 2*U;
275          }          }
276          roots = ROOT_SINGLE_DOUBLE;          roots = ROOT_SINGLE_DOUBLE;
277      }      }
278      else {      else {
279        /* a triple root! */        /* a triple root! */
280          eval->s0 = eval->s1 = eval->s2 = 0.0;          (*eval)s0 = (*eval)s1 = (*eval)s2 = 0.0;
281          roots = ROOT_TRIPLE;          roots = ROOT_TRIPLE;
282      }      }
283    
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->s0, M01, M02),                  VEC3(M00 - (*eval)s0, M01, M02),
288                  VEC3(M01, M11 - eval->s0, M12),                  VEC3(M01, M11 - (*eval)s0, M12),
289                  VEC3(M02, M12, M22 - eval->s0));                  VEC3(M02, M12, M22 - (*eval)s0));
290          evecs[1] = nullspace1 (          evecs[1] = nullspace1 (
291                  VEC3(M00 - eval->s1, M01, M02),                  VEC3(M00 - (*eval)s1, M01, M02),
292                  VEC3(M01, M11 - eval->s1, M12),                  VEC3(M01, M11 - (*eval)s1, M12),
293                  VEC3(M02, M12, M22 - eval->s1));                  VEC3(M02, M12, M22 - (*eval)s1));
294          evecs[2] = nullspace1 (          evecs[2] = nullspace1 (
295                  VEC3(M00 - eval->s2, M01, M02),                  VEC3(M00 - (*eval)s2, M01, M02),
296                  VEC3(M01, M11 - eval->s2, M12),                  VEC3(M01, M11 - (*eval)s2, M12),
297                  VEC3(M02, M12, M22 - eval->s2));                  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->s0, M01, M02),                      VEC3(M00 - (*eval)s0, M01, M02),
304                      VEC3(M01, M11 - eval->s0, M12),                      VEC3(M01, M11 - (*eval)s0, M12),
305                      VEC3(M02, M12, M22 - eval->s0));                      VEC3(M02, M12, M22 - (*eval)s0));
306              nullspace2 (evecs+1,              nullspace2 (evecs+1,
307                      VEC3(M00 - eval->s1, M01, M02),                      VEC3(M00 - (*eval)s1, M01, M02),
308                      VEC3(M01, M11 - eval->s1, M12),                      VEC3(M01, M11 - (*eval)s1, M12),
309                      VEC3(M02, M12, M22 - eval->s1));                      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->s0, M01, M02),                      VEC3(M00 - (*eval)s0, M01, M02),
315                      VEC3(M01, M11 - eval->s0, M12),                      VEC3(M01, M11 - (*eval)s0, M12),
316                      VEC3(M02, M12, M22 - eval->s0));                      VEC3(M02, M12, M22 - (*eval)s0));
317              evecs[2] = nullspace1 (              evecs[2] = nullspace1 (
318                      VEC3(M00 - eval->s2, M01, M02),                      VEC3(M00 - (*eval)s2, M01, M02),
319                      VEC3(M01, M11 - eval->s2, M12),                      VEC3(M01, M11 - (*eval)s2, M12),
320                      VEC3(M02, M12, M22 - eval->s2));                      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 (fabs(eval->s0) > fabs(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];
# Line 350  Line 350 
350         The next iteration of the code will do that */         The next iteration of the code will do that */
351    
352    /* multiply back by eigenvalue L2 norm */    /* multiply back by eigenvalue L2 norm */
353      eval->s0 /= rnorm;      (*eval)s0 /= rnorm;
354      eval->s1 /= rnorm;      (*eval)s1 /= rnorm;
355      eval->s2 /= rnorm;      (*eval)s2 /= rnorm;
356    
357    /* add back in the eigenvalue mean */    /* add back in the eigenvalue mean */
358      eval->s0 += mean;      (*eval)s0 += mean;
359      eval->s1 += mean;      (*eval)s1 += mean;
360      eval->s2 += mean;      (*eval)s2 += mean;
361    
362      return roots;      return roots;
363  }  }

Legend:
Removed from v.1874  
changed lines
  Added in v.1875

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