toon-members
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Toon-members] TooN Cholesky.h


From: Tom Drummond
Subject: [Toon-members] TooN Cholesky.h
Date: Wed, 08 Apr 2009 04:18:57 +0000

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Tom Drummond <twd20>    09/04/08 04:18:57

Modified files:
        .              : Cholesky.h 

Log message:
        added backsub for matrices, inverse and determinant

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.22&r2=1.23

Patches:
Index: Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Cholesky.h,v
retrieving revision 1.22
retrieving revision 1.23
diff -u -b -r1.22 -r1.23
--- Cholesky.h  8 Apr 2009 00:46:12 -0000       1.22
+++ Cholesky.h  8 Apr 2009 04:18:56 -0000       1.23
@@ -42,7 +42,6 @@
        template<class P2, class B2>
        Cholesky(const Matrix<Size, Size, P2, B2>& m)
                : my_cholesky(m) {
-               SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols());
                compute(m);
        }
        
@@ -56,7 +55,7 @@
                my_cholesky=m;
                int size=my_cholesky.num_rows();
                for(int col=0; col<size; col++){
-                       Precision inv_diag=1/my_cholesky(col,col);
+                       Precision inv_diag;
                        for(int row=col; row < size; row++){
                                // correct for the parts of cholesky already 
computed
                                Precision val = my_cholesky(row,col);
@@ -66,6 +65,7 @@
                                if(row==col){
                                        // this is the diagonal element so 
don't divide
                                        my_cholesky(row,col)=val;
+                                       inv_diag=1/val;
                                } else {
                                        // divide my the diagonal element for 
all others
                                        my_cholesky(row,col)=val*inv_diag;
@@ -107,492 +107,58 @@
                return result;
        }
 
-       Matrix<Size,Size,Precision> my_cholesky;
-};
-
-
-#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 {
-           template <class V1, class V2, class V3> static inline double 
eval(const V1& a, const V2& b, const V3& c) { return a[B]*b[B]*c[B] + 
Dot3<B+1,E>::eval(a,b,c); }
-       };
-       template <int B> struct Dot3<B,B> {
-           template <class V1, class V2, class V3> static inline double 
eval(const V1& a, const V2& b, const V3& c) { return a[B]*b[B]*c[B]; }
-       };
-
-       //
-       // Forward Sub
-       //
-       template <int N, int I=0> struct Forwardsub_L {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v, 
FixedVector<N,A3>& x) {
-               x[I] = v[I] - Dot<0,I-1>::eval(L[I], x);
-               Forwardsub_L<N,I+1>::eval(L, v, x);
-           }
-           template <class A1, class A2, class Vec> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const DynamicVector<A2>& v, Vec & x) {
-               x[I] = v[I] - Dot<0,I-1>::eval(L[I], x);
-               Forwardsub_L<N,I+1>::eval(L, v, x);
-           }
-       };
-       template <int N> struct Forwardsub_L<N,0> {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v, 
FixedVector<N,A3>& x) {
-               x[0] = v[0];
-               Forwardsub_L<N,1>::eval(L, v, x);
-           }
-           template <class A1, class A2, class Vec> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const DynamicVector<A2>& v, Vec & x) {
-               assert(v.size() == N && x.size() == N);
-               x[0] = v[0];
-               Forwardsub_L<N,1>::eval(L, v, x);
-           }
-       };
-       template <int N> struct Forwardsub_L<N,N> {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& /*L*/, const FixedVector<N,A2>& /*v*/, 
FixedVector<N,A3>& /*x*/) {}
-           template <class A1, class A2, class Vec> static inline void 
eval(const FixedMatrix<N,N,A1>& /*L*/, const DynamicVector<A2>& /*v*/, Vec & 
/*x*/) {}
-       };
-       
-       //
-       // Back Sub
-       //
-       template <int N, int I=N> struct Backsub_LT {
-           template <class A1, class A2, class A3, class A4> static inline 
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v, 
-                                                                               
      const FixedVector<N,A3>& invdiag, FixedVector<N,A4>& x) {
-               x[I-1] = v[I-1]*invdiag[I-1] - Dot<I,N-1>::eval(L.T()[I-1], x);
-               Backsub_LT<N,I-1>::eval(L, v, invdiag, x);
-           }
-       };
-       template <int N> struct Backsub_LT<N,N> {
-           template <class A1, class A2, class A3, class A4> static inline 
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v, 
-                                                                               
      const FixedVector<N,A3>& invdiag, FixedVector<N,A4>& x) {
-               x[N-1] = v[N-1]*invdiag[N-1];
-               Backsub_LT<N,N-1>::eval(L, v, invdiag, x);
-           }
-       };
-       template <int N> struct Backsub_LT<N,0> {
-           template <class A1, class A2, class A3, class A4> static inline 
void eval(const FixedMatrix<N,N,A1>& /*L*/, const FixedVector<N,A2>& /*v*/, 
-                                                                               
      const FixedVector<N,A3>& /*invdiag*/, FixedVector<N,A4>& /*x*/) {}
-       };
-
-       template <int N, class A1, class A2, class A3, class A4>
-       void cholesky_solve(const FixedMatrix<N,N,A1>& L, const 
FixedVector<N,A2>& invdiag, const FixedVector<N,A3>& v, FixedVector<N,A4>& x) 
-       {
-           Vector<N> t;
-           Forwardsub_L<N>::eval(L, v, t);
-           Backsub_LT<N>::eval(L, t, invdiag, x);
-       }
-
-       //
-       // Compute cholesky
-       //
-       template <int N, int I, int J=I+1> struct CholeskyInner {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L, const 
FixedVector<N,A3>& invdiag) {
-               const double a = M[I][J] - Dot<0,I-1>::eval(L[I], L.T()[J]);
-               L[I][J] = a;
-               L[J][I] = a * invdiag[I];
-               CholeskyInner<N, I, J+1>::eval(M, L, invdiag);
-           }
-       };
-
-       template <int N, int I> struct CholeskyInner<N,I,N> {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& /*M*/, FixedMatrix<N,N,A2>& /*L*/, const 
FixedVector<N,A3>& /*invdiag*/) {}
-       };
-       template <int N, int I=0> struct CholeskyOuter {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L, FixedVector<N,A3>& 
invdiag, int& rank) {
-               const double a = M[I][I] - Dot<0,I-1>::eval(L[I], L.T()[I]);
-               if (0 < a) {
-                   L[I][I] = a;
-                   invdiag[I] = 1.0/a;
-                   ++rank;
-               } else {
-                   L[I][I] = invdiag[I] = 0;
-                   return;
-               }
-               CholeskyInner<N,I>::eval(M,L,invdiag);
-               CholeskyOuter<N,I+1>::eval(M,L,invdiag,rank);
-           }
-       };
-       template <int N> struct CholeskyOuter<N,N> {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& /*M*/, FixedMatrix<N,N,A2>& /*L*/, 
FixedVector<N,A3>& /*invdiag*/, int& /*rank*/) {}
-       };
-
-
-       template <int N, class A1, class A2, class A3> int 
cholesky_compute(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L, 
FixedVector<N,A3>& invdiag)
-       {
-           int rank=0;
-           CholeskyOuter<N>::eval(M, L, invdiag, rank);
-           return rank;
-       }
-
-       //
-       // Compute inverse
-       //
-       template <int N, int Col, int I=Col+1> struct InverseForward {
-           template <class A1, class A2> static inline void eval(const 
FixedMatrix<N,N,A1>& L, FixedVector<N,A2>& t) {
-               t[I] = -(L[I][Col] + Dot<Col+1,I-1>::eval(L[I], t));
-               InverseForward<N,Col,I+1>::eval(L,t);
-           }
-       };
-       template <int N, int Col> struct InverseForward<N,Col,N> { 
-           template <class A1, class A2> static inline void eval(const 
FixedMatrix<N,N,A1>& L, FixedVector<N,A2>& t) {} 
-       };
-
-       template <int N, int Col, int I=N-1> struct InverseBack {
-           template <class A1, class A2, class A3, class A4> static inline 
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& t, 
-                                                                               
      const FixedVector<N,A3>& invdiag, FixedMatrix<N,N,A4>& Inv) {
-               Inv[I][Col] = Inv[Col][I] = invdiag[I]*t[I] - 
Dot<I+1,N-1>::eval(L.T()[I], Inv[Col]);
-               InverseBack<N,Col,I-1>::eval(L, t, invdiag, Inv);
-           }
-       };
-       template <int N, int Col> struct InverseBack<N,Col,Col> {
-           template <class A1, class A2, class A3, class A4> static inline 
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& t, 
-                                                                               
      const FixedVector<N,A3>& invdiag, FixedMatrix<N,N,A4>& Inv) {
-               Inv[Col][Col] = invdiag[Col] - Dot<Col+1,N-1>::eval(L.T()[Col], 
Inv[Col]);
-           }
-       };
-       
-       template <int N, int Col=0> struct InverseOuter {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& invdiag, 
FixedMatrix<N,N,A3>& Inv) {
-               Vector<N> t;
-               InverseForward<N,Col>::eval(L, t);
-               InverseBack<N,Col>::eval(L, t, invdiag, Inv);
-               InverseOuter<N,Col+1>::eval(L,invdiag,Inv);
-           }
-       };
-       template <int N> struct InverseOuter<N,N> {
-           template <class A1, class A2, class A3> static inline void 
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& invdiag, 
FixedMatrix<N,N,A3>& Inv) {}
-       };
-
-       template <int S, class A1, class A2, class A3> void 
cholesky_inverse(const FixedMatrix<S,S,A1>& L, const FixedVector<S,A2>& 
invdiag, FixedMatrix<S,S,A3>& I)
-       {
-           InverseOuter<S>::eval(L, invdiag, I);
-       }
-
-       //
-       // Cholesky update
-       //
-       template <int N, int I, int J=I-1> struct UpdateInner {
-           template <class LT, class PT, class FT, class ST> static inline 
void eval(LT& L, const PT& p, const FT& F, const ST sigma) {
-               const double old = L[I][J];
-               L[I][J] = old + (sigma * p[J]) * F[J];
-               UpdateInner<N,I,J-1>::eval(L, p, F, sigma + old * p[J]);
-           }
-       };
-
-       template <int N, int I> struct UpdateInner<N,I,-1> {
-           template <class LT, class PT, class FT, class ST> static inline 
void eval(LT& L, const PT& p, const FT& F, const ST sigma) {}
-       };
-           
-       template <int N, int I=1> struct UpdateOuter {
-           template <class LT, class PT, class FT> static inline void eval(LT& 
L, const PT& p, const FT& F) {
-               UpdateInner<N,I>::eval(L, p, F, p[I]);
-               UpdateOuter<N,I+1>::eval(L, p ,F);
-           }
-       };
-
-       template <int N> struct UpdateOuter<N,N> {
-           template <class LT, class PT, class FT> static inline void eval(LT& 
L, const PT& p, const FT& F) {}
-       };
-           
-       template <class P, int N, class A1, class A2> void 
cholesky_update(TooN::FixedMatrix<N,N,A1>& L, TooN::FixedVector<N,A2>& invdiag, 
const P& p)
-       {
-           double F[N-1];
-           double alpha = 1;
-           for (int i=0; i<N-1; ++i) {
-               const double p2 = p[i]*p[i];
-               const double di = L[i][i];
-               L[i][i] += alpha * p2;
-               invdiag[i] = 1.0/L[i][i];
-               const double a_over_d = alpha*invdiag[i];
-               F[i] = a_over_d;
-               alpha = di * a_over_d;
-           }
-           L[N-1][N-1] += alpha * (p[N-1]*p[N-1]);
-           invdiag[N-1] = 1.0/L[N-1][N-1];
-               
-           UpdateOuter<N>::eval(L, p, F);
-       }
-    }
-
-    template <int N=-1, typename Precision = double>
-    class Cholesky {
-    public:
-
-               Cholesky() {}
-
-               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<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; }
-               
-               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,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) {
-                               double a = A[i][i];
-                               for (int k=0; k<i; ++k)
-                                       a -= L[i][k]*L[i][k];
-                               if (0<a) {
-                                       L[i][i] = ::sqrt(a);
-                               } else {
-                                       Zero(L.slice(i,i,N-i,N-i));
-                                       return;
-                               }
-                               const double id = 1.0/L[i][i];
-                               for (int j=i+1; j<N; ++j) {
-                                       L[i][j] = 0;
-                                       double a = A[i][j];
-                                       for (int k=0; k<i; ++k)
-                                               a -= L[i][k]*L[j][k];
-                                       L[j][i] = a * id;
-                               }
-                       }
-               }
-               
-               template <class A1> static Matrix<N> sqrt(const 
FixedMatrix<N,N,A1>& A) { 
-                       Matrix<N> L;
-                       Cholesky<N>::sqrt(A,L);
-                       return L;
-               }
-
-               template <class A> void get_sqrt(FixedMatrix<N,N,A>& M) const {
-                       for (int i=0; i<N; ++i) {
-                               const double root_d = ::sqrt(L[i][i]);
-                               M[i][i] = root_d;
-                               for (int j=i+1; j<N; ++j) {
-                                       M[j][i] = L[j][i]*root_d;
-                                       M[i][j] = 0;
-                               }
-                       }
-               }
-               
-               Matrix<N> get_sqrt() const {
-                       Matrix<N> S;
-                       get_sqrt(S);
-                       return S;
-               }
-       
-               Matrix<N> get_L() const { return get_sqrt(); }
-       
-               template <class A> void get_inv_sqrt(FixedMatrix<N,N,A>& M) 
const {
-                       Vector<N> root;
-                       for (int i=0; i<N; ++i)
-                               root[i] = ::sqrt(invdiag[i]);
-                       for (int j=0; j<N; ++j) {
-                               for (int i=j+1; i<N; ++i) {
-                                       double sum = L[i][j];
-                                       for (int k=j+1; k<i; ++k)
-                                               sum += L[i][k]*M[j][k];
-                                       M[j][i] = -sum;
-                                       M[i][j] = 0;
-                               }
-                               M[j][j] = root[j];
-                               for (int i=j+1; i<N; ++i)
-                                       M[j][i] *= root[i];
-                       }
-               }
-       
-               template <class A> double mahalanobis(const FixedVector<N,A>& 
v) const {
-                       Vector<N> L_inv_v;
-                       util::Forwardsub_L<N>::eval(L, v, L_inv_v);
-                       return util::Dot3<0,N-1>::eval(L_inv_v, L_inv_v, 
invdiag);
-               }
-               
-               template <class F, int M, class A1, class A2> void 
transform_inverse(const FixedMatrix<M,N,A1>& J, FixedMatrix<M,M,A2>& T) const {
-                       Matrix<M,N> L_inv_JT;
-                       for (int i=0; i<M; ++i)
-                               util::Forwardsub_L<N>::eval(L, J[i], 
L_inv_JT[i]);
-                       for (int i=0; i<M; ++i) {               
-                               
F::eval(T[i][i],util::Dot3<0,N-1>::eval(L_inv_JT[i], L_inv_JT[i], invdiag));
-                               for (int j=i+1; j<M; ++j) {
-                                       const double x = 
util::Dot3<0,N-1>::eval(L_inv_JT[i], L_inv_JT[j], invdiag);
-                                       F::eval(T[i][j],x);
-                                       F::eval(T[j][i],x);
-                               }
-                       }
-               }
-               
-               template <class F, class A1, class M2> void 
transform_inverse(const DynamicMatrix<A1>& J, M2 & T) const {
-                       const int M = J.num_rows();
-                       assert( T.num_rows() == M && T.num_cols() == M);
-                       Matrix<> L_inv_JT(M,N);
-                       for (int i=0; i<M; ++i)
-                               util::Forwardsub_L<N>::eval(L, J[i], 
L_inv_JT[i]);
-                       for (int i=0; i<M; ++i) {
-                               
F::eval(T[i][i],util::Dot3<0,N-1>::eval(L_inv_JT[i], L_inv_JT[i], invdiag));
-                               for (int j=i+1; j<M; ++j) {
-                                       const double x = 
util::Dot3<0,N-1>::eval(L_inv_JT[i], L_inv_JT[j], invdiag);
-                                       F::eval(T[i][j],x);
-                                       F::eval(T[j][i],x);
-                               }
-                       }
-               }
                
-               template <int M, class A1, class A2> void 
transform_inverse(const FixedMatrix<M,N,A1>& J, FixedMatrix<M,M,A2>& T) const {
-                       transform_inverse<util::Assign>(J,T);
-               }
+       template<int Size2, int C2, class P2, class B2>
+       Matrix<Size, C2, Precision> backsub (const Matrix<Size2, C2, P2, B2>& 
m) {
+               int size=my_cholesky.num_rows();
+               SizeMismatch<Size,Size2>::test(size, m.num_rows());
 
-               template <class A1, class M2> void transform_inverse(const 
DynamicMatrix<A1>& J, M2 & T) const {
-                       transform_inverse<util::Assign>(J,T);
+               // first backsub through L
+               Matrix<Size, C2, Precision> y(size, m.num_cols());
+               for(int i=0; i<size; i++){
+                       Vector<C2, Precision> val = m[i];
+                       for(int j=0; j<i; j++){
+                               val -= my_cholesky(i,j)*y[j];
                }
-               
-               template <int M, class A> Matrix<M> transform_inverse(const 
FixedMatrix<M,N,A>& J) const {
-                       Matrix<M> T;
-                       transform_inverse(J,T);
-                       return T;
+                       y[i]=val;
                }
                
-               template <class A> Matrix<> transform_inverse(const 
DynamicMatrix<A>& J) const {
-                       Matrix<> T(J.num_rows(), J.num_rows());
-                       transform_inverse(J,T);
-                       return T;
+               // backsub through diagonal
+               for(int i=0; i<size; i++){
+                       y[i]*=(1/my_cholesky(i,i));
                }
 
-               template <class A1, class A2> inline 
-               void inverse_times(const FixedVector<N,A1>& v, 
FixedVector<N,A2>& x) const
-               {
-                       util::cholesky_solve(L, invdiag, v, x);
+               // backsub through L.T()
+               Matrix<Size,C2,Precision> result(size, m.num_cols());
+               for(int i=size-1; i>=0; i--){
+                       Vector<C2,Precision> val = y[i];
+                       for(int j=i+1; j<size; j++){
+                               val -= my_cholesky(j,i)*result[j];
                }
-               
-               template <class Accessor> inline 
-               Vector<N> inverse_times(const FixedVector<N,Accessor>& v) const
-               {
-                       Vector<N> x;
-                       inverse_times(v, x);
-                       return x;
+                       result[i]=val;
                }
-               
-               template <class Accessor> inline 
-               Vector<N> backsub(const FixedVector<N,Accessor>& v) const { 
return inverse_times(v); }
-               
-               template <class A, int M> inline Matrix<N,M> 
inverse_times(const FixedMatrix<N,M,A>& m)
-               {
-                       Matrix<N,M> result;
-                       for (int i=0; i<M; i++)
-                               inverse_times(m.T()[i], result.T()[i]);
                        return result;
                }
 
-               template <class A> void get_inverse(FixedMatrix<N,N,A>& M) 
const {
-                       util::cholesky_inverse(L, invdiag, M);
-               }
-               
-               Matrix<N> get_inverse() const {
-                       Matrix<N> M;
-                       get_inverse(M);
-                       return M;
-               }
-       
-               template <int M, class A>  void update(const 
FixedMatrix<N,M,A>& V) {
-                       for (int i=0; i<M; ++i) {
-                               Vector<N> p;
-                               util::Forwardsub_L<N>::eval(L, V.T()[i], p);
-                               util::cholesky_update(L, invdiag, p);
-                       }
-               }
                
-               template <class A>  void update(const FixedVector<N,A>& v) {
-                       Vector<N> p;
-                       util::Forwardsub_L<N>::eval(L, v, p);
-                       util::cholesky_update(L, invdiag, p);
+       // easy way to get inverse - could be made more efficient
+       Matrix<Size,Size,Precision> get_inverse(){
+               Matrix<Size,Size,Precision>I(Identity(my_cholesky.num_rows()));
+               return backsub(I);
                }
        
-    private:
-               Matrix<N, N, Precision> L;
-               Vector<N, Precision> invdiag;
-               int rank;
-    };
-  
-       // #endif
-
-template <int Size = -1, typename Precision = double>
-class Cholesky {
-public:
-
-       Cholesky(){}
-
-       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);
+       Precision determinant(){
+               Precision answer=my_cholesky(0,0);
+               for(int i=1; i<my_cholesky.num_rows(); i++){
+                       answer*=my_cholesky(i,i);
        }
-
-       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();
-               potrf_("L", &N, L.my_data, &N, &info);
-       }
-
-       int get_info() const { return info; }
-
-       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
-       {
-               SizeMismatch<S,Size>::test(v.size(), L.num_rows());
-               Vector<Size, Precision> x = v;
-               int N=L.num_rows();
-               int NRHS=1;
-               potrs_("L", &N, &NRHS, L.my_data, &N, x.my_data, &N, &info);    
    
-               return x;
-       }
-
-       template <int S, typename P, typename B>
-       Precision mahalanobis(const Vector<S,P,B> & v) const {
-               return v * backsub(v);
-       }
-
-       const Matrix<Size,Size,Precision>& get_L() const {
-               return L;
-       }
-
-       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;
-       }
-
-       Matrix<Size, Size, Precision> get_inverse() const {
-               Matrix<Size, Size, Precision> M = L;
-               int N = L.num_rows();
-               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;
+               return answer;
        }
 
 private:
-       Matrix<Size, Size, Precision> L;
-       mutable int info;
+       Matrix<Size,Size,Precision> my_cholesky;
 };
 
-#endif
 
 #ifndef TOON_NO_NAMESPACE
 }




reply via email to

[Prev in Thread] Current Thread [Next in Thread]