use gmpxx instead of mpfr++

CMakeLists.txt

@@ -0,0 +1,19 @@
+
+PROJECT(Cholesky)
+ 
+cmake_minimum_required(VERSION 2.8)
+FIND_PACKAGE(Boost 1.33.0)
+
+INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR})
+
+INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS})
+LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
+
+INCLUDE_DIRECTORIES("/usr/include/")
+LINK_DIRECTORIES("/usr/lib/")
+
+ADD_DEFINITIONS( -Wall -W -Wextra )
+
+ADD_EXECUTABLE(test test.cpp)
+TARGET_LINK_LIBRARIES(test ${Boost_LIBRARIES} gmpxx)
+

cholesky.h

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

test.cpp

@@ -28,24 +28,27 @@ Authors:
 #include <map>
 #include <vector>
 
+#include <gmpxx.h>
+
 #include "cholesky.h"
 
 
+template<typename T>
 void setformat( std::ostream& out )
 {
     out << std::right;
     out << std::setfill(' ');
-    out << std::setw( 5 + std::numeric_limits<double>::digits10);
-    out << std::setprecision(std::numeric_limits<double>::digits10);
+    out << std::setw( 5 + std::numeric_limits<T>::digits10);
+    out << std::setprecision(std::numeric_limits<T>::digits10);
     out << std::setiosflags(std::ios_base::showpoint);
 }
 
 
-template<typename MT>
+template<typename T, typename MT>
 std::string format(const MT& mat )
 {
     std::ostringstream out;
-    setformat(out);
+    setformat<T>(out);
 
     for( unsigned int i=0; i<mat.size1(); ++i) {
         for( unsigned int j=0; j<mat.size2(); ++j) {
@@ -67,8 +70,8 @@ T round( T val, T prec = 1.0 )
 }
 
 
-template<typename MT>
-bool equal( const MT& M1, const MT& M2, double prec /* = 1/std::numeric_limits<double>::digits10 ???*/ )
+template<typename T, typename MT>
+bool equal( const MT& M1, const MT& M2, T prec )
 {
     if( M1.size1() != M2.size1() || M1.size2() != M2.size2() ) {
         return false;
@@ -88,12 +91,12 @@ bool equal( const MT& M1, const MT& M2, double prec /* = 1/std::numeric_limits<d
 }
 
 
-template<typename MT >
+template<typename T, typename MT>
 MT error( const MT& M1, const MT& M2 )
 {
     assert( M1.size1() == M2.size1() && M1.size1() == M2.size2() );
 
-    MT Err = ublas::zero_matrix<double>(M1.size1(),M1.size2());
+    MT Err = ublas::zero_matrix<T>(M1.size1(),M1.size2());
 
     for( unsigned int i=0; i<M1.size1(); ++i ) {
         for( unsigned int j=0; j<M1.size2(); ++j ) {
@@ -105,23 +108,24 @@ MT error( const MT& M1, const MT& M2 )
 }
 
 
-template<typename MT >
-double trigsum( const MT& M )
+template<typename T, typename MT>
+T trigsum( const MT& M )
 {
-    double sum = 0;
+    T sum = 0;
     for( unsigned int i=0; i<M.size1(); ++i ) {
         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;
 }
 
 
-template<typename T>
-double sum( const T& c )
+template<typename T,typename V>
+T sum( const V& c )
 {
-     return std::accumulate(c.begin(), c.end(), 0);
+     return std::accumulate(c.begin(), c.end(), static_cast<T>(0) );
 }
 
 
@@ -145,10 +149,10 @@ void test( unsigned int M, unsigned int N, unsigned int F, unsigned int R, unsig
     // init data structures on the same keys than given algorithms
     std::map<std::string,unsigned int> fails;
     // triangular errors sum
-    std::map<std::string, std::vector<double> > errors;
+    std::map<std::string, std::vector<T> > errors;
     for( typename AlgoMap::iterator ialgo = algos.begin(); ialgo != algos.end(); ++ialgo ) {
          fails[ialgo->first] = 0;
-        errors[ialgo->first] = std::vector<double>();
+        errors[ialgo->first] = std::vector<T>();
     }
 
     for( unsigned int n=0; n<R; ++n ) {
@@ -164,12 +168,14 @@ void test( unsigned int M, unsigned int N, unsigned int F, unsigned int R, unsig
         // a variance-covariance matrix of size N*N
         // Note: a covariance matrix is necessarily semi-positive definite
         //       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 );
 #ifndef NDEBUG
         if( R == 1 ) {
             std::cout << std::endl << "Covariance matrix:" << std::endl;
-            std::cout << format(V) << std::endl;
+            std::cout << format<T>(V) << std::endl;
         }
 #endif
         for( typename AlgoMap::iterator ialgo = algos.begin(); ialgo != algos.end(); ++ialgo ) {
@@ -179,7 +185,7 @@ void test( unsigned int M, unsigned int N, unsigned int F, unsigned int R, unsig
             try {
                 L = (*ialgo->second)(V);
                 CovarMat Vn = ublas::prod( L, ublas::trans(L) );
-                errors[ialgo->first].push_back( trigsum(error(V,Vn)) );
+                errors[ialgo->first].push_back( trigsum<T>(error<T>(V,Vn)) );
 
             } catch( cholesky::NotDefinitePositive & error ) {
                 fails[ialgo->first]++;
@@ -201,7 +207,7 @@ void test( unsigned int M, unsigned int N, unsigned int F, unsigned int R, unsig
             std::cout << "NAN";
         } else {
             assert( errors[a].size() == R - fails[a] );
-            std::cout << sum(errors[a])/R;
+            std::cout << sum<T>(errors[a])/R;
         }
         std::cout << std::endl;
     } // for a in algos
@@ -237,6 +243,7 @@ int main(int argc, char** argv)
     std::clog << "Legend:" << std::endl;
     std::clog << "\tAlgo: (failures/runs)\tAverage error" << std::endl;
 
+    /*
     std::cout << std::endl << "FLOAT" << std::endl;
     test<float>(M,N,F,R,seed);
 
@@ -245,4 +252,10 @@ int main(int argc, char** argv)
 
     std::cout << std::endl << "LONG DOUBLE" << std::endl;
     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);
 }
+