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

SCM Repository

[diderot] Diff of /branches/pure-cfg/src/lib/c-target/eigen2x2.c
ViewVC logotype

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

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

branches/pure-cfg/src/lib/c-target/eigen2D.c revision 1543, Mon Oct 17 19:53:03 2011 UTC branches/pure-cfg/src/lib/c-target/eigen2x2.c revision 1544, Tue Oct 18 13:07:58 2011 UTC
# Line 37  Line 37 
37  #define ROOT_TWO    2  #define ROOT_TWO    2
38    
39  #ifdef DIDEROT_DOUBLE_PRECISION  #ifdef DIDEROT_DOUBLE_PRECISION
40    #  define EPSILON       1.0E-12
41  # define SQRT(x)        sqrt(x)  # define SQRT(x)        sqrt(x)
42  # define VEC(a,b)       vec2d(a,b)  # define VEC(a,b)       vec2d(a,b)
43  # define DOT(u,v)       dot2d(u,v)  # define DOT(u,v)       dot2d(u,v)
44  # define SCALE(s,v)     scale2d(s,v)  # define SCALE(s,v)     scale2d(s,v)
45  # define NORMALIZE(v)   normalize2d(v)  # define NORMALIZE(v)   normalize2d(v)
46  #else  #else
47    #  define EPSILON       1.0E-12         // FIXME: probably want a different value of epsilon for single-precision
48  # define SQRT(x)        sqrtf(x)  # define SQRT(x)        sqrtf(x)
49  # define VEC(a,b)       vec2f(a,b)  # define VEC(a,b)       vec2f(a,b)
50  # define DOT(u,v)       dot2f(u,v)  # define DOT(u,v)       dot2f(u,v)
51  # define SCALE(s,v)     scale2f(s,v)  # define SCALE(s,v)     scale2f(s,v)
52  # define NORMALIZE(v)   normalize2f(v)  # define NORMALIZE(v)   normalize2f(v)
53  #endif  #endif
54  #define SUB(v,i)       (((Diderot_union3_t)(v)).r[i])  #define SUB(v,i)       (((Diderot_union2_t)(v)).r[i])
55    
56  // OpenCL 1.0 does not specify the C representation of the host types  // OpenCL 1.0 does not specify the C representation of the host types
57  // for OpenCL vector types (e.g., cl_float4), so we have to handle this  // for OpenCL vector types (e.g., cl_float4), so we have to handle this
# Line 74  Line 76 
76  **  **
77  ** HEY: the numerical precision issues here merit some more scrutiny.  ** HEY: the numerical precision issues here merit some more scrutiny.
78  */  */
79  int Diderot_evals2D (  int Diderot_evals2x2 (
80          Diderot_real_t eval[2],          Diderot_real_t eval[2],
81          const Diderot_real_t _M00, const Diderot_real_t _M01,          const Diderot_real_t _M00, const Diderot_real_t _M01,
82          const Diderot_real_t _M11)          const Diderot_real_t _M11)
83  {  {
84      Diderot_real_t mean, Q, D, M00, M01, M11;      Diderot_real_t mean, Q, D, M00, M01, M11;
 // FIXME: probably want a different value of epsilon for single-precision  
     Diderot_real_t epsilon = 1.0E-12;  
85      int roots;      int roots;
86    
87    /* copy the given matrix elements */    /* copy the given matrix elements */
# Line 99  Line 99 
99    
100      Q = M00 - M11;      Q = M00 - M11;
101      D = 4.0*M01*M01 + Q*Q;      D = 4.0*M01*M01 + Q*Q;
102      if (D > epsilon) {      if (D > EPSILON) {
103        /* two distinct roots */        /* two distinct roots */
104        Diderot_real_t vv;        Diderot_real_t vv;
105        vv = SQRT(D)/2.0;        vv = SQRT(D)/2.0;
106        eval[0] = vv;        eval[0] = vv;
107        eval[1] = -vv;        eval[1] = -vv;
108        roots = ROOT_TWO;        roots = ROOT_TWO;
109      } else {      }
110        else {
111        /* double root */        /* double root */
112        eval[0] = eval[1] = 0.0;        eval[0] = eval[1] = 0.0;
113        roots = ROOT_DOUBLE;        roots = ROOT_DOUBLE;
# Line 119  Line 120 
120      return roots;      return roots;
121  }  }
122    
123  int Diderot_evals_evecs2D (  int Diderot_evecs2x2 (
124          Diderot_real_t eval[2], Diderot_vec2_t evec[2],          Diderot_real_t eval[2], Diderot_vec2_t evec[2],
125          const Diderot_real_t _M00, const Diderot_real_t _M01,          const Diderot_real_t _M00, const Diderot_real_t _M01,
126          const Diderot_real_t _M11)          const Diderot_real_t _M11)
127  {  {
128      Diderot_real_t mean, Q, D, M00, M01, M11;      Diderot_real_t mean, Q, D, M00, M01, M11;
 // FIXME: probably want a different value of epsilon for single-precision  
     Diderot_real_t epsilon = 1.0E-12;  
129      int roots;      int roots;
130    
131    /* copy the given matrix elements */    /* copy the given matrix elements */
# Line 144  Line 143 
143    
144      Q = M00 - M11;      Q = M00 - M11;
145      D = 4.0*M01*M01 + Q*Q;      D = 4.0*M01*M01 + Q*Q;
146      if (D > epsilon) {      if (D > EPSILON) {
147        /* two distinct roots */        /* two distinct roots */
148        Diderot_real_t vv;        Diderot_real_t vv;
149        Diderot_vec2_t r1, r2;        Diderot_vec2_t r1, r2;
# Line 160  Line 159 
159        r1 = VEC(M00 - vv, M01);        r1 = VEC(M00 - vv, M01);
160        r2 = VEC(M01, M11 - vv);        r2 = VEC(M01, M11 - vv);
161        if (DOT(r1,r2) > 0.0) {        if (DOT(r1,r2) > 0.0) {
162          evec[1] = VEC(r1[0]+r2[0], r1[1]+r2[1]);              evec[1] = VEC(SUB(r1,0)+SUB(r2,0), SUB(r1,1)+SUB(r2,1));
163        } else {          }
164          evec[1] = VEC(r1[0]-r2[0], r1[1]-r2[1]);          else {
165                evec[1] = VEC(SUB(r1,0)-SUB(r2,0), SUB(r1,1)-SUB(r2,1));
166        }        }
167        evec[1] = NORMALIZE(evec[1]);        evec[1] = NORMALIZE(evec[1]);
168        evec[0] = VEC(evec[1][1], -evec[1][0]);          evec[0] = VEC(SUB(evec[1],1), -SUB(evec[1],0));
169        evec[0] = NORMALIZE(evec[0]);        evec[0] = NORMALIZE(evec[0]);
170        roots = ROOT_TWO;        roots = ROOT_TWO;
171      } else {      }
172        else {
173        /* double root */        /* double root */
174        eval[0] = eval[1] = 0.0;        eval[0] = eval[1] = 0.0;
175        /* use any basis for eigenvectors */        /* use any basis for eigenvectors */

Legend:
Removed from v.1543  
changed lines
  Added in v.1544

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