cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

boostmultiprec.cpp (5731B)


      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #include <sstream>
     11 
     12 #ifdef EIGEN_TEST_MAX_SIZE
     13 #undef EIGEN_TEST_MAX_SIZE
     14 #endif
     15 
     16 #define EIGEN_TEST_MAX_SIZE 50
     17 
     18 #ifdef EIGEN_TEST_PART_1
     19 #include "cholesky.cpp"
     20 #endif
     21 
     22 #ifdef EIGEN_TEST_PART_2
     23 #include "lu.cpp"
     24 #endif
     25 
     26 #ifdef EIGEN_TEST_PART_3
     27 #include "qr.cpp"
     28 #endif
     29 
     30 #ifdef EIGEN_TEST_PART_4
     31 #include "qr_colpivoting.cpp"
     32 #endif
     33 
     34 #ifdef EIGEN_TEST_PART_5
     35 #include "qr_fullpivoting.cpp"
     36 #endif
     37 
     38 #ifdef EIGEN_TEST_PART_6
     39 #include "eigensolver_selfadjoint.cpp"
     40 #endif
     41 
     42 #ifdef EIGEN_TEST_PART_7
     43 #include "eigensolver_generic.cpp"
     44 #endif
     45 
     46 #ifdef EIGEN_TEST_PART_8
     47 #include "eigensolver_generalized_real.cpp"
     48 #endif
     49 
     50 #ifdef EIGEN_TEST_PART_9
     51 #include "jacobisvd.cpp"
     52 #endif
     53 
     54 #ifdef EIGEN_TEST_PART_10
     55 #include "bdcsvd.cpp"
     56 #endif
     57 
     58 #ifdef EIGEN_TEST_PART_11
     59 #include "simplicial_cholesky.cpp"
     60 #endif
     61 
     62 #include <Eigen/Dense>
     63 
     64 #undef min
     65 #undef max
     66 #undef isnan
     67 #undef isinf
     68 #undef isfinite
     69 #undef I
     70 
     71 #include <boost/serialization/nvp.hpp>
     72 #include <boost/multiprecision/cpp_dec_float.hpp>
     73 #include <boost/multiprecision/number.hpp>
     74 #include <boost/math/special_functions.hpp>
     75 #include <boost/math/complex.hpp>
     76 
     77 namespace mp = boost::multiprecision;
     78 typedef mp::number<mp::cpp_dec_float<100>, mp::et_on> Real;
     79 
     80 namespace Eigen {
     81   template<> struct NumTraits<Real> : GenericNumTraits<Real> {
     82     static inline Real dummy_precision() { return 1e-50; }
     83   };
     84 
     85   template<typename T1,typename T2,typename T3,typename T4,typename T5>
     86   struct NumTraits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> > : NumTraits<Real> {};
     87 
     88   template<>
     89   Real test_precision<Real>() { return 1e-50; }
     90 
     91   // needed in C++93 mode where number does not support explicit cast.
     92   namespace internal {
     93     template<typename NewType>
     94     struct cast_impl<Real,NewType> {
     95       static inline NewType run(const Real& x) {
     96         return x.template convert_to<NewType>();
     97       }
     98     };
     99 
    100     template<>
    101     struct cast_impl<Real,std::complex<Real> > {
    102       static inline std::complex<Real>  run(const Real& x) {
    103         return std::complex<Real>(x);
    104       }
    105     };
    106   }
    107 }
    108 
    109 namespace boost {
    110 namespace multiprecision {
    111   // to make ADL works as expected:
    112   using boost::math::isfinite;
    113   using boost::math::isnan;
    114   using boost::math::isinf;
    115   using boost::math::copysign;
    116   using boost::math::hypot;
    117 
    118   // The following is needed for std::complex<Real>:
    119   Real fabs(const Real& a) { return abs EIGEN_NOT_A_MACRO (a); }
    120   Real fmax(const Real& a, const Real& b) { using std::max; return max(a,b); }
    121 
    122   // some specialization for the unit tests:
    123   inline bool test_isMuchSmallerThan(const Real& a, const Real& b) {
    124     return internal::isMuchSmallerThan(a, b, test_precision<Real>());
    125   }
    126 
    127   inline bool test_isApprox(const Real& a, const Real& b) {
    128     return internal::isApprox(a, b, test_precision<Real>());
    129   }
    130 
    131   inline bool test_isApproxOrLessThan(const Real& a, const Real& b) {
    132     return internal::isApproxOrLessThan(a, b, test_precision<Real>());
    133   }
    134 
    135   Real get_test_precision(const Real&) {
    136     return test_precision<Real>();
    137   }
    138 
    139   Real test_relative_error(const Real &a, const Real &b) {
    140     using Eigen::numext::abs2;
    141     return sqrt(abs2<Real>(a-b)/Eigen::numext::mini<Real>(abs2(a),abs2(b)));
    142   }
    143 }
    144 }
    145 
    146 namespace Eigen {
    147 
    148 }
    149 
    150 EIGEN_DECLARE_TEST(boostmultiprec)
    151 {
    152   typedef Matrix<Real,Dynamic,Dynamic> Mat;
    153   typedef Matrix<std::complex<Real>,Dynamic,Dynamic> MatC;
    154 
    155   std::cout << "NumTraits<Real>::epsilon()         = " << NumTraits<Real>::epsilon() << std::endl;
    156   std::cout << "NumTraits<Real>::dummy_precision() = " << NumTraits<Real>::dummy_precision() << std::endl;
    157   std::cout << "NumTraits<Real>::lowest()          = " << NumTraits<Real>::lowest() << std::endl;
    158   std::cout << "NumTraits<Real>::highest()         = " << NumTraits<Real>::highest() << std::endl;
    159   std::cout << "NumTraits<Real>::digits10()        = " << NumTraits<Real>::digits10() << std::endl;
    160 
    161   // check stream output
    162   {
    163     Mat A(10,10);
    164     A.setRandom();
    165     std::stringstream ss;
    166     ss << A;
    167   }
    168   {
    169     MatC A(10,10);
    170     A.setRandom();
    171     std::stringstream ss;
    172     ss << A;
    173   }
    174 
    175   for(int i = 0; i < g_repeat; i++) {
    176     int s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
    177 
    178     CALL_SUBTEST_1( cholesky(Mat(s,s)) );
    179 
    180     CALL_SUBTEST_2( lu_non_invertible<Mat>() );
    181     CALL_SUBTEST_2( lu_invertible<Mat>() );
    182     CALL_SUBTEST_2( lu_non_invertible<MatC>() );
    183     CALL_SUBTEST_2( lu_invertible<MatC>() );
    184 
    185     CALL_SUBTEST_3( qr(Mat(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
    186     CALL_SUBTEST_3( qr_invertible<Mat>() );
    187 
    188     CALL_SUBTEST_4( qr<Mat>() );
    189     CALL_SUBTEST_4( cod<Mat>() );
    190     CALL_SUBTEST_4( qr_invertible<Mat>() );
    191 
    192     CALL_SUBTEST_5( qr<Mat>() );
    193     CALL_SUBTEST_5( qr_invertible<Mat>() );
    194 
    195     CALL_SUBTEST_6( selfadjointeigensolver(Mat(s,s)) );
    196 
    197     CALL_SUBTEST_7( eigensolver(Mat(s,s)) );
    198 
    199     CALL_SUBTEST_8( generalized_eigensolver_real(Mat(s,s)) );
    200 
    201     TEST_SET_BUT_UNUSED_VARIABLE(s)
    202   }
    203 
    204   CALL_SUBTEST_9(( jacobisvd(Mat(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
    205   CALL_SUBTEST_10(( bdcsvd(Mat(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
    206 
    207   CALL_SUBTEST_11(( test_simplicial_cholesky_T<Real,int,ColMajor>() ));
    208 }