diff --git a/config/Jamfile.v2 b/config/Jamfile.v2 index 55ead06d..271a3a5e 100644 --- a/config/Jamfile.v2 +++ b/config/Jamfile.v2 @@ -69,6 +69,7 @@ mp-run-simple has_tommath.cpp tommath : : : $(tommath_path) : has_tommath ; mp-run-simple has_float128.cpp quadmath : : : : has_float128 ; exe has_intel_quad : has_intel_quad.cpp : -Qoption,cpp,--extended_float_type ; +exe has_eigen : has_eigen.cpp ; explicit has_gmp ; explicit has_mpfr ; @@ -77,6 +78,6 @@ explicit has_tommath ; explicit has_float128 ; explicit has_intel_quad ; explicit has_mpc ; - +explicit has_eigen ; diff --git a/config/has_eigen.cpp b/config/has_eigen.cpp new file mode 100644 index 00000000..9412f3ee --- /dev/null +++ b/config/has_eigen.cpp @@ -0,0 +1,17 @@ +// Copyright John Maddock 2011. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + + +int main() +{ +#if EIGEN_VERSION_AT_LEAST(3,3,0) +#else +#error "Obsolete Eigen" +#endif + return 0; +} + diff --git a/include/boost/multiprecision/eigen.hpp b/include/boost/multiprecision/eigen.hpp new file mode 100644 index 00000000..c288507a --- /dev/null +++ b/include/boost/multiprecision/eigen.hpp @@ -0,0 +1,173 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MP_EIGEN_HPP +#define BOOST_MP_EIGEN_HPP + +#include +#include + +// +// Generic Eigen support code: +// +namespace Eigen +{ + template + struct NumTraits > + { + typedef boost::multiprecision::number self_type; + typedef typename boost::multiprecision::scalar_result_from_possible_complex::type Real; + typedef self_type NonInteger; // Not correct but we can't do much better?? + typedef double Literal; + typedef self_type Nested; + enum { + IsComplex = boost::multiprecision::number_category::value == boost::multiprecision::number_kind_complex, + IsInteger = boost::multiprecision::number_category::value == boost::multiprecision::number_kind_integer, + ReadCost = 1, + AddCost = 4, + MulCost = 8, + IsSigned = std::numeric_limits::is_specialized ? std::numeric_limits::is_signed : true, + RequireInitialization = 1, + }; + static Real epsilon() + { + return std::numeric_limits::epsilon(); + } + static Real dummy_precision() + { + return 1000 * epsilon(); + } + static Real highest() + { + return (std::numeric_limits::max)(); + } + static Real lowest() + { + return (std::numeric_limits::min)(); + } + static int digits10_imp(const boost::mpl::true_&) + { + return std::numeric_limits::digits10; + } + template + static int digits10_imp(const boost::mpl::bool_&) + { + return Real::default_precision(); + } + static int digits10() + { + return digits10_imp(boost::mpl::bool_::digits10 && (std::numeric_limits::digits10 != INT_MAX) ? true : false>()); + } + }; + template + struct NumTraits > : public NumTraits::result_type> + { + }; + +#define BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(A)\ + template\ + struct ScalarBinaryOpTraits, A, BinaryOp>\ + {\ + /*static_assert(boost::multiprecision::is_compatible_arithmetic_type >::value, "Interoperability with this arithmetic type is not supported.");*/\ + typedef boost::multiprecision::number ReturnType;\ + };\ + template\ + struct ScalarBinaryOpTraits, BinaryOp>\ + {\ + /*static_assert(boost::multiprecision::is_compatible_arithmetic_type >::value, "Interoperability with this arithmetic type is not supported.");*/\ + typedef boost::multiprecision::number ReturnType;\ + };\ + + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(float) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(double) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(long double) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(char) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned char) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(signed char) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(short) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned short) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(int) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned int) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(long) + BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned long) + +#if 0 + template + struct ScalarBinaryOpTraits, boost::multiprecision::number, BinaryOp> + { + static_assert( + boost::multiprecision::is_compatible_arithmetic_type, boost::multiprecision::number >::value + || boost::multiprecision::is_compatible_arithmetic_type, boost::multiprecision::number >::value, "Interoperability with this arithmetic type is not supported."); + typedef typename boost::mpl::if_c, boost::multiprecision::number >::value, + boost::multiprecision::number, boost::multiprecision::number >::type ReturnType; + }; + + template + struct ScalarBinaryOpTraits, boost::multiprecision::et_on>, boost::multiprecision::mpfr_float, BinaryOp> + { + typedef boost::multiprecision::number, boost::multiprecision::et_on> ReturnType; + }; + + template + struct ScalarBinaryOpTraits + { + typedef boost::multiprecision::number, boost::multiprecision::et_on> ReturnType; + }; + + template + struct ScalarBinaryOpTraits, boost::multiprecision::number, BinaryOp> + { + typedef boost::multiprecision::number ReturnType; + }; +#endif + + template + struct ScalarBinaryOpTraits, boost::multiprecision::detail::expression, BinaryOp> + { + static_assert(boost::is_convertible::result_type, boost::multiprecision::number >::value, "Interoperability with this arithmetic type is not supported."); + typedef boost::multiprecision::number ReturnType; + }; + + template + struct ScalarBinaryOpTraits, boost::multiprecision::number, BinaryOp> + { + static_assert(boost::is_convertible::result_type, boost::multiprecision::number >::value, "Interoperability with this arithmetic type is not supported."); + typedef boost::multiprecision::number ReturnType; + }; + + + namespace internal + { + template + struct conj_retval; + + template + struct conj_impl; + + + template + struct conj_retval > + { + typedef typename boost::multiprecision::detail::expression::result_type type; + }; + + template + struct conj_impl, true> + { + EIGEN_DEVICE_FUNC + static inline typename boost::multiprecision::detail::expression::result_type run(const typename boost::multiprecision::detail::expression& x) + { + return conj(x); + } + }; + + + + } + +} + + +#endif diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 416b8c94..ed7481bf 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -879,6 +879,21 @@ run test_hash.cpp : : : [ check-target-builds ../config//has_mpfi : TEST_MPFI gmp mpfr mpfi : ] [ check-target-builds ../config//has_tommath : TEST_TOMMATH tommath : ] ; +# +# Eigen interoperability: +# +run test_eigen_interop_cpp_int.cpp : : : release [ check-target-builds ../config//has_eigen : : no ] ; +run test_eigen_interop_cpp_dec_float.cpp : : : release [ check-target-builds ../config//has_eigen : : no ] ; +run test_eigen_interop_cpp_dec_float_2.cpp : : : release [ check-target-builds ../config//has_eigen : : no ] ; +run test_eigen_interop_cpp_dec_float_3.cpp : : : release [ check-target-builds ../config//has_eigen : : no ] ; +run test_eigen_interop_cpp_bin_float_1.cpp : : : release [ check-target-builds ../config//has_eigen : : no ] ; +run test_eigen_interop_cpp_bin_float_2.cpp : : : release [ check-target-builds ../config//has_eigen : : no ] ; +run test_eigen_interop_cpp_bin_float_3.cpp : : : release [ check-target-builds ../config//has_eigen : : no ] ; +run test_eigen_interop_mpfr_1.cpp mpfr gmp : : : release [ check-target-builds ../config//has_eigen : : no ] [ check-target-builds ../config//has_mpfr : : no ] ; +run test_eigen_interop_mpfr_2.cpp mpfr gmp : : : release [ check-target-builds ../config//has_eigen : : no ] [ check-target-builds ../config//has_mpfr : : no ] ; +run test_eigen_interop_mpfr_3.cpp mpfr gmp : : : release [ check-target-builds ../config//has_eigen : : no ] [ check-target-builds ../config//has_mpfr : : no ] ; +run test_eigen_interop_gmp.cpp gmp : : : release [ check-target-builds ../config//has_eigen : : no ] [ check-target-builds ../config//has_gmp : : no ] ; +run test_eigen_interop_mpc.cpp mpc mpfr gmp : : : release [ check-target-builds ../config//has_eigen : : no ] [ check-target-builds ../config//has_mpc : : no ] ; if $(enable-specfun) { diff --git a/test/eigen.hpp b/test/eigen.hpp new file mode 100644 index 00000000..9c442505 --- /dev/null +++ b/test/eigen.hpp @@ -0,0 +1,557 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include "test.hpp" + +using namespace Eigen; + +template +struct related_number +{ + typedef T type; +}; + +template +void example1() +{ + // expected results first: + Matrix r1, r2; + r1 << 3, 5, 4, 8; + r2 << -1, -1, 2, 0; + Matrix r3; + r3 << -1, -4, -6; + + + Matrix a; + a << 1, 2, 3, 4; + Matrix b(2, 2); + b << 2, 3, 1, 4; + std::cout << "a + b =\n" << a + b << std::endl; + BOOST_CHECK_EQUAL(a + b, r1); + std::cout << "a - b =\n" << a - b << std::endl; + BOOST_CHECK_EQUAL(a - b, r2); + std::cout << "Doing a += b;" << std::endl; + a += b; + std::cout << "Now a =\n" << a << std::endl; + Matrix v(1, 2, 3); + Matrix w(1, 0, 0); + std::cout << "-v + w - v =\n" << -v + w - v << std::endl; + BOOST_CHECK_EQUAL(-v + w - v, r3); +} + +template +void example2() +{ + Matrix a; + a << 1, 2, 3, 4; + Matrix v(1, 2, 3); + std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl; + std::cout << "0.1 * v =\n" << 0.1 * v << std::endl; + std::cout << "Doing v *= 2;" << std::endl; + v *= 2; + std::cout << "Now v =\n" << v << std::endl; + Num n(4); + std::cout << "Doing v *= Num;" << std::endl; + v *= n; + std::cout << "Now v =\n" << v << std::endl; + typedef typename related_number::type related_type; + related_type r(6); + std::cout << "Doing v *= RelatedType;" << std::endl; + v *= r; + std::cout << "Now v =\n" << v << std::endl; + std::cout << "RelatedType * v =\n" << r * v << std::endl; + std::cout << "Doing v *= RelatedType^2;" << std::endl; + v *= r * r; + std::cout << "Now v =\n" << v << std::endl; + std::cout << "RelatedType^2 * v =\n" << r * r * v << std::endl; +} + +template +void example3() +{ + using namespace std; + Matrix a = Matrix::Random(2, 2); + cout << "Here is the matrix a\n" << a << endl; + cout << "Here is the matrix a^T\n" << a.transpose() << endl; + cout << "Here is the conjugate of a\n" << a.conjugate() << endl; + cout << "Here is the matrix a^*\n" << a.adjoint() << endl; +} + +template +void example4() +{ + Matrix mat; + mat << 1, 2, + 3, 4; + Matrix u(-1, 1), v(2, 0); + std::cout << "Here is mat*mat:\n" << mat * mat << std::endl; + std::cout << "Here is mat*u:\n" << mat * u << std::endl; + std::cout << "Here is u^T*mat:\n" << u.transpose()*mat << std::endl; + std::cout << "Here is u^T*v:\n" << u.transpose()*v << std::endl; + std::cout << "Here is u*v^T:\n" << u * v.transpose() << std::endl; + std::cout << "Let's multiply mat by itself" << std::endl; + mat = mat * mat; + std::cout << "Now mat is mat:\n" << mat << std::endl; +} + +template +void example5() +{ + using namespace std; + Matrix v(1, 2, 3); + Matrix w(0, 1, 2); + cout << "Dot product: " << v.dot(w) << endl; + Num dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar + cout << "Dot product via a matrix product: " << dp << endl; + cout << "Cross product:\n" << v.cross(w) << endl; +} + +template +void example6() +{ + using namespace std; + Matrix mat; + mat << 1, 2, + 3, 4; + cout << "Here is mat.sum(): " << mat.sum() << endl; + cout << "Here is mat.prod(): " << mat.prod() << endl; + cout << "Here is mat.mean(): " << mat.mean() << endl; + cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; + cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; + cout << "Here is mat.trace(): " << mat.trace() << endl; +} + +template +void example7() +{ + using namespace std; + + Array m(2, 2); + + // assign some values coefficient by coefficient + m(0, 0) = 1.0; m(0, 1) = 2.0; + m(1, 0) = 3.0; m(1, 1) = m(0, 1) + m(1, 0); + + // print values to standard output + cout << m << endl << endl; + + // using the comma-initializer is also allowed + m << 1.0, 2.0, + 3.0, 4.0; + + // print values to standard output + cout << m << endl; +} + +template +void example8() +{ + using namespace std; + Array a(3, 3); + Array b(3, 3); + a << 1, 2, 3, + 4, 5, 6, + 7, 8, 9; + b << 1, 2, 3, + 1, 2, 3, + 1, 2, 3; + + // Adding two arrays + cout << "a + b = " << endl << a + b << endl << endl; + // Subtracting a scalar from an array + cout << "a - 2 = " << endl << a - 2 << endl; +} + +template +void example9() +{ + using namespace std; + Array a(2, 2); + Array b(2, 2); + a << 1, 2, + 3, 4; + b << 5, 6, + 7, 8; + cout << "a * b = " << endl << a * b << endl; +} + +template +void example10() +{ + using namespace std; + Array a = Array::Random(5); + a *= 2; + cout << "a =" << endl + << a << endl; + cout << "a.abs() =" << endl + << a.abs() << endl; + cout << "a.abs().sqrt() =" << endl + << a.abs().sqrt() << endl; + cout << "a.min(a.abs().sqrt()) =" << endl + << a.min(a.abs().sqrt()) << endl; +} + +template +void example11() +{ + using namespace std; + Matrix m(2, 2); + Matrix n(2, 2); + Matrix result(2, 2); + m << 1, 2, + 3, 4; + n << 5, 6, + 7, 8; + result = m * n; + cout << "-- Matrix m*n: --" << endl << result << endl << endl; + result = m.array() * n.array(); + cout << "-- Array m*n: --" << endl << result << endl << endl; + result = m.cwiseProduct(n); + cout << "-- With cwiseProduct: --" << endl << result << endl << endl; + result = m.array() + 4; + cout << "-- Array m + 4: --" << endl << result << endl << endl; +} + +template +void example12() +{ + using namespace std; + Matrix m(2, 2); + Matrix n(2, 2); + Matrix result(2, 2); + m << 1, 2, + 3, 4; + n << 5, 6, + 7, 8; + + result = (m.array() + 4).matrix() * m; + cout << "-- Combination 1: --" << endl << result << endl << endl; + result = (m.array() * n.array()).matrix() * m; + cout << "-- Combination 2: --" << endl << result << endl << endl; +} + +template +void example13() +{ + using namespace std; + Matrix m(4, 4); + m << 1, 2, 3, 4, + 5, 6, 7, 8, + 9, 10, 11, 12, + 13, 14, 15, 16; + cout << "Block in the middle" << endl; + cout << m.template block<2, 2>(1, 1) << endl << endl; + for (int i = 1; i <= 3; ++i) + { + cout << "Block of size " << i << "x" << i << endl; + cout << m.block(0, 0, i, i) << endl << endl; + } +} + +template +void example14() +{ + using namespace std; + Array m; + m << 1, 2, + 3, 4; + Array a = Array::Constant(0.6); + cout << "Here is the array a:" << endl << a << endl << endl; + a.template block<2, 2>(1, 1) = m; + cout << "Here is now a with m copied into its central 2x2 block:" << endl << a << endl << endl; + a.block(0, 0, 2, 3) = a.block(2, 1, 2, 3); + cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:" << endl << a << endl << endl; +} + +template +void example15() +{ + using namespace std; + Eigen::Matrix m(3, 3); + m << 1, 2, 3, + 4, 5, 6, + 7, 8, 9; + cout << "Here is the matrix m:" << endl << m << endl; + cout << "2nd Row: " << m.row(1) << endl; + m.col(2) += 3 * m.col(0); + cout << "After adding 3 times the first column into the third column, the matrix m is:\n"; + cout << m << endl; +} + +template +void example16() +{ + using namespace std; + Matrix m; + m << 1, 2, 3, 4, + 5, 6, 7, 8, + 9, 10, 11, 12, + 13, 14, 15, 16; + cout << "m.leftCols(2) =" << endl << m.leftCols(2) << endl << endl; + cout << "m.bottomRows<2>() =" << endl << m.template bottomRows<2>() << endl << endl; + m.topLeftCorner(1, 3) = m.bottomRightCorner(3, 1).transpose(); + cout << "After assignment, m = " << endl << m << endl; +} + +template +void example17() +{ + using namespace std; + Array v(6); + v << 1, 2, 3, 4, 5, 6; + cout << "v.head(3) =" << endl << v.head(3) << endl << endl; + cout << "v.tail<3>() = " << endl << v.template tail<3>() << endl << endl; + v.segment(1, 4) *= 2; + cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl; +} + +template +void example18() +{ + using namespace std; + Matrix mat; + mat << 1, 2, + 3, 4; + cout << "Here is mat.sum(): " << mat.sum() << endl; + cout << "Here is mat.prod(): " << mat.prod() << endl; + cout << "Here is mat.mean(): " << mat.mean() << endl; + cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; + cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; + cout << "Here is mat.trace(): " << mat.trace() << endl; + + BOOST_CHECK_EQUAL(mat.sum(), 10); + BOOST_CHECK_EQUAL(mat.prod(), 24); + BOOST_CHECK_EQUAL(mat.mean(), Num(5) / 2); + BOOST_CHECK_EQUAL(mat.minCoeff(), 1); + BOOST_CHECK_EQUAL(mat.maxCoeff(), 4); + BOOST_CHECK_EQUAL(mat.trace(), 5); +} + +template +void example18a() +{ + using namespace std; + Matrix mat; + mat << 1, 2, + 3, 4; + cout << "Here is mat.sum(): " << mat.sum() << endl; + cout << "Here is mat.prod(): " << mat.prod() << endl; + cout << "Here is mat.mean(): " << mat.mean() << endl; + //cout << "Here is mat.minCoeff(): " << mat.minCoeff() << endl; + //cout << "Here is mat.maxCoeff(): " << mat.maxCoeff() << endl; + cout << "Here is mat.trace(): " << mat.trace() << endl; +} + +template +void example19() +{ + using namespace std; + Matrix v(2); + Matrix m(2, 2), n(2, 2); + + v << -1, + 2; + + m << 1, -2, + -3, 4; + cout << "v.squaredNorm() = " << v.squaredNorm() << endl; + cout << "v.norm() = " << v.norm() << endl; + cout << "v.lpNorm<1>() = " << v.template lpNorm<1>() << endl; + cout << "v.lpNorm() = " << v.template lpNorm() << endl; + cout << endl; + cout << "m.squaredNorm() = " << m.squaredNorm() << endl; + cout << "m.norm() = " << m.norm() << endl; + cout << "m.lpNorm<1>() = " << m.template lpNorm<1>() << endl; + cout << "m.lpNorm() = " << m.template lpNorm() << endl; +} + +template +void example20() +{ + using namespace std; + Matrix A; + Matrix b; + A << 1, 2, 3, 4, 5, 6, 7, 8, 10; + b << 3, 3, 4; + cout << "Here is the matrix A:\n" << A << endl; + cout << "Here is the vector b:\n" << b << endl; + Matrix x = A.colPivHouseholderQr().solve(b); + cout << "The solution is:\n" << x << endl; +} + +template +void example21() +{ + using namespace std; + Matrix A, b; + A << 2, -1, -1, 3; + b << 1, 2, 3, 1; + cout << "Here is the matrix A:\n" << A << endl; + cout << "Here is the right hand side b:\n" << b << endl; + Matrix x = A.ldlt().solve(b); + cout << "The solution is:\n" << x << endl; +} + +template +void example22() +{ + using namespace std; + Matrix A = Matrix::Random(100, 100); + Matrix b = Matrix::Random(100, 50); + Matrix x = A.fullPivLu().solve(b); + Matrix axmb = A * x - b; + double relative_error = static_cast(abs(axmb.norm() / b.norm())); // norm() is L2 norm + cout << "norm1 = " << axmb.norm() << endl; + cout << "norm2 = " << b.norm() << endl; + cout << "The relative error is:\n" << relative_error << endl; +} + +template +void example23() +{ + using namespace std; + Matrix A; + A << 1, 2, 2, 3; + cout << "Here is the matrix A:\n" << A << endl; + SelfAdjointEigenSolver > eigensolver(A); + if (eigensolver.info() != Success) + { + std::cout << "Eigenvalue solver failed!" << endl; + } + else + { + cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl; + cout << "Here's a matrix whose columns are eigenvectors of A \n" + << "corresponding to these eigenvalues:\n" + << eigensolver.eigenvectors() << endl; + } +} + +template +void example24() +{ + using namespace std; + Matrix A; + A << 1, 2, 1, + 2, 1, 0, + -1, 1, 2; + cout << "Here is the matrix A:\n" << A << endl; + cout << "The determinant of A is " << A.determinant() << endl; + cout << "The inverse of A is:\n" << A.inverse() << endl; +} + +template +void test_integer_type() +{ + example1(); + //example2(); + example18(); +} + +template +void test_float_type() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example1(); + example2(); + example4(); + example5(); + example6(); + example7(); + example8(); + example9(); + example10(); + example11(); + example12(); + example13(); + example14(); + example15(); + example16(); + example17(); + /* + example18(); + example19(); + example20(); + example21(); + example22(); + example23(); + example24(); + */ +} + +template +void test_float_type_2() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example18(); + example19(); + example20(); + example21(); + + //example22(); + //example23(); + //example24(); +} + +template +void test_float_type_3() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example22(); + example23(); + example24(); +} + +template +void test_complex_type() +{ + std::cout << "Epsilon = " << Eigen::NumTraits::epsilon() << std::endl; + std::cout << "Dummy Prec = " << Eigen::NumTraits::dummy_precision() << std::endl; + std::cout << "Highest = " << Eigen::NumTraits::highest() << std::endl; + std::cout << "Lowest = " << Eigen::NumTraits::lowest() << std::endl; + std::cout << "Digits10 = " << Eigen::NumTraits::digits10() << std::endl; + + example1(); + example2(); + example3(); + example4(); + example5(); + example7(); + example8(); + example9(); + example11(); + example12(); + example13(); + example14(); + example15(); + example16(); + example17(); + example18a(); + example19(); + example20(); + example21(); + example22(); + // example23(); //requires comparisons. + example24(); +} + diff --git a/test/test_eigen_interop.cpp b/test/test_eigen_interop.cpp index 3a2fb344..a675d9c3 100644 --- a/test/test_eigen_interop.cpp +++ b/test/test_eigen_interop.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include @@ -90,7 +91,7 @@ namespace Eigen BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned int) BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(long) BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned long) - +#if 0 template struct ScalarBinaryOpTraits, boost::multiprecision::number, BinaryOp> { @@ -101,16 +102,16 @@ namespace Eigen boost::multiprecision::number, boost::multiprecision::number >::type ReturnType; }; - template - struct ScalarBinaryOpTraits, typename NumTraits >::IsComplex, boost::multiprecision::number >::type>::Real, BinaryOp> + template + struct ScalarBinaryOpTraits, boost::multiprecision::et_on>, boost::multiprecision::mpfr_float, BinaryOp> { - typedef boost::multiprecision::number ReturnType; + typedef boost::multiprecision::number, boost::multiprecision::et_on> ReturnType; }; - template - struct ScalarBinaryOpTraits >::IsComplex, boost::multiprecision::number >::type>::Real, boost::multiprecision::number, BinaryOp> + template + struct ScalarBinaryOpTraits { - typedef boost::multiprecision::number ReturnType; + typedef boost::multiprecision::number, boost::multiprecision::et_on> ReturnType; }; template @@ -118,7 +119,7 @@ namespace Eigen { typedef boost::multiprecision::number ReturnType; }; - +#endif template struct ScalarBinaryOpTraits, boost::multiprecision::detail::expression, BinaryOp> { @@ -135,11 +136,13 @@ namespace Eigen } + template struct related_number { typedef T type; }; +/* template <> struct related_number { @@ -149,7 +152,7 @@ template <> struct related_number { typedef boost::multiprecision::cpp_dec_float_50 type; -}; +};*/ template void example1() @@ -541,10 +544,47 @@ void example22() Matrix A = Matrix::Random(100, 100); Matrix b = Matrix::Random(100, 50); Matrix x = A.fullPivLu().solve(b); - double relative_error = static_cast(abs((A*x - b).norm() / b.norm())); // norm() is L2 norm + Matrix axmb = A * x - b; + double relative_error = static_cast(abs(axmb.norm() / b.norm())); // norm() is L2 norm + cout << "norm1 = " << axmb.norm() << endl; + cout << "norm2 = " << b.norm() << endl; cout << "The relative error is:\n" << relative_error << endl; } +template +void example23() +{ + using namespace std; + Matrix A; + A << 1, 2, 2, 3; + cout << "Here is the matrix A:\n" << A << endl; + SelfAdjointEigenSolver > eigensolver(A); + if (eigensolver.info() != Success) + { + std::cout << "Eigenvalue solver failed!" << endl; + } + else + { + cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl; + cout << "Here's a matrix whose columns are eigenvectors of A \n" + << "corresponding to these eigenvalues:\n" + << eigensolver.eigenvectors() << endl; + } +} + +template +void example24() +{ + using namespace std; + Matrix A; + A << 1, 2, 1, + 2, 1, 0, + -1, 1, 2; + cout << "Here is the matrix A:\n" << A << endl; + cout << "The determinant of A is " << A.determinant() << endl; + cout << "The inverse of A is:\n" << A.inverse() << endl; +} + template void test_integer_type() { @@ -583,6 +623,8 @@ void test_float_type() example20(); example21(); example22(); + example23(); + example24(); } template @@ -614,13 +656,29 @@ void test_complex_type() example20(); example21(); example22(); + // example23(); //requires comparisons. + example24(); +} + +namespace boost { + namespace multiprecision { + + template + inline void log_postfix_event(const mpc_complex_backend& val, const char* event_description) + { + if (mpfr_nan_p(mpc_realref(val.data()))) + { + std::cout << "Found a NaN! " << event_description << std::endl; + } + } + + } } int main() { - test_float_type(); - test_float_type(); - test_complex_type(); + using namespace boost::multiprecision; + test_complex_type(); #if 0 test_integer_type(); @@ -628,6 +686,8 @@ int main() test_complex_type >(); test_float_type(); + test_float_type(); + test_float_type(); test_integer_type(); test_integer_type(); diff --git a/test/test_eigen_interop_cpp_bin_float_1.cpp b/test/test_eigen_interop_cpp_bin_float_1.cpp new file mode 100644 index 00000000..80cedcdb --- /dev/null +++ b/test/test_eigen_interop_cpp_bin_float_1.cpp @@ -0,0 +1,16 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_float_type(); + test_float_type(); + return 0; +} diff --git a/test/test_eigen_interop_cpp_bin_float_2.cpp b/test/test_eigen_interop_cpp_bin_float_2.cpp new file mode 100644 index 00000000..13c576af --- /dev/null +++ b/test/test_eigen_interop_cpp_bin_float_2.cpp @@ -0,0 +1,16 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_float_type_2(); + test_float_type_2(); + return 0; +} diff --git a/test/test_eigen_interop_cpp_bin_float_3.cpp b/test/test_eigen_interop_cpp_bin_float_3.cpp new file mode 100644 index 00000000..67fce053 --- /dev/null +++ b/test/test_eigen_interop_cpp_bin_float_3.cpp @@ -0,0 +1,16 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_float_type_3(); + test_float_type_3(); + return 0; +} diff --git a/test/test_eigen_interop_cpp_dec_float.cpp b/test/test_eigen_interop_cpp_dec_float.cpp new file mode 100644 index 00000000..8723ce31 --- /dev/null +++ b/test/test_eigen_interop_cpp_dec_float.cpp @@ -0,0 +1,23 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +template <> +struct related_number +{ + typedef boost::multiprecision::cpp_dec_float_50 type; +}; + + +int main() +{ + using namespace boost::multiprecision; + test_float_type(); + test_float_type(); + return 0; +} diff --git a/test/test_eigen_interop_cpp_dec_float_2.cpp b/test/test_eigen_interop_cpp_dec_float_2.cpp new file mode 100644 index 00000000..89ef555e --- /dev/null +++ b/test/test_eigen_interop_cpp_dec_float_2.cpp @@ -0,0 +1,23 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +template <> +struct related_number +{ + typedef boost::multiprecision::cpp_dec_float_50 type; +}; + + +int main() +{ + using namespace boost::multiprecision; + //test_float_type_2(); + test_float_type_2(); + return 0; +} diff --git a/test/test_eigen_interop_cpp_dec_float_3.cpp b/test/test_eigen_interop_cpp_dec_float_3.cpp new file mode 100644 index 00000000..76a99c90 --- /dev/null +++ b/test/test_eigen_interop_cpp_dec_float_3.cpp @@ -0,0 +1,23 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +template <> +struct related_number +{ + typedef boost::multiprecision::cpp_dec_float_50 type; +}; + + +int main() +{ + using namespace boost::multiprecision; + //test_float_type_2(); + test_float_type_3(); + return 0; +} diff --git a/test/test_eigen_interop_cpp_int.cpp b/test/test_eigen_interop_cpp_int.cpp new file mode 100644 index 00000000..fba52e06 --- /dev/null +++ b/test/test_eigen_interop_cpp_int.cpp @@ -0,0 +1,17 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_integer_type(); + test_integer_type(); + test_integer_type(); + return 0; +} diff --git a/test/test_eigen_interop_gmp.cpp b/test/test_eigen_interop_gmp.cpp new file mode 100644 index 00000000..4cb2385b --- /dev/null +++ b/test/test_eigen_interop_gmp.cpp @@ -0,0 +1,16 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_integer_type(); + test_integer_type(); + return 0; +} diff --git a/test/test_eigen_interop_mpc.cpp b/test/test_eigen_interop_mpc.cpp new file mode 100644 index 00000000..6ac9966c --- /dev/null +++ b/test/test_eigen_interop_mpc.cpp @@ -0,0 +1,15 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_complex_type(); + return 0; +} diff --git a/test/test_eigen_interop_mpfr_1.cpp b/test/test_eigen_interop_mpfr_1.cpp new file mode 100644 index 00000000..7ddbbcc5 --- /dev/null +++ b/test/test_eigen_interop_mpfr_1.cpp @@ -0,0 +1,15 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_float_type(); + return 0; +} diff --git a/test/test_eigen_interop_mpfr_2.cpp b/test/test_eigen_interop_mpfr_2.cpp new file mode 100644 index 00000000..27cd9ddb --- /dev/null +++ b/test/test_eigen_interop_mpfr_2.cpp @@ -0,0 +1,15 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_float_type_2(); + return 0; +} diff --git a/test/test_eigen_interop_mpfr_3.cpp b/test/test_eigen_interop_mpfr_3.cpp new file mode 100644 index 00000000..d4541098 --- /dev/null +++ b/test/test_eigen_interop_mpfr_3.cpp @@ -0,0 +1,15 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2018 John Maddock. Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include + +#include "eigen.hpp" + +int main() +{ + using namespace boost::multiprecision; + test_float_type_3(); + return 0; +}