Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

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;
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 */

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