[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN Cholesky.h helpers.h
From: |
Gerhard Reitmayr |
Subject: |
[Toon-members] TooN Cholesky.h helpers.h |
Date: |
Mon, 14 May 2007 21:04:41 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Gerhard Reitmayr <gerhard> 07/05/14 21:04:41
Modified files:
. : Cholesky.h helpers.h
Log message:
additional versions of various functions with support for DynamicMatrix
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.15&r2=1.16
http://cvs.savannah.gnu.org/viewcvs/TooN/helpers.h?cvsroot=toon&r1=1.18&r2=1.19
Patches:
Index: Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Cholesky.h,v
retrieving revision 1.15
retrieving revision 1.16
diff -u -b -r1.15 -r1.16
--- Cholesky.h 11 Feb 2007 13:20:11 -0000 1.15
+++ Cholesky.h 14 May 2007 21:04:40 -0000 1.16
@@ -48,15 +48,25 @@
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) {}
};
//
@@ -225,6 +235,8 @@
class Cholesky {
public:
+ Cholesky() {}
+
template<class Accessor>
Cholesky(const FixedMatrix<N,N,Accessor>& m){
compute(m);
@@ -329,15 +341,43 @@
}
}
}
+
+ 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 <class A1, class M2> void transform_inverse(const
DynamicMatrix<A1>& J, M2 & T) const {
+ transform_inverse<util::Assign>(J,T);
+ }
+
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;
}
+ 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;
+ }
+
template <class A1, class A2> inline
void inverse_times(const FixedVector<N,A1>& v, FixedVector<N,A2>& x)
const
{
@@ -397,6 +437,8 @@
class Cholesky<-1> {
public:
+ Cholesky(){}
+
template<class Accessor>
Cholesky(const DynamicMatrix<Accessor>& m) {
compute(m);
Index: helpers.h
===================================================================
RCS file: /cvsroot/toon/TooN/helpers.h,v
retrieving revision 1.18
retrieving revision 1.19
diff -u -b -r1.18 -r1.19
--- helpers.h 30 Jan 2007 15:57:11 -0000 1.18
+++ helpers.h 14 May 2007 21:04:40 -0000 1.19
@@ -284,7 +284,7 @@
}
}
- template <class A1, class A2, class MatM> inline void
transformCovariance(const DynamicMatrix<A1>& A, const DynamicMatrix<A2>& B,
MatM& M)
+ template <class F, class A1, class M2, class M3> inline void
transformCovarianceUpper(const DynamicMatrix<A1>& A, const M2 & B, M3 & M)
{
const int R = A.num_rows();
const int N = A.num_cols();
@@ -294,17 +294,35 @@
B.num_cols() == N);
for (int i=0; i<R; ++i) {
const Vector<> ABi = B * A[i];
- M[i][i] = ABi * A[i];
- for (int j=i+1; j<R; ++j)
- M[j][i] = M[i][j] = ABi * A[j];
+ for (int j=i; j<R; ++j)
+ F::eval(M[i][j], ABi * A[j]);
+ }
+ }
+
+ template <class F, class A1, class M2, class MatM> inline void
transformCovariance(const DynamicMatrix<A1> & A, const M2 & B, MatM& M)
+ {
+ const int R = A.num_rows();
+ const int N = A.num_cols();
+ assert(M.num_rows() == R &&
+ M.num_cols() == R &&
+ B.num_rows() == N &&
+ B.num_cols() == N);
+ for (int i=0; i<R; ++i) {
+ const Vector<> ABi = B * A[i];
+ F::eval(M[i][i], ABi * A[i]);
+ for (int j=i+1; j<R; ++j){
+ const double v = ABi * A[j];
+ F::eval(M[j][i], v);
+ F::eval(M[i][j], v);
+ }
}
}
}
- template <class A1, class A2> Matrix<> inline transformCovariance(const
DynamicMatrix<A1>& A, const DynamicMatrix<A2>& B)
+ template <class A1, class A2> Matrix<> inline transformCovariance(const
DynamicMatrix<A1> & A, const A2 & B)
{
Matrix<> M(A.num_rows(), A.num_rows());
- util::transformCovariance(A,B,M);
+ util::transformCovariance<util::Assign>(A,B,M);
return M;
}
@@ -349,6 +367,19 @@
}
}
+ template <class A, class B, class Vec> inline void add_product(const
DynamicMatrix<A>& Ma, const DynamicVector<B>& Mb, Vec & r)
+ {
+ const int M=Ma.num_rows();
+ const int N=Ma.num_cols();
+ assert(N == Mb.size());
+ for (int i=0; i<M; ++i) {
+ double sum = 0;
+ for (int k=0; k<N; ++k)
+ sum += Ma(i,k)*Mb[k];
+ r[i] += sum;
+ }
+ }
+
template <int M, int N, int C, class A1, class A2, class Mat> void
matrix_multiply(const FixedMatrix<M,N,A1>& A, const FixedMatrix<N,C,A2>& B,
Mat& R)
{
util::matrix_multiply<M,N,C>(A,B,R);
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN Cholesky.h helpers.h,
Gerhard Reitmayr <=