Compare commits

..

3 commits

3 changed files with 8 additions and 26 deletions

View file

@ -9,11 +9,8 @@ INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR})
INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS}) INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS})
LINK_DIRECTORIES(${Boost_LIBRARY_DIRS}) LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
INCLUDE_DIRECTORIES("/usr/include/")
LINK_DIRECTORIES("/usr/lib/")
ADD_DEFINITIONS( -Wall -W -Wextra ) ADD_DEFINITIONS( -Wall -W -Wextra )
ADD_EXECUTABLE(test test.cpp) ADD_EXECUTABLE(test test.cpp)
TARGET_LINK_LIBRARIES(test ${Boost_LIBRARIES} gmpxx) TARGET_LINK_LIBRARIES(test ${Boost_LIBRARIES})

View file

@ -55,9 +55,6 @@ public:
// FIXME check if triangular types behaviour is like having 0 // FIXME check if triangular types behaviour is like having 0
typedef ublas::matrix< T > FactorMat; typedef ublas::matrix< T > FactorMat;
//! The type of the scalars
typedef T AtomType;
/** Instanciate without computing anything, you are responsible of /** Instanciate without computing anything, you are responsible of
* calling the algorithm and getting the result with operator() * calling the algorithm and getting the result with operator()
* */ * */
@ -208,8 +205,7 @@ public:
inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const T& sum ) const inline virtual T L_i_i( const typename Cholesky<T>::CovarMat& V, const unsigned int& i, const T& sum ) const
{ {
/***** ugly hack *****/ /***** ugly hack *****/
//return sqrt( fabs( V(i,i) - sum) ); return sqrt( fabs( V(i,i) - sum) );
return sqrt( abs( V(i,i) - sum) );
} }
}; };
@ -326,8 +322,8 @@ public:
T sum = 0.0; T sum = 0.0;
// for( unsigned int k=1; k<i-1; ++k) { // for( unsigned int k=1; k<i-1; ++k) {
// for( unsigned int k=1; k<i; ++k) { // for( unsigned int k=1; k<i; ++k) {
// for( unsigned int k=0; k<i; ++k) { for( unsigned int k=0; k<i; ++k) {
for( unsigned int k=0; k<i-1; ++k) { // for( unsigned int k=0; k<i-1; ++k) {
sum += this->_L(i, k) * this->_L(i, k); sum += this->_L(i, k) * this->_L(i, k);
} }
this->_L(i,i) = this->_L(i,i) - sum; this->_L(i,i) = this->_L(i,i) - sum;
@ -342,8 +338,8 @@ public:
sum = 0.0; sum = 0.0;
// for( unsigned int k = 1; k < i-1; ++k ) { // for( unsigned int k = 1; k < i-1; ++k ) {
// for( unsigned int k = 1; k < i; ++k ) { // for( unsigned int k = 1; k < i; ++k ) {
// for( unsigned int k = 0; k < i; ++k ) { for( unsigned int k = 0; k < i; ++k ) {
for( unsigned int k = 0; k < i-1; ++k ) { // for( unsigned int k = 0; k < i-1; ++k ) {
sum += this->_L(j, k) * this->_L(i, k); sum += this->_L(j, k) * this->_L(i, k);
} }
this->_L(j,i) = this->_L(j,i) - sum; this->_L(j,i) = this->_L(j,i) - sum;

View file

@ -28,8 +28,6 @@ Authors:
#include <map> #include <map>
#include <vector> #include <vector>
#include <gmpxx.h>
#include "cholesky.h" #include "cholesky.h"
@ -114,8 +112,7 @@ T trigsum( const MT& M )
T sum = 0; T sum = 0;
for( unsigned int i=0; i<M.size1(); ++i ) { for( unsigned int i=0; i<M.size1(); ++i ) {
for( unsigned int j=i; j<M.size2(); ++j ) { // triangular browsing for( unsigned int j=i; j<M.size2(); ++j ) { // triangular browsing
//sum += fabs( M(i,j) ); // absolute deviation sum += fabs( M(i,j) ); // absolute deviation
sum += abs( M(i,j) ); // absolute deviation
} }
} }
return sum; return sum;
@ -168,9 +165,7 @@ void test( unsigned int M, unsigned int N, unsigned int F, unsigned int R, unsig
// a variance-covariance matrix of size N*N // a variance-covariance matrix of size N*N
// Note: a covariance matrix is necessarily semi-positive definite // Note: a covariance matrix is necessarily semi-positive definite
// thus, any failure in the Cholesky factorization is due to round-off errors // thus, any failure in the Cholesky factorization is due to round-off errors
//CovarMat V = ublas::prod( ublas::trans(S), S ); CovarMat V = ublas::prod( ublas::trans(S), S );
CovarMat ST = ublas::trans(S);
CovarMat V = ublas::prod( ST, S );
assert( V.size1() == N && V.size2() == N ); assert( V.size1() == N && V.size2() == N );
#ifndef NDEBUG #ifndef NDEBUG
if( R == 1 ) { if( R == 1 ) {
@ -243,7 +238,6 @@ int main(int argc, char** argv)
std::clog << "Legend:" << std::endl; std::clog << "Legend:" << std::endl;
std::clog << "\tAlgo: (failures/runs)\tAverage error" << std::endl; std::clog << "\tAlgo: (failures/runs)\tAverage error" << std::endl;
/*
std::cout << std::endl << "FLOAT" << std::endl; std::cout << std::endl << "FLOAT" << std::endl;
test<float>(M,N,F,R,seed); test<float>(M,N,F,R,seed);
@ -252,10 +246,5 @@ int main(int argc, char** argv)
std::cout << std::endl << "LONG DOUBLE" << std::endl; std::cout << std::endl << "LONG DOUBLE" << std::endl;
test<long double>(M,N,F,R,seed); test<long double>(M,N,F,R,seed);
*/
std::cout << std::endl << "MPF 128" << std::endl;
mpf_set_default_prec(128);
test<mpf_class>(M,N,F,R,seed);
} }