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

SCM Repository

[diderot] Diff of /trunk/src/lib/c-target/eigen.c
ViewVC logotype

Diff of /trunk/src/lib/c-target/eigen.c

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

revision 1115, Thu May 5 04:42:18 2011 UTC revision 1232, Mon May 16 23:37:52 2011 UTC
# Line 1  Line 1 
1    /*! \file eigen.c
2     *
3     * \author Gordon Kindlmann
4     */
5    
6    /*
7     * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
8     * All rights reserved.
9     */
10    
11  /*  /*
12    Teem: Tools to process and visualize scientific data and images    Teem: Tools to process and visualize scientific data and images
13    Copyright (C) 2011, 2010, 2009, University of Chicago    Copyright (C) 2011, 2010, 2009, University of Chicago
# Line 21  Line 31 
31    51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA    51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
32  */  */
33    
34    #ifdef PORTED
35    
36  #include "../ten.h"  #include <Diderot/diderot.h>
   
 char *info = ("tests tenEigensolve_d and new stand-alone function.");  
37    
38  #define ROOT_TRIPLE 2           /* ell_cubic_root_triple */  #define ROOT_TRIPLE 2           /* ell_cubic_root_triple */
39  #define ROOT_SINGLE_DOUBLE 3    /* ell_cubic_root_single_double */  #define ROOT_SINGLE_DOUBLE 3    /* ell_cubic_root_single_double */
40  #define ROOT_THREE 4            /* ell_cubic_root_three */  #define ROOT_THREE 4            /* ell_cubic_root_three */
41    
42  #define ABS(a) (((a) > 0.0f ? (a) : -(a)))  // from http://en.wikipedia.org/wiki/Square_root_of_3
43  #define VEC_SET(v, a, b, c) \  #define M_SQRT3 1.732050807568877293527446341506
44    ((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))  
45  #define VEC_DOT(v1, v2) \  #ifdef DIDEROT_DOUBLE_PRECISION
46    ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])  # define ABS(x)         fabs(x)
47  #define VEC_CROSS(v3, v1, v2) \  # define SQRT(x)        sqrt(x)
48    ((v3)[0] = (v1)[1]*(v2)[2] - (v1)[2]*(v2)[1], \  # define SIN(x)         sin(x)
49     (v3)[1] = (v1)[2]*(v2)[0] - (v1)[0]*(v2)[2], \  # define COS(x)         cos(x)
50     (v3)[2] = (v1)[0]*(v2)[1] - (v1)[1]*(v2)[0])  # define ATAN2(x,y)     atan2(x,y)
51  #define VEC_ADD(v1, v2)  \  # define CBRT(x)        cbrt(x)
52    ((v1)[0] += (v2)[0], \  #else
53     (v1)[1] += (v2)[1], \  # define ABS(x)         fabsf(x)
54     (v1)[2] += (v2)[2])  # define SQRT(x)        sqrtf(x)
55  #define VEC_SUB(v1, v2)  \  # define SIN(x)         sinf(x)
56    ((v1)[0] -= (v2)[0], \  # define COS(x)         cosf(x)
57     (v1)[1] -= (v2)[1], \  # define ATAN2(x,y)     atan2f(x,y)
58     (v1)[2] -= (v2)[2])  # define CBRT(x)        cbrtf(x)
59  #define VEC_SCL(v1, s)  \  #endif
   ((v1)[0] *= (s),      \  
    (v1)[1] *= (s),      \  
    (v1)[2] *= (s))  
 #define VEC_LEN(v) (sqrt(VEC_DOT(v,v)))  
 #define VEC_NORM(v, len) ((len) = VEC_LEN(v), VEC_SCL(v, 1.0/len))  
 #define VEC_SCL_SUB(v1, s, v2) \  
   ((v1)[0] -= (s)*(v2)[0],      \  
    (v1)[1] -= (s)*(v2)[1], \  
    (v1)[2] -= (s)*(v2)[2])  
 #define VEC_COPY(v1, v2)  \  
   ((v1)[0] = (v2)[0], \  
    (v1)[1] = (v2)[1], \  
    (v1)[2] = (v2)[2])  
60    
61  /*  /*
62  ** All the three given vectors span only a 2D space, and this finds  ** All the three given vectors span only a 2D space, and this finds
# Line 68  Line 64 
64  ** cross-products to get a good estimate.  Trick is getting the cross  ** cross-products to get a good estimate.  Trick is getting the cross
65  ** products to line up before summing.  ** products to line up before summing.
66  */  */
67  void  STATIC_INLINE Diderot_vec3_t nullspace1 (
68  nullspace1(double ret[3],      const Diderot_vec3_t r0,
69             const double r0[3], const double r1[3], const double r2[3]) {      const Diderot_vec3_t r1,
70    double crs[3];      const Diderot_vec3_t r2)
71    {
72    /* ret = r0 x r1 */      Diderot_vec3_t ret, crs;
73    VEC_CROSS(ret, r0, r1);  
74    /* crs = r1 x r2 */      ret = cross(r0, r1);
75    VEC_CROSS(crs, r1, r2);      crs = cross(r1, r2);
76    /* ret += crs or ret -= crs; whichever makes ret longer */    /* ret += crs or ret -= crs; whichever makes ret longer */
77    if (VEC_DOT(ret, crs) > 0) {      if (dot3(ret, crs) > 0.0)
78      VEC_ADD(ret, crs);          ret += crs;
79    } else {      else
80      VEC_SUB(ret, crs);          ret -= crs;
81    }  
82    /* crs = r0 x r2 */      crs = cross(r0, r2);
   VEC_CROSS(crs, r0, r2);  
83    /* ret += crs or ret -= crs; whichever makes ret longer */    /* ret += crs or ret -= crs; whichever makes ret longer */
84    if (VEC_DOT(ret, crs) > 0) {      if (dot3(ret, crs) > 0.0)
85      VEC_ADD(ret, crs);          ret += crs;
86    } else {      else
87      VEC_SUB(ret, crs);          ret -= crs;
   }  
88    
89    return;      return ret;
90  }  }
91    
92  /*  /*
# Line 100  Line 94 
94  ** mutually vectors perpendicular to that span  ** mutually vectors perpendicular to that span
95  */  */
96  void  void
97  nullspace2(double reta[3], double retb[3],  nullspace2(Diderot_vec3_t rets[2],
98             const double r0[3], const double r1[3], const double r2[3]) {      const Diderot_vec3_t r0,
99    double sqr[3], sum[3];      const Diderot_vec3_t r1,
100        const Diderot_vec3_t r2)
101    {
102      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;
112      VEC_ADD(sum, r2);      else
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      VEC_SET(reta, 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) {
125      VEC_SET(reta, -sum[1], sum[0] - sum[2], sum[1]);        rets[0] = vec3(-sum[1], sum[0] - sum[2], sum[1]);
126    } else {    } else {
127      VEC_SET(reta, -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    VEC_CROSS(retb, reta, sum);      rets[1] = cross(rets[0], sum);
131    return;    return;
132  }  }
133    
# Line 146  Line 138 
138  **  M01  M11  M12  **  M01  M11  M12
139  **  M02  M12  M22  **  M02  M12  M22
140  **  **
141  ** Must be passed eval[3] vector, and will compute eigenvalues  ** Must be passed pointer to eval vector, and will compute eigenvalues
142  ** only if evec[9] is non-NULL.  Computed eigenvectors are at evec+0,  ** only if evec[9] is non-NULL.  Computed eigenvectors are at evec+0,
143  ** evec+3, and evec+6.  ** evec+3, and evec+6.
144  **  **
145  ** Return value indicates something about the eigenvalue solution to  ** Return value indicates something about the eigenvalue solution to
146  ** the cubic characteristic equation; see ROOT_ #defines above  ** the cubic characteristic equation; see ROOT_ #defines above
147  **  **
 ** Relies on the ABS and VEC_* macros above, as well as math functions  
 ** atan2(), sin(), cos(), sqrt(), and airCbrt(), defined as:  
   
   double  
   airCbrt(double v) {  
   #if defined(_WIN32) || defined(__STRICT_ANSI__)  
     return (v < 0.0 ? -pow(-v,1.0/3.0) : pow(v,1.0/3.0));  
   #else  
     return cbrt(v);  
   #endif  
   }  
   
 **  
148  ** HEY: the numerical precision issues here are very subtle, and  ** HEY: the numerical precision issues here are very subtle, and
149  ** merit some more scrutiny.  With evals (1.000001, 1, 1), for example,  ** merit some more scrutiny.  With evals (1.000001, 1, 1), for example,
150  ** whether it comes back as a single/double root, vs three distinct roots,  ** whether it comes back as a single/double root, vs three distinct roots,
# Line 173  Line 152 
152  ** setting of epsilon seems pretty arbitrary at this point...  ** setting of epsilon seems pretty arbitrary at this point...
153  **  **
154  */  */
155  int  int Diderot_evals (
156  evals(double eval[3],          Diderot_vec3_t eval[3],
157        const double _M00, const double _M01, const double _M02,          const Diderot_real_t _M00, const Diderot_real_t _M01, const Diderot_real_t _M02,
158        const double _M11, const double _M12,          const Diderot_real_t _M11, const Diderot_real_t _M12,
159        const double _M22) {          const Diderot_real_t _M22)
160    {
161        Diderot_real_t mean, norm, rnorm, Q, R, QQQ, D, theta, M00, M01, M02, M11, M12, M22;
162    // FIXME: probably want a different value of epsilon for single-precision
163        Diderot_real_t epsilon = 1.0E-12;
164        int roots;
165    
166  #include "teigen-evals-A.c"    /* copy the given matrix elements */
167        M00 = _M00;
168        M01 = _M01;
169        M02 = _M02;
170        M11 = _M11;
171        M12 = _M12;
172        M22 = _M22;
173    
174  #include "teigen-evals-B.c"    /*
175      ** subtract out the eigenvalue mean (will add back to evals later);
176      ** helps with numerical stability
177      */
178        mean = (M00 + M11 + M22)/3.0;
179        M00 -= mean;
180        M11 -= mean;
181        M22 -= mean;
182    
183      /*
184      ** divide out L2 norm of eigenvalues (will multiply back later);
185      ** this too seems to help with stability
186      */
187        norm = SQRT(M00*M00 + 2*M01*M01 + 2*M02*M02 +
188                    M11*M11 + 2*M12*M12 +
189                    M22*M22);
190    // QUESTION: what if norm is very close to zero?
191        rnorm = norm ? 1.0/norm : 1.0;
192        M00 *= rnorm;
193        M01 *= rnorm;
194        M02 *= rnorm;
195        M11 *= rnorm;
196        M12 *= rnorm;
197        M22 *= rnorm;
198    
199      /* this code is a mix of prior Teem code and ideas from Eberly's
200       * "Eigensystems for 3 x 3 Symmetric Matrices (Revisited)"
201       */
202        Q = (M01*M01 + M02*M02 + M12*M12 - M00*M11 - M00*M22 - M11*M22)/3.0;
203        QQQ = Q*Q*Q;
204        R = (M00*M11*M22 + M02*(2*M01*M12 - M02*M11)
205             - M00*M12*M12 - M01*M01*M22)/2.0;
206        D = QQQ - R*R;
207        if (D > epsilon) {
208          /* three distinct roots- this is the most common case */
209            double mm, ss, cc;
210            theta = ATAN2(SQRT(D), R)/3.0;
211            mm = SQRT(Q);
212            ss = SIN(theta);
213            cc = COS(theta);
214            eval[0] = 2*mm*cc;
215            eval[1] = mm*(-cc + M_SQRT3*ss);
216            eval[2] = mm*(-cc - M_SQRT3*ss);
217            roots = ROOT_THREE;
218        }
219      /* else D is near enough to zero */
220        else if (R < -epsilon || epsilon < R) {
221            double U;
222          /* one double root and one single root */
223            U = CBRT(R); /* cube root function */
224            if (U > 0) {
225                eval[0] = 2*U;
226                eval[1] = -U;
227                eval[2] = -U;
228            }
229            else {
230                eval[0] = -U;
231                eval[1] = -U;
232                eval[2] = 2*U;
233            }
234            roots = ROOT_SINGLE_DOUBLE;
235        }
236        else {
237          /* a triple root! */
238            eval[0] = eval[1] = eval[2] = 0.0;
239            roots = ROOT_TRIPLE;
240        }
241    
242      /* multiply back by eigenvalue L2 norm */
243        eval[0] /= rnorm;
244        eval[1] /= rnorm;
245        eval[2] /= rnorm;
246    
247      /* add back in the eigenvalue mean */
248        eval[0] += mean;
249        eval[1] += mean;
250        eval[2] += mean;
251    
252    return roots;    return roots;
253  }  }
254    
255  int  int evals_evecs(
256  evals_evecs(double eval[3], double evec[9],          Diderot_real_t eval[3], Diderot_vec3_t evecs[3],
257              const double _M00, const double _M01, const double _M02,          const Diderot_real_t _M00, const Diderot_real_t _M01, const Diderot_real_t _M02,
258              const double _M11, const double _M12,          const Diderot_real_t _M11, const Diderot_real_t _M12,
259              const double _M22) {          const Diderot_real_t _M22)
260    double r0[3], r1[3], r2[3], crs[3], len, dot;  {
261        Diderot_vec3_t r0, r1, r2, crs;
262        Diderot_real_t len, dot;
263    
264  #include "teigen-evals-A.c"  #include "teigen-evals-A.c"
265    
266    /* r0, r1, r2 are the vectors we manipulate to    /* r0, r1, r2 are the vectors we manipulate to
267       find the nullspaces of M - lambda*I */       find the nullspaces of M - lambda*I */
268    VEC_SET(r0, 0.0, M01, M02);      r0 = vec3(0.0, M01, M02);
269    VEC_SET(r1, M01, 0.0, M12);      r1 = vec3(M01, 0.0, M12);
270    VEC_SET(r2, M02, M12, 0.0);      r2 = vec3(M02, M12, 0.0);
271        Diderot_vec3_t ev = vec3(eval[0], eval[1], eval[2]);
272    if (ROOT_THREE == roots) {    if (ROOT_THREE == roots) {
273      r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];          evec[0] = nullspace1 (
274      nullspace1(evec+0, r0, r1, r2);                  vec3(M00 - eval[0], M01, M02);
275      r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];                  vec3(M01, M11 - eval[0], M12);
276      nullspace1(evec+3, r0, r1, r2);                  vec3(M02, M12, M22 - eval[0]));
277      r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];          evec[1] = nullspace1 (
278      nullspace1(evec+6, r0, r1, r2);                  vec3(M00 - eval[1], M01, M02);
279    } else if (ROOT_SINGLE_DOUBLE == roots) {                  vec3(M01, M11 - eval[1], M12);
280                    vec3(M02, M12, M22 - eval[1]));
281            evec[2] = nullspace1 (
282                    vec3(M00 - eval[2], M01, M02);
283                    vec3(M01, M11 - eval[2], M12);
284                    vec3(M02, M12, M22 - eval[2]));
285        }
286        else if (ROOT_SINGLE_DOUBLE == roots) {
287      if (eval[1] == eval[2]) {      if (eval[1] == eval[2]) {
288        /* one big (eval[0]) , two small (eval[1,2]) */        /* one big (eval[0]) , two small (eval[1,2]) */
289        r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];              evec[0] = nullspace1 (
290        nullspace1(evec+0, r0, r1, r2);                      vec3(M00 - eval[0], M01, M02);
291        r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];                      vec3(M01, M11 - eval[0], M12);
292        nullspace2(evec+3, evec+6, r0, r1, r2);                      vec3(M02, M12, M22 - eval[0]));
293                nullspace2 (evec+1,
294                        vec3(M00 - eval[1], M01, M02);
295                        vec3(M01, M11 - eval[1], M12);
296                        vec3(M02, M12, M22 - eval[1]));
297      }      }
298      else {      else {
299        /* two big (eval[0,1]), one small (eval[2]) */        /* two big (eval[0,1]), one small (eval[2]) */
300        r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];              nullspace2 (evec,
301        nullspace2(evec+0, evec+3, r0, r1, r2);                      vec3(M00 - eval[0], M01, M02);
302        r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];                      vec3(M01, M11 - eval[0], M12);
303        nullspace1(evec+6, r0, r1, r2);                      vec3(M02, M12, M22 - eval[0]));
304                evec[2] = nullspace1 (
305                        vec3(M00 - eval[2], M01, M02);
306                        vec3(M01, M11 - eval[2], M12);
307                        vec3(M02, M12, M22 - eval[2]));
308      }      }
309    } else {      }
310        else {
311      /* ROOT_TRIPLE == roots; use any basis for eigenvectors */      /* ROOT_TRIPLE == roots; use any basis for eigenvectors */
312      VEC_SET(evec+0, 1, 0, 0);          evec[0] = vec3(1, 0, 0);
313      VEC_SET(evec+3, 0, 1, 0);          evec[1] = vec3(0, 1, 0);
314      VEC_SET(evec+6, 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  #include "teigen-evals-B.c"    /* multiply back by eigenvalue L2 norm */
346        eval[0] /= rnorm;
347        eval[1] /= rnorm;
348        eval[2] /= rnorm;
349    
350      /* add back in the eigenvalue mean */
351        eval[0] += mean;
352        eval[1] += mean;
353        eval[2] += mean;
354    
355    return roots;    return roots;
356  }  }
357    
358    #endif /* PORTED */
 void  
 testeigen(double tt[7], double eval[3], double evec[9]) {  
   double mat[9], dot[3], cross[3];  
   unsigned int ii;  
   
   TEN_T2M(mat, tt);  
   printf("evals %g %g %g\n", eval[0], eval[1], eval[2]);  
   printf("evec0 (%g) %g %g %g\n",  
          ELL_3V_LEN(evec + 0), evec[0], evec[1], evec[2]);  
   printf("evec1 (%g) %g %g %g\n",  
          ELL_3V_LEN(evec + 3), evec[3], evec[4], evec[5]);  
   printf("evec2 (%g) %g %g %g\n",  
          ELL_3V_LEN(evec + 6), evec[6], evec[7], evec[8]);  
   printf("Mv - lv: (len) X Y Z (should be ~zeros)\n");  
   for (ii=0; ii<3; ii++) {  
     double uu[3], vv[3], dd[3];  
     ELL_3MV_MUL(uu, mat, evec + 3*ii);  
     ELL_3V_SCALE(vv, eval[ii], evec + 3*ii);  
     ELL_3V_SUB(dd, uu, vv);  
     printf("%d: (%g) %g %g %g\n", ii, ELL_3V_LEN(dd), dd[0], dd[1], dd[2]);  
   }  
   dot[0] = ELL_3V_DOT(evec + 0, evec + 3);  
   dot[1] = ELL_3V_DOT(evec + 0, evec + 6);  
   dot[2] = ELL_3V_DOT(evec + 3, evec + 6);  
   printf("pairwise dots: (%g) %g %g %g\n",  
          ELL_3V_LEN(dot), dot[0], dot[1], dot[2]);  
   ELL_3V_CROSS(cross, evec+0, evec+3);  
   printf("right-handed: %g\n", ELL_3V_DOT(evec+6, cross));  
   return;  
 }  
   
 int  
 main(int argc, char *argv[]) {  
   char *me;  
   hestOpt *hopt=NULL;  
   airArray *mop;  
   
   double _tt[6], tt[7], ss, pp[3], qq[4], rot[9], mat1[9], mat2[9], tmp,  
     evalA[3], evecA[9], evalB[3], evecB[9];  
   int roots;  
   
   mop = airMopNew();  
   me = argv[0];  
   hestOptAdd(&hopt, NULL, "m00 m01 m02 m11 m12 m22",  
              airTypeDouble, 6, 6, _tt, NULL, "symmtric matrix coeffs");  
   hestOptAdd(&hopt, "p", "vec", airTypeDouble, 3, 3, pp, "0 0 0",  
              "rotation as P vector");  
   hestOptAdd(&hopt, "s", "scl", airTypeDouble, 1, 1, &ss, "1.0",  
              "scaling");  
   hestParseOrDie(hopt, argc-1, argv+1, NULL,  
                  me, info, AIR_TRUE, AIR_TRUE, AIR_TRUE);  
   airMopAdd(mop, hopt, (airMopper)hestOptFree, airMopAlways);  
   airMopAdd(mop, hopt, (airMopper)hestParseFree, airMopAlways);  
   
   ELL_6V_COPY(tt + 1, _tt);  
   tt[0] = 1.0;  
   TEN_T_SCALE(tt, ss, tt);  
   
   ELL_4V_SET(qq, 1, pp[0], pp[1], pp[2]);  
   ELL_4V_NORM(qq, qq, tmp);  
   ell_q_to_3m_d(rot, qq);  
   printf("%s: rot\n", me);  
   printf("  %g %g %g\n", rot[0], rot[1], rot[2]);  
   printf("  %g %g %g\n", rot[3], rot[4], rot[5]);  
   printf("  %g %g %g\n", rot[6], rot[7], rot[8]);  
   
   TEN_T2M(mat1, tt);  
   ell_3m_mul_d(mat2, rot, mat1);  
   ELL_3M_TRANSPOSE_IP(rot, tmp);  
   ell_3m_mul_d(mat1, mat2, rot);  
   TEN_M2T(tt, mat1);  
   
   printf("input matrix = \n %g %g %g\n %g %g\n %g\n",  
           tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);  
   
   printf("================== tenEigensolve_d ==================\n");  
   roots = tenEigensolve_d(evalA, evecA, tt);  
   printf("%s roots\n", airEnumStr(ell_cubic_root, roots));  
   testeigen(tt, evalA, evecA);  
   
   printf("================== new eigensolve ==================\n");  
   roots = evals(evalB, tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);  
   printf("%s roots: %g %g %g\n", airEnumStr(ell_cubic_root, roots),  
          evalB[0], evalB[1], evalB[2]);  
   roots = evals_evecs(evalB, evecB,  
                       tt[1], tt[2], tt[3], tt[4], tt[5], tt[6]);  
   printf("%s roots\n", airEnumStr(ell_cubic_root, roots));  
   testeigen(tt, evalB, evecB);  
   
   airMopOkay(mop);  
   return 0;  
 }  

Legend:
Removed from v.1115  
changed lines
  Added in v.1232

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