[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN Cholesky.h lapack.h test/chol.cc
From: |
Gerhard Reitmayr |
Subject: |
[Toon-members] TooN Cholesky.h lapack.h test/chol.cc |
Date: |
Fri, 03 Apr 2009 18:28:46 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Gerhard Reitmayr <gerhard> 09/04/03 18:28:45
Modified files:
. : Cholesky.h lapack.h
Added files:
test : chol.cc
Log message:
first stab at porting Cholesky to TooN2, general LAPACK-based
implementation for static and dynamic matrices works
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.19&r2=1.20
http://cvs.savannah.gnu.org/viewcvs/TooN/lapack.h?cvsroot=toon&r1=1.9&r2=1.10
http://cvs.savannah.gnu.org/viewcvs/TooN/test/chol.cc?cvsroot=toon&rev=1.1
Patches:
Index: Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Cholesky.h,v
retrieving revision 1.19
retrieving revision 1.20
diff -u -b -r1.19 -r1.20
--- Cholesky.h 10 Mar 2009 16:05:26 -0000 1.19
+++ Cholesky.h 3 Apr 2009 18:28:45 -0000 1.20
@@ -20,8 +20,6 @@
#ifndef CHOLESKY_H
#define CHOLESKY_H
-#include <iostream>
-
#include <TooN/lapack.h>
#include <TooN/TooN.h>
@@ -31,6 +29,8 @@
#ifndef TOON_NO_NAMESPACE
namespace TooN {
#endif
+
+#if 0 // should be removed and possibly replaced with wls_cholesky
implementation for small fixed sizes
namespace util {
template <int B, int E> struct Dot3 {
@@ -231,31 +231,40 @@
}
}
- template <int N=-1>
+ template <int N=-1, typename Precision = double>
class Cholesky {
public:
Cholesky() {}
- template<class Accessor>
- Cholesky(const FixedMatrix<N,N,Accessor>& m){
+ template<int R, int C, typename P, typename A>
+ Cholesky(const Matrix<R,C,P,A>& m) : L(m.num_rows(),
m.num_cols()), invdiag(m.num_cols()) {
compute(m);
}
- template<class Accessor>
- void compute(const FixedMatrix<N,N,Accessor>& m){
- rank = util::cholesky_compute(m,L,invdiag);
+ template<int R, int C, typename P, typename A>
+ void compute(const Matrix<R,C,P,A>& m){
+ SizeMismatch<R,C>::test(m.num_rows(), m.num_cols());
+ SizeMismatch<R,N>::test(m.num_rows(), L.num_rows());
+ // rank = util::cholesky_compute(m,L,invdiag);
+ L = m;
+ int Size = L.num_rows();
+ int info;
+ dpotrf_("L", &Size, L.get_data_ptr(), &Size, &info);
+ assert(info >= 0);
+ if (info > 0)
+ rank = info-1;
}
int get_rank() const { return rank; }
- double get_determinant() const {
- double det = L[0][0];
- for (int i=1; i<N; i++)
- det *= L[i][i];
+ Precision get_determinant() const {
+ Precision det = L(0,0);
+ for (int i=1; i<L.num_rows(); ++i)
+ det *= L(i,i);
return det;
}
- const Matrix<N>& get_L_D() const { return L; }
+ const Matrix<N,N,Precision> & get_L_D() const { return L; }
template <class A1, class A2> static void sqrt(const
FixedMatrix<N,N,A1>& A, FixedMatrix<N,N,A2>& L) {
for (int i=0; i<N; ++i) {
@@ -428,94 +437,83 @@
}
private:
- Matrix<N> L;
- Vector<N> invdiag;
+ Matrix<N, N, Precision> L;
+ Vector<N, Precision> invdiag;
int rank;
};
+#endif
- template <>
- class Cholesky<-1> {
- public:
+template <int Size = -1, typename Precision = double>
+class Cholesky {
+public:
Cholesky(){}
- template<class Accessor>
- Cholesky(const DynamicMatrix<Accessor>& m) {
+ template<int R, int C, typename P, typename B>
+ Cholesky(const Matrix<R,C,P,B> & m) : L(m.num_rows(), m.num_cols()) {
compute(m);
}
- template<class Accessor>
- void compute(const DynamicMatrix<Accessor>& m){
- assert(m.num_rows() == m.num_cols());
- L.assign(m);
+ template<int R, int C, typename P, typename B>
+ void compute(const Matrix<R,C,P,B> & m){
+ SizeMismatch<R,C>::test(m.num_rows(), m.num_cols());
+ SizeMismatch<R,Size>::test(m.num_rows(), L.num_rows());
+ L = m;
int N = L.num_rows();
- int info;
- dpotrf_("L", &N, L.get_data_ptr(), &N, &info);
- assert(info >= 0);
- if (info > 0)
- rank = info-1;
+ potrf_("L", &N, L.my_data, &N, &info);
}
- int get_rank() const { return rank; }
- template <class V> inline
- Vector<> inverse_times(const V& v) const { return backsub(v); }
+ int get_info() const { return info; }
- template <class V> inline
- Vector<> backsub(const V& v) const
+ template <int S, typename P, typename B>
+ Vector<Size, Precision> inverse_times(const Vector<S,P,B> & v) const {
return backsub(v); }
+
+ template <int S, typename P, typename B>
+ Vector<Size, Precision> backsub(const Vector<S,P,B> & v) const
{
- assert(v.size() == L.num_rows());
- Vector<> x = v;
+ SizeMismatch<S,Size>::test(v.size(), L.num_rows());
+ Vector<Size, Precision> x = v;
int N=L.num_rows();
int NRHS=1;
- int info;
- dpotrs_("L", &N, &NRHS, L.get_data_ptr(), &N,
x.get_data_ptr(), &N, &info);
- assert(info==0);
+ potrs_("L", &N, &NRHS, L.my_data, &N, x.my_data, &N, &info);
return x;
}
- template <class V> double mahalanobis(const V& v) const {
+ template <int S, typename P, typename B>
+ Precision mahalanobis(const Vector<S,P,B> & v) const {
return v * backsub(v);
}
- const Matrix<>& get_L() const {
+ const Matrix<Size,Size,Precision>& get_L() const {
return L;
}
- double get_determinant() const {
- double det = L[0][0];
- for (int i=1; i<L.num_rows(); i++)
- det *= L[i][i];
+ Precision get_determinant() const {
+ Precision det = L(0,0);
+ for (int i=1; i<L.num_rows(); ++i)
+ det *= L(i,i);
return det*det;
}
- template <class Mat> void get_inverse(Mat& M) const {
- M = L;
+ Matrix<Size, Size, Precision> get_inverse() const {
+ Matrix<Size, Size, Precision> M = L;
int N = L.num_rows();
- int info;
- dpotri_("L", &N, M.get_data_ptr(), &N, &info);
- assert(info == 0);
- TooN::Symmetrize(M);
- }
-
- Matrix<> get_inverse() const {
- Matrix<> M(L.num_rows(), L.num_rows());
- get_inverse(M);
+ potri_("L", &N, M.my_data, &N, &info);
+ for(int i = 1; i < M.num_rows(); ++i)
+ for(int j = 0; j < i; ++j)
+ M(i,j) = M(j,i);
return M;
}
- private:
- Matrix<> L;
- int rank;
- };
+private:
+ Matrix<Size, Size, Precision> L;
+ mutable int info;
+};
#ifndef TOON_NO_NAMESPACE
}
#endif
-
-
-
-
#endif
Index: lapack.h
===================================================================
RCS file: /cvsroot/toon/TooN/lapack.h,v
retrieving revision 1.9
retrieving revision 1.10
diff -u -b -r1.9 -r1.10
--- lapack.h 20 Mar 2009 16:49:25 -0000 1.9
+++ lapack.h 3 Apr 2009 18:28:45 -0000 1.10
@@ -50,12 +50,15 @@
// Cholesky decomposition
void dpotrf_(const char* UPLO, const int* N, double* A, const int* LDA,
int* INFO);
+ void spotrf_(const char* UPLO, const int* N, float* A, const int* LDA,
int* INFO);
// Cholesky solve AX=B given decomposition
void dpotrs_(const char* UPLO, const int* N, const int* NRHS, const
double* A, const int* LDA, double* B, const int* LDB, int* INFO);
+ void spotrs_(const char* UPLO, const int* N, const int* NRHS, const float*
A, const int* LDA, float* B, const int* LDB, int* INFO);
// Cholesky inverse given decomposition
void dpotri_(const char* UPLO, const int* N, double* A, const int* LDA,
int* INFO);
+ void spotri_(const char* UPLO, const int* N, float* A, const int* LDA,
int* INFO);
}
void getrf_(int* M, int *N, float* A, int* lda, int* IPIV, int* INFO){
@@ -82,5 +85,30 @@
sgetri_(N, A, lda, IPIV, WORK, lwork, INFO);
}
+void potrf_(const char * UPLO, const int* N, double* A, const int* LDA, int*
INFO){
+ dpotrf_(UPLO, N, A, LDA, INFO);
+}
+
+void potrf_(const char * UPLO, const int* N, float* A, const int* LDA, int*
INFO){
+ spotrf_(UPLO, N, A, LDA, INFO);
+}
+// Cholesky solve AX=B given decomposition
+void potrs_(const char* UPLO, const int* N, const int* NRHS, const double* A,
const int* LDA, double* B, const int* LDB, int* INFO){
+ dpotrs_(UPLO, N, NRHS, A, LDA, B, LDB, INFO);
+}
+
+void potrs_(const char* UPLO, const int* N, const int* NRHS, const float* A,
const int* LDA, float* B, const int* LDB, int* INFO){
+ spotrs_(UPLO, N, NRHS, A, LDA, B, LDB, INFO);
+}
+
+// Cholesky inverse given decomposition
+void potri_(const char* UPLO, const int* N, double* A, const int* LDA, int*
INFO){
+ dpotri_(UPLO, N, A, LDA, INFO);
+}
+
+void potri_(const char* UPLO, const int* N, float* A, const int* LDA, int*
INFO){
+ spotri_(UPLO, N, A, LDA, INFO);
+}
+
}
#endif
Index: test/chol.cc
===================================================================
RCS file: test/chol.cc
diff -N test/chol.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ test/chol.cc 3 Apr 2009 18:28:45 -0000 1.1
@@ -0,0 +1,37 @@
+#include <TooN/Cholesky.h>
+
+#include <iostream>
+
+using namespace std;
+
+int main(int, char ** ){
+
+ 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::Cholesky<3> chol(t);
+ cout << chol.get_determinant() << endl;
+
+ cout << t << "\n" << chol.get_inverse() << "\n" << t * chol.get_inverse()
<< endl;
+ cout << chol.get_L() << 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::Cholesky<-1,float> chol2(t2);
+
+ cout << t2 << "\n" << chol2.get_inverse() << "\n" << t2 *
chol2.get_inverse() << endl;
+ cout << chol2.get_L() << endl;
+
+ TooN::Vector<3> bla = TooN::makeVector(1,2,3);
+
+ cout << chol.backsub(bla) << endl;
+ cout << chol.get_inverse() * bla << endl;
+ cout << chol.mahalanobis(bla) - bla * chol.backsub(bla) << endl;
+
+ return 0;
+}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN Cholesky.h lapack.h test/chol.cc,
Gerhard Reitmayr <=