[Toon-members] TooN Cholesky.h lapack.h test/chol.cc

[Gerhard Reitmayr]

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