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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1115 - (download) (as text) (annotate)
Thu May 5 04:42:18 2011 UTC (10 years ago) by jhr
File size: 11501 byte(s)
  More merging of pure-cfg back into trunk
/*
  Teem: Tools to process and visualize scientific data and images              
  Copyright (C) 2011, 2010, 2009, University of Chicago
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public License
  (LGPL) as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  The terms of redistributing and/or modifying this software also
  include exceptions to the LGPL that facilitate static linking.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with this library; if not, write to Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/


#include "../ten.h"

char *info = ("tests tenEigensolve_d and new stand-alone function.");

#define ROOT_TRIPLE 2           /* ell_cubic_root_triple */
#define ROOT_SINGLE_DOUBLE 3    /* ell_cubic_root_single_double */
#define ROOT_THREE 4            /* ell_cubic_root_three */

#define ABS(a) (((a) > 0.0f ? (a) : -(a)))
#define VEC_SET(v, a, b, c) \
  ((v)[0] = (a), (v)[1] = (b), (v)[2] = (c))
#define VEC_DOT(v1, v2) \
  ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
#define VEC_CROSS(v3, v1, v2) \
  ((v3)[0] = (v1)[1]*(v2)[2] - (v1)[2]*(v2)[1], \
   (v3)[1] = (v1)[2]*(v2)[0] - (v1)[0]*(v2)[2], \
   (v3)[2] = (v1)[0]*(v2)[1] - (v1)[1]*(v2)[0])
#define VEC_ADD(v1, v2)  \
  ((v1)[0] += (v2)[0], \
   (v1)[1] += (v2)[1], \
   (v1)[2] += (v2)[2])
#define VEC_SUB(v1, v2)  \
  ((v1)[0] -= (v2)[0], \
   (v1)[1] -= (v2)[1], \
   (v1)[2] -= (v2)[2])
#define VEC_SCL(v1, s)  \
  ((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])

/*
** All the three given vectors span only a 2D space, and this finds
** the normal to that plane.  Simply sums up all the pair-wise
** cross-products to get a good estimate.  Trick is getting the cross
** products to line up before summing.
*/
void
nullspace1(double ret[3], 
	   const double r0[3], const double r1[3], const double r2[3]) {
  double crs[3];
  
  /* ret = r0 x r1 */
  VEC_CROSS(ret, r0, r1);
  /* crs = r1 x r2 */
  VEC_CROSS(crs, r1, r2);
  /* ret += crs or ret -= crs; whichever makes ret longer */
  if (VEC_DOT(ret, crs) > 0) {
    VEC_ADD(ret, crs);
  } else {
    VEC_SUB(ret, crs);
  }
  /* crs = r0 x r2 */
  VEC_CROSS(crs, r0, r2);
  /* ret += crs or ret -= crs; whichever makes ret longer */
  if (VEC_DOT(ret, crs) > 0) {
    VEC_ADD(ret, crs);
  } else {
    VEC_SUB(ret, crs);
  }

  return;
}

/*
** All vectors are in the same 1D space, we have to find two 
** mutually vectors perpendicular to that span
*/
void
nullspace2(double reta[3], double retb[3],
	   const double r0[3], const double r1[3], const double r2[3]) {
  double sqr[3], sum[3];
  int idx;

  VEC_COPY(sum, r0);
  if (VEC_DOT(sum, r1) > 0) {
    VEC_ADD(sum, r1);
  } else {
    VEC_SUB(sum, r1);
  }
  if (VEC_DOT(sum, r2) > 0) {
    VEC_ADD(sum, r2);
  } else {
    VEC_SUB(sum, r2);
  }
  /* 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];
  idx = 0;
  if (sqr[0] < sqr[1])
    idx = 1;
  if (sqr[idx] < sqr[2])
    idx = 2;
  /* reta will be perpendicular to sum */
  if (0 == idx) {
    VEC_SET(reta, sum[1] - sum[2], -sum[0], sum[0]);
  } else if (1 == idx) {
    VEC_SET(reta, -sum[1], sum[0] - sum[2], sum[1]);
  } else {
    VEC_SET(reta, -sum[2], sum[2], sum[0] - sum[1]);
  }
  /* and now retb will be perpendicular to both reta and sum */
  VEC_CROSS(retb, reta, sum);
  return;
}

/*
** Eigensolver for symmetric 3x3 matrix:
**
**  M00  M01  M02
**  M01  M11  M12
**  M02  M12  M22
**
** Must be passed eval[3] vector, and will compute eigenvalues
** only if evec[9] is non-NULL.  Computed eigenvectors are at evec+0,
** evec+3, and evec+6.
**
** Return value indicates something about the eigenvalue solution to
** the cubic characteristic equation; see ROOT_ #defines above
**
** 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
  }

** 
** HEY: the numerical precision issues here are very subtle, and
** merit some more scrutiny.  With evals (1.000001, 1, 1), for example,
** whether it comes back as a single/double root, vs three distinct roots,
** is determines by the comparison between "D" and "epsilon", and the 
** setting of epsilon seems pretty arbitrary at this point...
** 
*/
int
evals(double eval[3],
      const double _M00, const double _M01, const double _M02, 
      const double _M11, const double _M12, 
      const double _M22) {

#include "teigen-evals-A.c"

#include "teigen-evals-B.c"

  return roots;
}

int
evals_evecs(double eval[3], double evec[9],
	    const double _M00, const double _M01, const double _M02, 
	    const double _M11, const double _M12, 
	    const double _M22) {
  double r0[3], r1[3], r2[3], crs[3], len, dot;

#include "teigen-evals-A.c"

  /* r0, r1, r2 are the vectors we manipulate to 
     find the nullspaces of M - lambda*I */
  VEC_SET(r0, 0.0, M01, M02);
  VEC_SET(r1, M01, 0.0, M12);
  VEC_SET(r2, M02, M12, 0.0);
  if (ROOT_THREE == roots) {
    r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
    nullspace1(evec+0, r0, r1, r2);
    r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
    nullspace1(evec+3, r0, r1, r2);
    r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
    nullspace1(evec+6, r0, r1, r2);
  } else if (ROOT_SINGLE_DOUBLE == roots) {
    if (eval[1] == eval[2]) {
      /* one big (eval[0]) , two small (eval[1,2]) */
      r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
      nullspace1(evec+0, r0, r1, r2);
      r0[0] = M00 - eval[1]; r1[1] = M11 - eval[1]; r2[2] = M22 - eval[1];
      nullspace2(evec+3, evec+6, r0, r1, r2);
    }
    else {
      /* two big (eval[0,1]), one small (eval[2]) */
      r0[0] = M00 - eval[0]; r1[1] = M11 - eval[0]; r2[2] = M22 - eval[0];
      nullspace2(evec+0, evec+3, r0, r1, r2);
      r0[0] = M00 - eval[2]; r1[1] = M11 - eval[2]; r2[2] = M22 - eval[2];
      nullspace1(evec+6, r0, r1, r2);
    }
  } else {
    /* ROOT_TRIPLE == roots; use any basis for eigenvectors */
    VEC_SET(evec+0, 1, 0, 0);
    VEC_SET(evec+3, 0, 1, 0);
    VEC_SET(evec+6, 0, 0, 1);
  }
  /* we always make sure its really orthonormal; keeping fixed the
     eigenvector associated with the largest-magnitude eigenvalue */
  if (ABS(eval[0]) > ABS(eval[2])) {
    /* normalize evec+0 but don't move it */
    VEC_NORM(evec+0, len);
    dot = VEC_DOT(evec+0, evec+3); VEC_SCL_SUB(evec+3, dot, evec+0);
    VEC_NORM(evec+3, len);
    dot = VEC_DOT(evec+0, evec+6); VEC_SCL_SUB(evec+6, dot, evec+0);
    dot = VEC_DOT(evec+3, evec+6); VEC_SCL_SUB(evec+6, dot, evec+3);
    VEC_NORM(evec+6, len);
  } else {
    /* normalize evec+6 but don't move it */
    VEC_NORM(evec+6, len);
    dot = VEC_DOT(evec+6, evec+3); VEC_SCL_SUB(evec+3, dot, evec+6);
    VEC_NORM(evec+3, len);
    dot = VEC_DOT(evec+3, evec+0); VEC_SCL_SUB(evec+0, dot, evec+3);
    dot = VEC_DOT(evec+6, evec+0); VEC_SCL_SUB(evec+0, dot, evec+6);
    VEC_NORM(evec+0, len);
  }
    
  /* to be nice, make it right-handed */
  VEC_CROSS(crs, evec+0, evec+3);
  if (0 > VEC_DOT(crs, evec+6)) {
    VEC_SCL(evec+6, -1);
  }

#include "teigen-evals-B.c"

  return roots;
}


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;
}

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