1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2016 Gael Guennebaud <[email protected]>
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> {
dummy_precisionEigen::NumTraits82 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<>
test_precision()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> {
runEigen::internal::cast_impl95 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> > {
runEigen::internal::cast_impl102 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>:
fabs(const Real & a)119 Real fabs(const Real& a) { return abs EIGEN_NOT_A_MACRO (a); }
fmax(const Real & a,const Real & b)120 Real fmax(const Real& a, const Real& b) { using std::max; return max(a,b); }
121
122 // some specialization for the unit tests:
test_isMuchSmallerThan(const Real & a,const Real & b)123 inline bool test_isMuchSmallerThan(const Real& a, const Real& b) {
124 return internal::isMuchSmallerThan(a, b, test_precision<Real>());
125 }
126
test_isApprox(const Real & a,const Real & b)127 inline bool test_isApprox(const Real& a, const Real& b) {
128 return internal::isApprox(a, b, test_precision<Real>());
129 }
130
test_isApproxOrLessThan(const Real & a,const Real & b)131 inline bool test_isApproxOrLessThan(const Real& a, const Real& b) {
132 return internal::isApproxOrLessThan(a, b, test_precision<Real>());
133 }
134
get_test_precision(const Real &)135 Real get_test_precision(const Real&) {
136 return test_precision<Real>();
137 }
138
test_relative_error(const Real & a,const Real & b)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
EIGEN_DECLARE_TEST(boostmultiprec)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 }
209