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#include <Eigen/Dense>
59
60#undef min
61#undef max
62#undef isnan
63#undef isinf
64#undef isfinite
65
66#include <boost/multiprecision/cpp_dec_float.hpp>
67#include <boost/multiprecision/number.hpp>
68#include <boost/math/special_functions.hpp>
69#include <boost/math/complex.hpp>
70
71namespace mp = boost::multiprecision;
72typedef mp::number<mp::cpp_dec_float<100>, mp::et_on> Real;
73
74namespace Eigen {
75  template<> struct NumTraits<Real> : GenericNumTraits<Real> {
76    static inline Real dummy_precision() { return 1e-50; }
77  };
78
79  template<typename T1,typename T2,typename T3,typename T4,typename T5>
80  struct NumTraits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> > : NumTraits<Real> {};
81
82  template<>
83  Real test_precision<Real>() { return 1e-50; }
84
85  // needed in C++93 mode where number does not support explicit cast.
86  namespace internal {
87    template<typename NewType>
88    struct cast_impl<Real,NewType> {
89      static inline NewType run(const Real& x) {
90        return x.template convert_to<NewType>();
91      }
92    };
93
94    template<>
95    struct cast_impl<Real,std::complex<Real> > {
96      static inline std::complex<Real>  run(const Real& x) {
97        return std::complex<Real>(x);
98      }
99    };
100  }
101}
102
103namespace boost {
104namespace multiprecision {
105  // to make ADL works as expected:
106  using boost::math::isfinite;
107  using boost::math::isnan;
108  using boost::math::isinf;
109  using boost::math::copysign;
110  using boost::math::hypot;
111
112  // The following is needed for std::complex<Real>:
113  Real fabs(const Real& a) { return abs EIGEN_NOT_A_MACRO (a); }
114  Real fmax(const Real& a, const Real& b) { using std::max; return max(a,b); }
115
116  // some specialization for the unit tests:
117  inline bool test_isMuchSmallerThan(const Real& a, const Real& b) {
118    return internal::isMuchSmallerThan(a, b, test_precision<Real>());
119  }
120
121  inline bool test_isApprox(const Real& a, const Real& b) {
122    return internal::isApprox(a, b, test_precision<Real>());
123  }
124
125  inline bool test_isApproxOrLessThan(const Real& a, const Real& b) {
126    return internal::isApproxOrLessThan(a, b, test_precision<Real>());
127  }
128
129  Real get_test_precision(const Real&) {
130    return test_precision<Real>();
131  }
132
133  Real test_relative_error(const Real &a, const Real &b) {
134    using Eigen::numext::abs2;
135    return sqrt(abs2<Real>(a-b)/Eigen::numext::mini<Real>(abs2(a),abs2(b)));
136  }
137}
138}
139
140namespace Eigen {
141
142}
143
144void test_boostmultiprec()
145{
146  typedef Matrix<Real,Dynamic,Dynamic> Mat;
147  typedef Matrix<std::complex<Real>,Dynamic,Dynamic> MatC;
148
149  std::cout << "NumTraits<Real>::epsilon()         = " << NumTraits<Real>::epsilon() << std::endl;
150  std::cout << "NumTraits<Real>::dummy_precision() = " << NumTraits<Real>::dummy_precision() << std::endl;
151  std::cout << "NumTraits<Real>::lowest()          = " << NumTraits<Real>::lowest() << std::endl;
152  std::cout << "NumTraits<Real>::highest()         = " << NumTraits<Real>::highest() << std::endl;
153  std::cout << "NumTraits<Real>::digits10()        = " << NumTraits<Real>::digits10() << std::endl;
154
155  // chekc stream output
156  {
157    Mat A(10,10);
158    A.setRandom();
159    std::stringstream ss;
160    ss << A;
161  }
162  {
163    MatC A(10,10);
164    A.setRandom();
165    std::stringstream ss;
166    ss << A;
167  }
168
169  for(int i = 0; i < g_repeat; i++) {
170    int s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
171
172    CALL_SUBTEST_1( cholesky(Mat(s,s)) );
173
174    CALL_SUBTEST_2( lu_non_invertible<Mat>() );
175    CALL_SUBTEST_2( lu_invertible<Mat>() );
176    CALL_SUBTEST_2( lu_non_invertible<MatC>() );
177    CALL_SUBTEST_2( lu_invertible<MatC>() );
178
179    CALL_SUBTEST_3( qr(Mat(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
180    CALL_SUBTEST_3( qr_invertible<Mat>() );
181
182    CALL_SUBTEST_4( qr<Mat>() );
183    CALL_SUBTEST_4( cod<Mat>() );
184    CALL_SUBTEST_4( qr_invertible<Mat>() );
185
186    CALL_SUBTEST_5( qr<Mat>() );
187    CALL_SUBTEST_5( qr_invertible<Mat>() );
188
189    CALL_SUBTEST_6( selfadjointeigensolver(Mat(s,s)) );
190
191    CALL_SUBTEST_7( eigensolver(Mat(s,s)) );
192
193    CALL_SUBTEST_8( generalized_eigensolver_real(Mat(s,s)) );
194
195    TEST_SET_BUT_UNUSED_VARIABLE(s)
196  }
197
198  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))) ));
199  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))) ));
200}
201
202