toon-members
[Top][All Lists]
Advanced

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

[Toon-members] TooN Cholesky.h linoperators.hh maccessor.hh va...


From: Ethan Eade
Subject: [Toon-members] TooN Cholesky.h linoperators.hh maccessor.hh va...
Date: Sun, 23 Jul 2006 10:17:11 +0000

CVSROOT:        /sources/toon
Module name:    TooN
Changes by:     Ethan Eade <ethaneade>  06/07/23 10:17:11

Modified files:
        .              : Cholesky.h linoperators.hh maccessor.hh 
                         vaccessor.hh 

Log message:
        - Removed old cruft from linoperators and accessors
        - Improved performance of Cholesky on dynamic matrices

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.11&r2=1.12
http://cvs.savannah.gnu.org/viewcvs/TooN/linoperators.hh?cvsroot=toon&r1=1.10&r2=1.11
http://cvs.savannah.gnu.org/viewcvs/TooN/maccessor.hh?cvsroot=toon&r1=1.8&r2=1.9
http://cvs.savannah.gnu.org/viewcvs/TooN/vaccessor.hh?cvsroot=toon&r1=1.10&r2=1.11

Patches:
Index: Cholesky.h
===================================================================
RCS file: /sources/toon/TooN/Cholesky.h,v
retrieving revision 1.11
retrieving revision 1.12
diff -u -b -r1.11 -r1.12
--- Cholesky.h  22 Jul 2006 12:29:28 -0000      1.11
+++ Cholesky.h  23 Jul 2006 10:17:11 -0000      1.12
@@ -52,6 +52,36 @@
            }
        }
        
+       inline double dynamic_dot(const double* a, const double* b, int size)
+       {
+           double sum = 0;
+           const double* const end = a+size;
+           for (;a != end; ++a, ++b)
+               sum += *a * *b; 
+           return sum;
+       }
+
+
+       template <class A1, class A2, class V, class A4>
+       void cholesky_backsub(const DynamicMatrix<A1>& L, const 
DynamicVector<A2>& invdiag, const V& v, DynamicVector<A4>& x) 
+       {
+           const int Size = L.num_rows();
+           Vector<> t(Size);
+           // forward substitution
+           for (int i=0; i<Size; i++) {
+               const double* Li = &L(i,0);
+               t[i] = invdiag[i] * (v[i] - dynamic_dot(Li, t.get_data_ptr(), 
i));
+           }
+           // back substitution
+           for (int i=Size-1; i >=0; i--) {
+               double b = t[i];
+               const double* Lji = &L(i+1,i);
+               for (int j=i+1; j<Size;j++, Lji += Size)
+                   b -= *Lji * x[j];
+               x[i] = b*invdiag[i];
+           }
+       }
+       
        template <int Size, class A1, class A2, class A3> inline int 
cholesky_compute(const FixedMatrix<Size,Size,A1>& M, FixedMatrix<Size,Size,A2>& 
L, FixedVector<Size,A3>& invdiag) 
        {
            int rank = Size;
@@ -78,6 +108,32 @@
            return rank;
        }
 
+       template <class A1, class A2, class A3> inline int 
cholesky_compute(const DynamicMatrix<A1>& M, DynamicMatrix<A2>& L, 
DynamicVector<A3>& invdiag) 
+       {
+           const int Size = M.num_rows();
+           int rank = Size;
+           const double eps = std::numeric_limits<double>::min();
+           for(int i=0;i<Size;i++) {
+               double* const Li = &L(i,0);
+               double a = M(i,i);
+               for (int j=0; j<i; ++j) {
+                   double sum = invdiag[j] * (M(j,i) - dynamic_dot(Li, 
&L(j,0), j));
+                   Li[j] = sum;
+                   a -= sum*sum;
+               }
+               if (eps < a) {
+                   const double s = sqrt(a);
+                   Li[i]= s;
+                   invdiag[i] = 1.0/s;
+               } else {
+                   --rank;
+                   Li[i] = invdiag[i] = 0;
+               }
+           }
+           return rank;
+       }
+
+
        template <int S, class A1, class A2, class A3> inline void 
cholesky_inverse(const FixedMatrix<S,S,A1>& L, const FixedVector<S,A2>& 
invdiag, FixedMatrix<S,S,A3>& I)
        {
            for (int col = 0; col<S; col++) {
@@ -104,30 +160,29 @@
 
        template <class A1, class A2, class A3> inline void 
cholesky_inverse(const DynamicMatrix<A1>& L, const DynamicVector<A2>& invdiag, 
DynamicMatrix<A3>& I)
        {
-           assert(L.num_rows() == L.num_cols() && 
-                  L.num_rows() == invdiag.size() && 
-                  L.num_rows() == I.num_rows() && 
-                  I.num_rows() == I.num_cols());
            int S = L.num_rows();
            for (int col = 0; col<S; col++) {
                Vector<> t(S),x(S);
                t[col] = invdiag[col];
                for (int i=col+1; i<S; i++) {
+                   const double* Li = &L(i,0);
                    double psum = 0;
                    for (int j=col; j<i; j++)
-                       psum -= L[i][j]*t[j];               
+                       psum -= Li[j]*t[j];
                    t[i] = invdiag[i] * psum;
                }
                for (int i=S-1; i>col; i--) {
+                   const double* Lji = &L(i+1,i);
                    double psum = t[i];
-                   for (int j=i+1; j<S; j++)
-                       psum -= L[j][i]*x[j];
-                   I[i][col] = I[col][i] = x[i] = invdiag[i] * psum;
+                   for (int j=i+1; j<S; j++, Lji += S)
+                       psum -= *Lji * x[j];
+                   I(i,col) = I(col,i) = x[i] = invdiag[i] * psum;
                }
                double psum = t[col];
-               for (int j=col+1; j<S; j++)
-                   psum -= L[j][col]*x[j];             
-               I[col][col] = invdiag[col]*psum;
+               const double* Ljc = &L(col+1,col);
+               for (int j=col+1; j<S; j++, Ljc += S)
+                   psum -= *Ljc * x[j];                
+               I(col,col) = invdiag[col]*psum;
            }
        }
 
@@ -204,39 +259,17 @@
 
        template<class Accessor>
        void compute(const DynamicMatrix<Accessor>& m){
-           int Size = m.num_rows();
-           rank = Size;
-           for(int i=0;i<Size;i++) {
-               double a=m[i][i];
-               for(int k=0;k<i;k++) 
-                   a-=L[i][k]*L[i][k];
-               if (a < std::numeric_limits<double>::epsilon()) {
-                   --rank;
-                   L[i][i] = 0;
-               } else
-                   L[i][i]=sqrt(a);
-               invdiag[i] = 1.0/L[i][i];
-               for(int j=i+1;j<Size;j++) {
-                   a=m[i][j];
-                   for(int k=0;k<i;k++) 
-                       a-=L[i][k]*L[j][k];
-                   L[j][i]=a*invdiag[i];
-               }
-           }
+           rank = util::cholesky_compute(m,L,invdiag);
        }
        int get_rank() const { return rank; }
-       template <class Accessor> inline 
-       Vector<> backsub(const DynamicVector<Accessor>& v) const
-       {
-           assert(v.size() == L.num_rows());
-           return backsub(v.get_data_ptr());
-       }
 
-       template <int N, class Accessor> inline
-       Vector<> backsub(const FixedVector<N, Accessor>& v) const
+       template <class V>
+       Vector<> backsub(const V& v) const
        {
-           assert(N == L.num_rows());
-           return backsub(v.get_data_ptr());
+           assert(v.size() == L.num_rows());
+           Vector<> x(L.num_rows());
+           util::cholesky_backsub(L, invdiag, v, x);
+           return x;
        }
     
        const Matrix<>& get_L() const {
@@ -262,27 +295,6 @@
        }
 
     private:
-       Vector<> backsub(const double* v) const
-       {
-           int Size = L.num_rows();
-           Vector<> t(Size);
-           // forward substitution
-           for (int i=0; i<Size; i++) {
-               double b = v[i];
-               for (int j=0; j<i;j++)
-                   b -= L[i][j]*t[j];
-               t[i] = b*invdiag[i];
-           }
-           // back substitution
-           Vector<> x(Size);
-           for (int i=Size-1; i >=0; i--) {
-               double b = t[i];
-               for (int j=i+1; j<Size;j++)
-                   b -= L[j][i]*x[j];
-               x[i] = b*invdiag[i];
-           }
-           return x;
-       }
        Matrix<> L;
        Vector<> invdiag;
        int rank;

Index: linoperators.hh
===================================================================
RCS file: /sources/toon/TooN/linoperators.hh,v
retrieving revision 1.10
retrieving revision 1.11
diff -u -b -r1.10 -r1.11
--- linoperators.hh     24 May 2006 17:04:33 -0000      1.10
+++ linoperators.hh     23 Jul 2006 10:17:11 -0000      1.11
@@ -183,26 +183,6 @@
   return lhs;
 }
 
-#if 0
-template <int Size, class LH, class RHAccessor>
-inline NonConstVector<LH> operator+= (NonConstVector<LH> lhs, const 
FixedVector<Size,RHAccessor>& rhs){
-  assert(lhs.size()==Size);
-  for(int i=0; i<Size; i++){
-    lhs[i]+=rhs[i];
-  }
-  return lhs;
-}
-
-template <class LH, class RHAccessor>
-inline NonConstVector<LH> operator+= (NonConstVector<LH> lhs, const 
DynamicVector<RHAccessor>& rhs){
-  assert(lhs.size()==rhs.size());
-  for(int i=0; i<lhs.size(); i++){
-    lhs[i]+=rhs[i];
-  }
-  return lhs;
-}
-#endif
-
 
 ////////////////
 //            //
@@ -388,16 +368,6 @@
   return lhs;
 }
 
-#if 0
-template <class LH>
-inline NonConstVector<LH> operator*= (NonConstVector<LH> lhs, const double& 
rhs){
-  assert(lhs.size()==Size);
-  for(int i=0; i<Size; i++){
-    lhs[i]*=rhs;
-  }
-  return lhs;
-}
-#endif
 /////////////////////
 // operator /      //
 // Vector / double //
@@ -431,17 +401,6 @@
   return lhs*=(1/rhs);
 }
 
-#if 0
-template <class LH>
-inline NonConstVector<LH> operator/= (NonConstVector<LH> lhs, const double& 
rhs){
-  assert(lhs.size()==Size);
-  for(int i=0; i<Size; i++){
-    lhs[i]/=rhs;
-  }
-  return lhs;
-}
-#endif
-
 ///////////////////
 //  operator *   //
 // (dot product) //
@@ -668,17 +627,6 @@
   }
 }
 
-#if 0
-template <class Accessor>
-inline void operator*=(NonConstMatrix<Accessor> lhs, double rhs){
-  for(int r=0; r<lhs.num_rows(); r++){
-    for(int c=0; c<lhs.num_cols(); c++){
-      lhs[r][c]*=rhs;
-    }
-  }
-}
-#endif
-
 /////////////////////
 // operator /      //
 // Matrix / double //
@@ -705,13 +653,6 @@
   lhs*=(1.0/rhs);
 }
 
-#if 0
-template <class Accessor>
-inline void operator/=(NonConstMatrix<Accessor> lhs, double rhs){
-  lhs*=(1.0/rhs);
-}
-#endif
-
 /////////////////////
 // operator +      //
 // Matrix + Matrix //
@@ -946,33 +887,6 @@
   return lhs;
 }
 
-#if 0
-template <int Rows, int Cols, class MAccessor1, class MAccessor2>
-NonConstMatrix<MAccessor1> operator -= ( NonConstMatrix<MAccessor1> lhs,
-                                        const 
FixedMatrix<Rows,Cols,MAccessor2>& rhs){
-  assert(lhs.num_rows() == Rows && lhs.num_cols() == Cols);
-  for(int r=0; r<Rows; r++){
-    for(int c=0; c<Cols; c++){
-      lhs(r,c)-=rhs(r,c);
-    }
-  }
-  return lhs;
-}
-
-template <class MAccessor1, class MAccessor2>
-NonConstMatrix<MAccessor1> operator -= ( NonConstMatrix<MAccessor1> lhs,
-                                        const DynamicMatrix<MAccessor2>& rhs){
-  assert(lhs.num_rows() == rhs.num_cols() && lhs.num_cols() == rhs.num_cols());
-  for(int r=0; r<lhs.num_rows; r++){
-    for(int c=0; c<lhs.num_cols(); c++){
-      lhs(r,c)-=rhs(r,c);
-    }
-  }
-  return lhs;
-}
-
-#endif
-
 /////////////////////
 // operator *      //
 // Matrix * Matrix //

Index: maccessor.hh
===================================================================
RCS file: /sources/toon/TooN/maccessor.hh,v
retrieving revision 1.8
retrieving revision 1.9
diff -u -b -r1.8 -r1.9
--- maccessor.hh        3 Jul 2006 10:03:31 -0000       1.8
+++ maccessor.hh        23 Jul 2006 10:17:11 -0000      1.9
@@ -41,7 +41,7 @@
 template <>
 class RefSkipMAccessor<RowMajor> {
  public:
-  typedef NonConst<DynamicMatrix<RefSkipMAccessor<RowMajor> > > 
RefSkipMatrixRM;
+  typedef DynamicMatrix<RefSkipMAccessor<RowMajor> > RefSkipMatrixRM;
   RefSkipMAccessor(){};
   
   const RefVector operator[](int r) const TOON_THROW{
@@ -116,7 +116,7 @@
 template <>
 class RefSkipMAccessor<ColMajor> {
  public:
-  typedef NonConst<DynamicMatrix<RefSkipMAccessor<ColMajor> > > 
RefSkipMatrixCM;
+  typedef DynamicMatrix<RefSkipMAccessor<ColMajor> > RefSkipMatrixCM;
 
   RefSkipMAccessor(){};
   
@@ -328,8 +328,8 @@
   double* my_values;
 };
 
-typedef NonConst<DynamicMatrix<DynamicMAccessor<RowMajor> > > RefMatrixRM;
-typedef NonConst<DynamicMatrix<DynamicMAccessor<ColMajor> > > RefMatrixCM;
+typedef DynamicMatrix<DynamicMAccessor<RowMajor> > RefMatrixRM;
+typedef DynamicMatrix<DynamicMAccessor<ColMajor> > RefMatrixCM;
 inline RefMatrixRM makeRefMatrixRM(int nr, int nc, double* v) { RefMatrixRM 
ret; ret.set(nr,nc,v); return ret; }
 inline RefMatrixCM makeRefMatrixCM(int nr, int nc, double* v) { RefMatrixCM 
ret; ret.set(nr,nc,v); return ret; }
 

Index: vaccessor.hh
===================================================================
RCS file: /sources/toon/TooN/vaccessor.hh,v
retrieving revision 1.10
retrieving revision 1.11
diff -u -b -r1.10 -r1.11
--- vaccessor.hh        3 Jul 2006 10:03:31 -0000       1.10
+++ vaccessor.hh        23 Jul 2006 10:17:11 -0000      1.11
@@ -34,7 +34,7 @@
 class DynamicVAccessor {
   friend class VSizer;
 public:
-  typedef NonConst<DynamicVector<DynamicVAccessor> > RefVector;
+  typedef DynamicVector<DynamicVAccessor> RefVector;
     
   template<int Start, int Length>
   inline FixedVector<Length, FixedVAccessor<Length, Stack<Length> > >& slice() 
@@ -78,8 +78,8 @@
     return my_size;
   }
 
-  inline NonConst<DynamicMatrix<DynamicMAccessor<RowMajor> > > as_row(); // 
implemented in linoperators.hh
-  inline NonConst<DynamicMatrix<DynamicMAccessor<ColMajor> > > as_col(); //
+  inline DynamicMatrix<DynamicMAccessor<RowMajor> > as_row(); // implemented 
in linoperators.hh
+  inline DynamicMatrix<DynamicMAccessor<ColMajor> > as_col(); //
   inline void set(int size, double* values) { my_size = size;  my_values = 
values; }
   
  protected:  
@@ -92,7 +92,7 @@
 
 class DynamicSkipAccessor{
 public:
-  typedef NonConst<DynamicVector<DynamicSkipAccessor> > RefSkipVector;
+  typedef DynamicVector<DynamicSkipAccessor> RefSkipVector;
   
   //CHECK THIS
   template<int Start, int Length>
@@ -153,8 +153,8 @@
     return my_size;
   }
 
-  inline NonConst<DynamicMatrix<RefSkipMAccessor<ColMajor> > > as_row(); // 
implemented in linoperators.hh
-  inline NonConst<DynamicMatrix<RefSkipMAccessor<RowMajor> > > as_col(); //
+  inline DynamicMatrix<RefSkipMAccessor<ColMajor> > as_row(); // implemented 
in linoperators.hh
+  inline DynamicMatrix<RefSkipMAccessor<RowMajor> > as_col(); //
 
   inline void set(int size, int skip, double* values) { my_size = size; 
my_skip = skip; my_values = values; }
 




reply via email to

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