[Top][All Lists]
[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; }
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN Cholesky.h linoperators.hh maccessor.hh va...,
Ethan Eade <=