[Toon-members] TooN SymEigen.h

[Ethan Eade]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Ethan Eade <ethaneade>  08/02/22 16:01:53

Modified files:
        .              : SymEigen.h 

Log message:
        Added specialization for SymEigen<2> which does not use LAPACK, but 
yields
        the same results.

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/SymEigen.h?cvsroot=toon&r1=1.10&r2=1.11

Patches:
Index: SymEigen.h
===================================================================
RCS file: /cvsroot/toon/TooN/SymEigen.h,v
retrieving revision 1.10
retrieving revision 1.11
diff -u -b -r1.10 -r1.11
--- SymEigen.h  16 Mar 2007 19:29:18 -0000      1.10
+++ SymEigen.h  22 Feb 2008 16:01:53 -0000      1.11
@@ -21,6 +21,7 @@
 #define __SYMEIGEN_H
 
 #include <iostream>
+#include <cassert>
 #include <TooN/lapack.h>
 
 #include <TooN/TooN.h>
@@ -31,20 +32,11 @@
 
 static const double symeigen_condition_no=1e9;
 
-template <int Size=-1>
-class SymEigen {
-public:
-  inline SymEigen(){}
-
-
-  template<class Accessor>
-  inline SymEigen(const FixedMatrix<Size,Size,Accessor>& m){
-    compute(m);
-  }
+    template <int Size> struct ComputeSymEigen {
 
   template<class Accessor>
-  inline void compute(const FixedMatrix<Size,Size,Accessor>& m){
-    my_evectors = m;
+       static inline void compute(const FixedMatrix<Size,Size,Accessor>& m, 
Matrix<Size,Size,RowMajor>& evectors, Vector<Size>& evalues) {
+           evectors = m;
     int N = Size;
     int lda = Size;
     int info;
@@ -52,13 +44,13 @@
     double size;
 
     // find out how much space fortran needs
-    
dsyev_("V","U",&N,my_evectors.get_data_ptr(),&lda,my_evalues.get_data_ptr(),
+           
dsyev_((char*)"V",(char*)"U",&N,evectors.get_data_ptr(),&lda,evalues.get_data_ptr(),
           &size,&lwork,&info);
     lwork = int(size);
     double* WORK = new double[lwork];
 
     // now compute the decomposition
-    
dsyev_("V","U",&N,my_evectors.get_data_ptr(),&lda,my_evalues.get_data_ptr(),
+           
dsyev_((char*)"V",(char*)"U",&N,evectors.get_data_ptr(),&lda,evalues.get_data_ptr(),
           WORK,&lwork,&info);
     delete [] WORK;
     if(info!=0){
@@ -67,6 +59,50 @@
                << "M = " << m << std::endl;
     }
   }
+    };
+
+    template <> struct ComputeSymEigen<2> {
+       
+       template<class Accessor>
+       static inline void compute(const FixedMatrix<2,2,Accessor>& m, 
Matrix<2,2,RowMajor>& eig, Vector<2>& ev) {
+           double trace = m[0][0] + m[1][1];
+           double det = m[0][0]*m[1][1] - m[0][1]*m[1][0];
+           double disc = trace*trace - 4 * det;
+           assert(disc>=0);
+           double root_disc = sqrt(disc);
+           ev[0] = 0.5 * (trace - root_disc);
+           ev[1] = 0.5 * (trace + root_disc);
+           double a = m[0][0] - ev[0];
+           double b = m[0][1];
+           double magsq = a*a + b*b;
+           if (magsq == 0) {
+               eig[0][0] = 1.0;
+               eig[0][1] = 0;
+           } else {
+               eig[0][0] = -b;
+               eig[0][1] = a;
+               eig[0] *= 1.0/sqrt(magsq);
+           }
+           eig[1][0] = -eig[0][1];
+           eig[1][1] = eig[0][0];
+       }       
+    };
+    
+template <int Size=-1>
+class SymEigen {
+public:
+  inline SymEigen(){}
+
+
+  template<class Accessor>
+  inline SymEigen(const FixedMatrix<Size,Size,Accessor>& m){
+    compute(m);
+  }
+
+  template<class Accessor>
+  inline void compute(const FixedMatrix<Size,Size,Accessor>& m){
+      ComputeSymEigen<Size>::compute(m, my_evectors, my_evalues);
+  }
   
   template <class Accessor>
   Vector<Size> backsub(const FixedVector<Size,Accessor>& rhs){