[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN SymEigen.h lapack.h test/sym.cc
From: |
Gerhard Reitmayr |
Subject: |
[Toon-members] TooN SymEigen.h lapack.h test/sym.cc |
Date: |
Tue, 14 Apr 2009 14:56:12 +0000 |
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;
+}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN SymEigen.h lapack.h test/sym.cc,
Gerhard Reitmayr <=