[Toon-members] TooN SymEigen.h lapack.h test/sym.cc

[Gerhard Reitmayr]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Gerhard Reitmayr <gerhard>      09/04/14 14:56:12

Modified files:
        .              : SymEigen.h lapack.h 
Added files:
        test           : sym.cc 

Log message:
        ported SymEigen to TooN2

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/SymEigen.h?cvsroot=toon&r1=1.12&r2=1.13
http://cvs.savannah.gnu.org/viewcvs/TooN/lapack.h?cvsroot=toon&r1=1.11&r2=1.12
http://cvs.savannah.gnu.org/viewcvs/TooN/test/sym.cc?cvsroot=toon&rev=1.1

Patches:
Index: SymEigen.h
===================================================================
RCS file: /cvsroot/toon/TooN/SymEigen.h,v
retrieving revision 1.12
retrieving revision 1.13
diff -u -b -r1.12 -r1.13
--- SymEigen.h  10 Apr 2009 06:23:54 -0000      1.12
+++ SymEigen.h  14 Apr 2009 14:56:11 -0000      1.13
@@ -36,45 +36,41 @@
 
 #include <TooN/TooN.h>
 
-#ifndef TOON_NO_NAMESPACE
 namespace TooN {
-#endif 
 
 static const double symeigen_condition_no=1e9;
 
-    template <int Size> struct ComputeSymEigen {
+template <int Size> struct ComputeSymEigen {
        
-       template<class Accessor>
-       static inline void compute(const FixedMatrix<Size,Size,Accessor>& m, 
Matrix<Size,Size,RowMajor>& evectors, Vector<Size>& evalues) {
+       template<typename P, typename B>
+       static inline void compute(const Matrix<Size,Size,P, B>& m, 
Matrix<Size,Size,P> & evectors, Vector<Size, P>& evalues) {
            evectors = m;
-           int N = Size;
-           int lda = Size;
+               int N = evalues.size();
+               int lda = evalues.size();
            int info;
            int lwork=-1;
-           double size;
+               P size;
            
            // find out how much space fortran needs
-           
dsyev_((char*)"V",(char*)"U",&N,evectors.get_data_ptr(),&lda,evalues.get_data_ptr(),
-                  &size,&lwork,&info);
+               
syev_((char*)"V",(char*)"U",&N,&evectors[0][0],&lda,&evalues[0], 
&size,&lwork,&info);
            lwork = int(size);
-           double* WORK = new double[lwork];
+               Vector<Dynamic, P> WORK(lwork);
            
            // now compute the decomposition
-           
dsyev_((char*)"V",(char*)"U",&N,evectors.get_data_ptr(),&lda,evalues.get_data_ptr(),
-                  WORK,&lwork,&info);
-           delete [] WORK;
+               
syev_((char*)"V",(char*)"U",&N,&evectors[0][0],&lda,&evalues[0], 
&WORK[0],&lwork,&info);
+
            if(info!=0){
                std::cerr << "In SymEigen<"<<Size<<">: " << info 
                          << " off-diagonal elements of an intermediate 
tridiagonal form did not converge to zero." << std::endl
                          << "M = " << m << std::endl;
            }
        }       
-    };
+};
 
-    template <> struct ComputeSymEigen<2> {
+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) {
+       template<typename P, typename B>
+       static inline void compute(const Matrix<2,2,P,B>& m, Matrix<2,2,P>& 
eig, Vector<2, P>& 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;
@@ -96,79 +92,59 @@
            eig[1][0] = -eig[0][1];
            eig[1][1] = eig[0][0];
        }       
-    };
+};
     
-template <int Size=-1>
+template <int Size=Dynamic, typename Precision = double>
 class SymEigen {
 public:
   inline SymEigen(){}
 
-
-  template<class Accessor>
-  inline SymEigen(const FixedMatrix<Size,Size,Accessor>& m){
+  template<int R, int C, typename B>
+  inline SymEigen(const Matrix<R, C, Precision, B>& m) : 
my_evectors(m.num_rows(), m.num_cols()), my_evalues(m.num_rows()) {
     compute(m);
   }
 
-  template<class Accessor>
-  inline void compute(const FixedMatrix<Size,Size,Accessor>& m){
+  template<int R, int C, typename B>
+  inline void compute(const Matrix<R,C,Precision,B>& m){
+               SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
+               SizeMismatch<R, Size>::test(m.num_rows(), 
my_evectors.num_rows());
       ComputeSymEigen<Size>::compute(m, my_evectors, my_evalues);
   }
   
-  template <class Accessor>
-  Vector<Size> backsub(const FixedVector<Size,Accessor>& rhs){
-    Vector<Size> invdiag;
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
-  }
-
-  template <class Accessor>
-  Vector<> backsub(const DynamicVector<Accessor>& rhs){
-    Vector<Size> invdiag;
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
-  }
-
-
-  template <int NRHS, class Accessor>
-  Matrix<Size,NRHS> backsub(const FixedMatrix<Size,NRHS,Accessor>& rhs){
-    Vector<Size> invdiag;
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
+       template <int S, typename P, typename B>
+       Vector<Size, Precision> backsub(const Vector<S,P,B>& rhs) const {
+               return (my_evectors.T() * 
diagmult(get_inv_diag(symeigen_condition_no),(my_evectors * rhs)));
   }
 
-  template <class Accessor>
-  Matrix<> backsub(const DynamicMatrix<Accessor>& rhs){
-    Vector<Size> invdiag;
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
+       template <int R, int C, typename P, typename B>
+       Matrix<Size,C, Precision> backsub(const Matrix<R,C,P,B>& rhs) const {
+               return (my_evectors.T() * 
diagmult(get_inv_diag(symeigen_condition_no),(my_evectors * rhs)));
   }
 
-
-  Matrix<Size> get_pinv(const double condition=symeigen_condition_no){
-    Vector<Size> invdiag;
-    get_inv_diag(invdiag,condition);
-    return my_evectors.T() * diagmult(invdiag,my_evectors);
+       Matrix<Size, Size, Precision> get_pinv(const double 
condition=symeigen_condition_no) const {
+               return my_evectors.T() * 
diagmult(get_inv_diag(condition),my_evectors);
   }
 
-
-  void get_inv_diag(Vector<Size>& invdiag, double condition){
-    double max_diag = -my_evalues[0] > my_evalues[Size-1] ? 
-my_evalues[0]:my_evalues[Size-1];
-    for(int i=0; i<Size; i++){
+       Vector<Size, Precision> get_inv_diag(const double condition) const {
+               Precision max_diag = -my_evalues[0] > 
my_evalues[my_evalues.size()-1] ? 
-my_evalues[0]:my_evalues[my_evalues.size()-1];
+               Vector<Size, Precision> invdiag(my_evalues.size());
+               for(int i=0; i<my_evalues.size(); i++){
       if(fabs(my_evalues[i]) * condition > max_diag) {
        invdiag[i] = 1/my_evalues[i];
       } else {
        invdiag[i]=0;
       }
     }
+               return invdiag;
   }
   
-  inline Matrix<Size,Size,RowMajor>& get_evectors() {return my_evectors;}
-  inline const Matrix<Size,Size,RowMajor>& get_evectors() const {return 
my_evectors;}
-  inline Vector<Size>& get_evalues() {return my_evalues;}
-  inline const Vector<Size>& get_evalues() const {return my_evalues;}
+       Matrix<Size,Size,Precision>& get_evectors() {return my_evectors;}
+       const Matrix<Size,Size,Precision>& get_evectors() const {return 
my_evectors;}
+       Vector<Size, Precision>& get_evalues() {return my_evalues;}
+       const Vector<Size, Precision>& get_evalues() const {return my_evalues;}
   
   bool is_posdef() const {
-      for (int i = 0; i < Size; ++i) {
+         for (int i = 0; i < my_evalues.size(); ++i) {
           if (my_evalues[i] <= 0.0)
               return false;
       }
@@ -176,16 +152,16 @@
   }
   
   bool is_negdef() const {
-      for (int i = 0; i < Size; ++i) {
+         for (int i = 0; i < my_evalues.size(); ++i) {
           if (my_evalues[i] >= 0.0)
               return false;
       }
       return true;
   }
 
-  double get_determinant () const {
-         double det = 1.0;
-         for (int i = 0; i < Size; ++i) {
+       Precision get_determinant () const {
+         Precision det = 1.0;
+         for (int i = 0; i < my_evalues.size(); ++i) {
                  det *= my_evalues[i];
          }
          return det;
@@ -193,121 +169,11 @@
 
 private:
   // eigen vectors laid out row-wise so evectors[i] is the ith evector
-  Matrix<Size,Size,RowMajor> my_evectors;
+       Matrix<Size,Size,Precision> my_evectors;
 
-  Vector<Size> my_evalues;
+       Vector<Size, Precision> my_evalues;
 };
 
-
-template <>
-class SymEigen<> {
-public:
-
-  inline SymEigen(int size) : 
-    my_evectors(size,size),
-    my_evalues(size) {}
-
-  template<class Accessor>
-  inline SymEigen(const DynamicMatrix<Accessor>& m) : 
-    my_evectors(m.num_rows(), m.num_cols()),
-    my_evalues(m.num_rows())
-  {
-    assert(m.num_rows() == m.num_cols());
-    compute(m);
-  }
-
-  template<class Accessor>
-  inline void compute(const DynamicMatrix<Accessor>& m){
-    my_evectors = m;
-    int N = m.num_cols();
-    int lda = m.num_cols();
-    int info;
-    int lwork=-1;
-    double size;
-
-    // find out how much space fortran needs
-    
dsyev_("V","U",&N,my_evectors.get_data_ptr(),&lda,my_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(),
-          WORK,&lwork,&info);
-    delete [] WORK;
-    if(info!=0){
-      std::cerr << "info was not zero in SymEigen - it was " << info << 
std::endl;
-    }
-  }
-  
-  template <int Size, class Accessor>
-  Vector<Size> backsub(const FixedVector<Size,Accessor>& rhs){
-    Vector<> invdiag(my_evalues.size());
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
-  }
-
-  template <class Accessor>
-  Vector<> backsub(const DynamicVector<Accessor>& rhs){
-    Vector<> invdiag(my_evalues.size());
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
-  }
-
-
-  template <int NRHS, class Accessor, int Size>
-  Matrix<Size,NRHS> backsub(const FixedMatrix<Size,NRHS,Accessor>& rhs){
-    Vector<> invdiag(my_evalues.size());
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
-  }
-
-  template <class Accessor>
-  Matrix<> backsub(const DynamicMatrix<Accessor>& rhs){
-    Vector<> invdiag(my_evalues.size());
-    get_inv_diag(invdiag,symeigen_condition_no);
-    return (my_evectors.T() * diagmult(invdiag,(my_evectors * rhs)));
-  }
-
-
-
-
-  Matrix<> get_pinv(const double condition=symeigen_condition_no){
-    Vector<> invdiag(my_evalues.size());
-    get_inv_diag(invdiag,condition);
-    return my_evectors.T() * diagmult(invdiag,my_evectors);
-  }
-
-
-  void get_inv_diag(Vector<>& invdiag, double condition){
-    double max_diag = -my_evalues[0] > my_evalues[invdiag.size()-1] ? 
-my_evalues[0]:my_evalues[invdiag.size()-1];
-    for(int i=0; i<invdiag.size(); i++){
-      if(fabs(my_evalues[i]) * condition > max_diag) {
-       invdiag[i] = 1/my_evalues[i];
-      } else {
-       invdiag[i]=0;
-      }
-    }
-  }
-  
-
-  
-  inline Matrix<-1,-1,RowMajor>& get_evectors() {return my_evectors;}
-  inline const Matrix<-1,-1,RowMajor>& get_evectors() const {return 
my_evectors;}
-  inline Vector<>& get_evalues() {return my_evalues;}
-  inline const Vector<>& get_evalues() const {return my_evalues;}
-
-private:
-  // eigen vectors laid out row-wise so evectors[i] is the ith evector
-  Matrix<-1,-1,RowMajor> my_evectors;
-  Vector<> my_evalues;
-};
-
-
-#ifndef TOON_NO_NAMESPACE
 }
-#endif 
-
-
 
 #endif

Index: lapack.h
===================================================================
RCS file: /cvsroot/toon/TooN/lapack.h,v
retrieving revision 1.11
retrieving revision 1.12
diff -u -b -r1.11 -r1.12
--- lapack.h    10 Apr 2009 05:47:12 -0000      1.11
+++ lapack.h    14 Apr 2009 14:56:11 -0000      1.12
@@ -54,9 +54,9 @@
               double* S, double *U, int* ldu, double* VT, int* ldvt,
               double* WORK, int* lwork, int* INFO);
 
-  // Eigen decomposition of a symmetric matrix of doubles
-  void dsyev_(const char* JOBZ, const char* UPLO, int* N, double* A, int* lda, 
double* W,
-             double* WORK, int* LWORK, int* INFO);
+    // Eigen decomposition of a symmetric matrix
+    void dsyev_(const char* JOBZ, const char* UPLO, int* N, double* A, int* 
lda, double* W, double* WORK, int* LWORK, int* INFO);
+    void ssyev_(const char* JOBZ, const char* UPLO, int* N, float* A, int* 
lda, float* W, float* WORK, int* LWORK, int* INFO);
 
     // Cholesky decomposition
     void dpotrf_(const char* UPLO, const int* N, double* A, const int* LDA, 
int* INFO);
@@ -120,5 +120,12 @@
        spotri_(UPLO, N, A, LDA, INFO);
 }
 
+void syev_(const char* JOBZ, const char* UPLO, int* N, double* A, int* lda, 
double* W, double* WORK, int* LWORK, int* INFO){
+       dsyev_(JOBZ, UPLO, N, A, lda, W, WORK, LWORK, INFO);
+}
+void syev_(const char* JOBZ, const char* UPLO, int* N, float* A, int* lda, 
float* W, float* WORK, int* LWORK, int* INFO){
+       ssyev_(JOBZ, UPLO, N, A, lda, W, WORK, LWORK, INFO);
+}
+
 }
 #endif

Index: test/sym.cc
===================================================================
RCS file: test/sym.cc
diff -N test/sym.cc
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ test/sym.cc 14 Apr 2009 14:56:11 -0000      1.1
@@ -0,0 +1,47 @@
+
+#include <iostream>
+
+using namespace std;
+
+#include <TooN/SymEigen.h>
+
+int main(int, char **){
+
+    TooN::Matrix<> mt(3,3);
+    cout << mt * TooN::makeVector(0,1,2) << endl;
+
+    TooN::Matrix<2> tt;
+    tt[0] = TooN::makeVector(1,2);
+    tt[1] = TooN::makeVector(2,1);
+    
+    TooN::SymEigen<2> symt(tt);
+    cout << symt.get_evectors() << symt.get_evalues() << endl;
+
+    TooN::Matrix<3> t;
+
+    t[0] = TooN::makeVector(1,0.5, 0.5);
+    t[1] = TooN::makeVector(0.5, 2, 0.7);
+    t[2] = TooN::makeVector(0.5,0.7, 3);
+    
+    TooN::SymEigen<3> sym(t);
+    cout << sym.backsub(TooN::makeVector(0,1,2)) << endl;
+    
+    TooN::Matrix<> t2(3,3);
+
+    t2[0] = TooN::makeVector(1,0.5, 0.5);
+    t2[1] = TooN::makeVector(0.5, 2, 0.7);
+    t2[2] = TooN::makeVector(0.5,0.7, 3);
+    
+    TooN::SymEigen<> sym2(t2);
+    cout << sym2.backsub(TooN::makeVector(0,1,2)) << endl;
+
+    
+    TooN::Matrix<-1, -1, float> t2f(3,3);
+
+    t2f[0] = TooN::makeVector(1,0.5, 0.5);
+    t2f[1] = TooN::makeVector(0.5, 2, 0.7);
+    t2f[2] = TooN::makeVector(0.5,0.7, 3);
+    
+    TooN::SymEigen<TooN::Dynamic,float> symf(t2f);
+    cout << symf.backsub(TooN::makeVector(0,1,2)) << endl;
+}