diff --git a/edo/src/utils/edoCholesky.h b/edo/src/utils/edoCholesky.h index be17f5fcd..e0437d009 100644 --- a/edo/src/utils/edoCholesky.h +++ b/edo/src/utils/edoCholesky.h @@ -48,26 +48,26 @@ public: /** Instanciate without computing anything, you are responsible of * calling the algorithm and getting the result with operator() * */ - Cholesky( size_t s1 = 1, size_t s2 = 1 ) : + CholeskyBase( size_t s1 = 1, size_t s2 = 1 ) : _L(ublas::zero_matrix(s1,s2)) {} /** Computation is made at instanciation and then cached in a member variable, * use decomposition() to get the result. */ - Cholesky(const CovarMat& V) : + CholeskyBase(const CovarMat& V) : _L(ublas::zero_matrix(V.size1(),V.size2())) { (*this)( V ); } /** Compute the factorization and cache the result */ - virtual void operator()( const CovarMat& V ) = 0; + virtual void factorize( const CovarMat& V ) = 0; /** Compute the factorization and return the result */ virtual const FactorMat& operator()( const CovarMat& V ) { - (*this)( V ); + this->factorize( V ); return decomposition(); } @@ -130,16 +130,16 @@ template< typename T > class CholeskyLLT : public CholeskyBase { public: - virtual void operator()( const CovarMat& V ) + virtual void factorize( const typename CholeskyBase::CovarMat& V ) { unsigned int N = assert_properties( V ); unsigned int i=0, j=0, k; - _L(0, 0) = sqrt( V(0, 0) ); + this->_L(0, 0) = sqrt( V(0, 0) ); // end of the column for ( j = 1; j < N; ++j ) { - _L(j, 0) = V(0, j) / _L(0, 0); + this->_L(j, 0) = V(0, j) / this->_L(0, 0); } // end of the matrix @@ -147,19 +147,19 @@ public: // diagonal double sum = 0.0; for ( k = 0; k < i; ++k) { - sum += _L(i, k) * _L(i, k); + sum += this->_L(i, k) * this->_L(i, k); } - _L(i,i) = L_i_i( V, i, sum ); + this->_L(i,i) = L_i_i( V, i, sum ); for ( j = i + 1; j < N; ++j ) { // rows // one element sum = 0.0; for ( k = 0; k < i; ++k ) { - sum += _L(j, k) * _L(i, k); + sum += this->_L(j, k) * this->_L(i, k); } - _L(j, i) = (V(j, i) - sum) / _L(i, i); + this->_L(j, i) = (V(j, i) - sum) / this->_L(i, i); } // for j in ]i,N[ } // for i in [1,N[ @@ -167,7 +167,7 @@ public: /** The step of the standard LLT algorithm where round off errors may appear */ - inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const + inline virtual T L_i_i( const typename CholeskyBase::CovarMat& V, const unsigned int& i, const double& sum ) const { // round-off errors may appear here assert( V(i,i) - sum >= 0 ); @@ -188,7 +188,7 @@ template< typename T > class CholeskyLLTabs : public CholeskyLLT { public: - inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const + inline virtual T L_i_i( const typename CholeskyBase::CovarMat& V, const unsigned int& i, const double& sum ) const { /***** ugly hack *****/ return sqrt( fabs( V(i,i) - sum) ); @@ -207,7 +207,7 @@ template< typename T > class CholeskyLLTzero : public CholeskyLLT { public: - inline virtual L_i_i( const CovarMat& V, const unsigned int& i, const double& sum ) const + inline virtual T L_i_i( const typename CholeskyBase::CovarMat& V, const unsigned int& i, const double& sum ) const { T Lii; if( V(i,i) - sum >= 0 ) { @@ -231,15 +231,15 @@ template< typename T > class CholeskyLDLT : public CholeskyBase { public: - virtual void operator()( const CovarMat& V ) + virtual void factorize( const typename CholeskyBase::CovarMat& V ) { // use "int" everywhere, because of the "j-1" operation int N = assert_properties( V ); // example of an invertible matrix whose decomposition is undefined assert( V(0,0) != 0 ); - FactorMat L = ublas::zero_matrix(N,N); - FactorMat D = ublas::zero_matrix(N,N); + typename CholeskyBase::FactorMat L = ublas::zero_matrix(N,N); + typename CholeskyBase::FactorMat D = ublas::zero_matrix(N,N); D(0,0) = V(0,0); for( int j=0; j_L = root( L, D ); } - inline FactorMat root( FactorMat& L, FactorMat& D ) + inline typename CholeskyBase::FactorMat root( typename CholeskyBase::FactorMat& L, typename CholeskyBase::FactorMat& D ) { - // now compute the final symetric matrix: _L = L D^1/2 + // now compute the final symetric matrix: this->_L = L D^1/2 // remember that V = ( L D^1/2) ( L D^1/2)^T // fortunately, the square root of a diagonal matrix is the square // root of all its elements - FactorMat sqrt_D = D; - for(int i=0; i::FactorMat sqrt_D = D; + for( int i=0; i < D.size1(); ++i) { sqrt_D(i,i) = sqrt(D(i,i)); } - // the factorization is thus _L*D^1/2 + // the factorization is thus this->_L*D^1/2 return ublas::prod( L, sqrt_D ); } };