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/eigen.c
ViewVC logotype

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

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

revision 1106, Wed May 4 21:01:25 2011 UTC revision 1172, Tue May 10 15:57:48 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 */
# Line 259  Line 268 
268    return roots;    return roots;
269  }  }
270    
271    #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.1106  
changed lines
  Added in v.1172

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