diff options
| author |  Benoit Steiner <benoit.steiner.goog@gmail.com> | 2016-10-05 18:48:55 -0700 |
|---|---|---|
| committer | Benoit Steiner <benoit.steiner.goog@gmail.com> | 2016-10-05 18:48:55 -0700 |
| commit | 78b569f68540c5609388864bd805dcf21dd6a187 (patch) | |
| tree | 0a5757bb11834d0109f99310f4493dfd63579901 /test | |
| parent | 9c2b6c049be19fd4c571b0df537169d277b26291 (diff) | |
| parent | 4387433acf9cd2eab3713349163cd1e8905b5854 (diff) | |
Merged latest updates from trunk
Diffstat (limited to 'test')
64 files changed, 1783 insertions, 342 deletions
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7bed6a45c..e17985107 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -147,7 +147,7 @@ ei_add_test(nomalloc) ei_add_test(first_aligned) ei_add_test(nullary) ei_add_test(mixingtypes) -ei_add_test(packetmath) +ei_add_test(packetmath "-DEIGEN_FAST_MATH=1") ei_add_test(unalignedassert) ei_add_test(vectorization_logic) ei_add_test(basicstuff) @@ -258,6 +258,11 @@ ei_add_test(rvalue_types) ei_add_test(dense_storage) ei_add_test(ctorleak) ei_add_test(mpl2only) +ei_add_test(inplace_decomposition) +ei_add_test(half_float) +ei_add_test(array_of_string) + +add_executable(bug1213 bug1213.cpp bug1213_main.cpp) check_cxx_compiler_flag("-ffast-math" COMPILER_SUPPORT_FASTMATH) if(COMPILER_SUPPORT_FASTMATH) @@ -324,6 +329,16 @@ if(EIGEN_TEST_EIGEN2) message(WARNING "The Eigen2 test suite has been removed") endif() +# boost MP unit test +find_package(Boost) +if(Boost_FOUND) + include_directories(${Boost_INCLUDE_DIRS}) + ei_add_test(boostmultiprec "" "${Boost_LIBRARIES}") + ei_add_property(EIGEN_TESTED_BACKENDS "Boost.Multiprecision, ") +else() + ei_add_property(EIGEN_MISSING_BACKENDS "Boost.Multiprecision, ") +endif() + # CUDA unit tests option(EIGEN_TEST_CUDA "Enable CUDA support in unit tests" OFF) @@ -340,7 +355,7 @@ if(CUDA_FOUND) set(CUDA_PROPAGATE_HOST_FLAGS OFF) if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") - set(CUDA_NVCC_FLAGS "-ccbin /usr/bin/clang" CACHE STRING "nvcc flags" FORCE) + set(CUDA_NVCC_FLAGS "-ccbin ${CMAKE_C_COMPILER}" CACHE STRING "nvcc flags" FORCE) endif() if(EIGEN_TEST_CUDA_CLANG) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 --cuda-gpu-arch=sm_30") diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 9c895e0ac..bdea51c10 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -169,7 +169,7 @@ void test_adjoint() // test a large static matrix only once CALL_SUBTEST_7( adjoint(Matrix<float, 100, 100>()) ); -#ifdef EIGEN_TEST_PART_4 +#ifdef EIGEN_TEST_PART_13 { MatrixXcf a(10,10), b(10,10); VERIFY_RAISES_ASSERT(a = a.transpose()); @@ -187,6 +187,13 @@ void test_adjoint() a.transpose() = a.adjoint(); a.transpose() += a.adjoint(); a.transpose() += a.adjoint() + b; + + // regression tests for check_for_aliasing + MatrixXd c(10,10); + c = 1.0 * MatrixXd::Ones(10,10) + c; + c = MatrixXd::Ones(10,10) * 1.0 + c; + c = c + MatrixXd::Ones(10,10) .cwiseProduct( MatrixXd::Zero(10,10) ); + c = MatrixXd::Ones(10,10) * MatrixXd::Zero(10,10); } #endif } diff --git a/test/array.cpp b/test/array.cpp index beaa62221..15c3266a9 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -13,6 +13,7 @@ template<typename ArrayType> void array(const ArrayType& m) { typedef typename ArrayType::Index Index; typedef typename ArrayType::Scalar Scalar; + typedef typename ArrayType::RealScalar RealScalar; typedef Array<Scalar, ArrayType::RowsAtCompileTime, 1> ColVectorType; typedef Array<Scalar, 1, ArrayType::ColsAtCompileTime> RowVectorType; @@ -72,7 +73,7 @@ template<typename ArrayType> void array(const ArrayType& m) VERIFY_IS_MUCH_SMALLER_THAN(abs(m1.rowwise().sum().sum() - m1.sum()), m1.abs().sum()); if (!internal::isMuchSmallerThan(abs(m1.sum() - (m1+m2).sum()), m1.abs().sum(), test_precision<Scalar>())) VERIFY_IS_NOT_APPROX(((m1+m2).rowwise().sum()).sum(), m1.sum()); - VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar>())); + VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar,Scalar>())); // vector-wise ops m3 = m1; @@ -102,6 +103,22 @@ template<typename ArrayType> void array(const ArrayType& m) FixedArrayType f4(f1.data()); VERIFY_IS_APPROX(f4, f1); + // pow + VERIFY_IS_APPROX(m1.pow(2), m1.square()); + VERIFY_IS_APPROX(pow(m1,2), m1.square()); + VERIFY_IS_APPROX(m1.pow(3), m1.cube()); + VERIFY_IS_APPROX(pow(m1,3), m1.cube()); + VERIFY_IS_APPROX((-m1).pow(3), -m1.cube()); + VERIFY_IS_APPROX(pow(2*m1,3), 8*m1.cube()); + ArrayType exponents = ArrayType::Constant(rows, cols, RealScalar(2)); + VERIFY_IS_APPROX(Eigen::pow(m1,exponents), m1.square()); + VERIFY_IS_APPROX(m1.pow(exponents), m1.square()); + VERIFY_IS_APPROX(Eigen::pow(2*m1,exponents), 4*m1.square()); + VERIFY_IS_APPROX((2*m1).pow(exponents), 4*m1.square()); + VERIFY_IS_APPROX(Eigen::pow(m1,2*exponents), m1.square().square()); + VERIFY_IS_APPROX(m1.pow(2*exponents), m1.square().square()); + VERIFY_IS_APPROX(Eigen::pow(m1(0,0), exponents), ArrayType::Constant(rows,cols,m1(0,0)*m1(0,0))); + // Check possible conflicts with 1D ctor typedef Array<Scalar, Dynamic, 1> OneDArrayType; OneDArrayType o1(rows); @@ -217,12 +234,7 @@ template<typename ArrayType> void array_real(const ArrayType& m) VERIFY_IS_APPROX(m1.sinh(), sinh(m1)); VERIFY_IS_APPROX(m1.cosh(), cosh(m1)); VERIFY_IS_APPROX(m1.tanh(), tanh(m1)); -#ifdef EIGEN_HAS_C99_MATH - VERIFY_IS_APPROX(m1.lgamma(), lgamma(m1)); - VERIFY_IS_APPROX(m1.digamma(), digamma(m1)); - VERIFY_IS_APPROX(m1.erf(), erf(m1)); - VERIFY_IS_APPROX(m1.erfc(), erfc(m1)); -#endif // EIGEN_HAS_C99_MATH + VERIFY_IS_APPROX(m1.arg(), arg(m1)); VERIFY_IS_APPROX(m1.round(), round(m1)); VERIFY_IS_APPROX(m1.floor(), floor(m1)); @@ -243,7 +255,9 @@ template<typename ArrayType> void array_real(const ArrayType& m) m3 = m1.abs(); VERIFY_IS_APPROX(m3.sqrt(), sqrt(abs(m1))); VERIFY_IS_APPROX(m3.rsqrt(), Scalar(1)/sqrt(abs(m1))); + VERIFY_IS_APPROX(rsqrt(m3), Scalar(1)/sqrt(abs(m1))); VERIFY_IS_APPROX(m3.log(), log(m3)); + VERIFY_IS_APPROX(m3.log1p(), log1p(m3)); VERIFY_IS_APPROX(m3.log10(), log10(m3)); @@ -275,27 +289,12 @@ template<typename ArrayType> void array_real(const ArrayType& m) // shift argument of logarithm so that it is not zero Scalar smallNumber = NumTraits<Scalar>::dummy_precision(); VERIFY_IS_APPROX((m3 + smallNumber).log() , log(abs(m1) + smallNumber)); + VERIFY_IS_APPROX((m3 + smallNumber + 1).log() , log1p(abs(m1) + smallNumber)); VERIFY_IS_APPROX(m1.exp() * m2.exp(), exp(m1+m2)); VERIFY_IS_APPROX(m1.exp(), exp(m1)); VERIFY_IS_APPROX(m1.exp() / m2.exp(),(m1-m2).exp()); - VERIFY_IS_APPROX(m1.pow(2), m1.square()); - VERIFY_IS_APPROX(pow(m1,2), m1.square()); - VERIFY_IS_APPROX(m1.pow(3), m1.cube()); - VERIFY_IS_APPROX(pow(m1,3), m1.cube()); - VERIFY_IS_APPROX((-m1).pow(3), -m1.cube()); - VERIFY_IS_APPROX(pow(2*m1,3), 8*m1.cube()); - - ArrayType exponents = ArrayType::Constant(rows, cols, RealScalar(2)); - VERIFY_IS_APPROX(Eigen::pow(m1,exponents), m1.square()); - VERIFY_IS_APPROX(m1.pow(exponents), m1.square()); - VERIFY_IS_APPROX(Eigen::pow(2*m1,exponents), 4*m1.square()); - VERIFY_IS_APPROX((2*m1).pow(exponents), 4*m1.square()); - VERIFY_IS_APPROX(Eigen::pow(m1,2*exponents), m1.square().square()); - VERIFY_IS_APPROX(m1.pow(2*exponents), m1.square().square()); - VERIFY_IS_APPROX(pow(m1(0,0), exponents), ArrayType::Constant(rows,cols,m1(0,0)*m1(0,0))); - VERIFY_IS_APPROX(m3.pow(RealScalar(0.5)), m3.sqrt()); VERIFY_IS_APPROX(pow(m3,RealScalar(0.5)), m3.sqrt()); @@ -310,122 +309,6 @@ template<typename ArrayType> void array_real(const ArrayType& m) m1 += ArrayType::Constant(rows,cols,Scalar(tiny)); VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse()); -#ifdef EIGEN_HAS_C99_MATH - // check special functions (comparing against numpy implementation) - if (!NumTraits<Scalar>::IsComplex) { - VERIFY_IS_APPROX(numext::digamma(Scalar(1)), RealScalar(-0.5772156649015329)); - VERIFY_IS_APPROX(numext::digamma(Scalar(1.5)), RealScalar(0.03648997397857645)); - VERIFY_IS_APPROX(numext::digamma(Scalar(4)), RealScalar(1.2561176684318)); - VERIFY_IS_APPROX(numext::digamma(Scalar(-10.5)), RealScalar(2.398239129535781)); - VERIFY_IS_APPROX(numext::digamma(Scalar(10000.5)), RealScalar(9.210340372392849)); - VERIFY_IS_EQUAL(numext::digamma(Scalar(0)), - std::numeric_limits<RealScalar>::infinity()); - VERIFY_IS_EQUAL(numext::digamma(Scalar(-1)), - std::numeric_limits<RealScalar>::infinity()); - - // Check the zeta function against scipy.special.zeta - VERIFY_IS_APPROX(numext::zeta(Scalar(1.5), Scalar(2)), RealScalar(1.61237534869)); - VERIFY_IS_APPROX(numext::zeta(Scalar(4), Scalar(1.5)), RealScalar(0.234848505667)); - VERIFY_IS_APPROX(numext::zeta(Scalar(10.5), Scalar(3)), RealScalar(1.03086757337e-5)); - VERIFY_IS_APPROX(numext::zeta(Scalar(10000.5), Scalar(1.0001)), RealScalar(0.367879440865)); - VERIFY_IS_APPROX(numext::zeta(Scalar(3), Scalar(-2.5)), RealScalar(0.054102025820864097)); - VERIFY_IS_EQUAL(numext::zeta(Scalar(1), Scalar(1.2345)), // The second scalar does not matter - std::numeric_limits<RealScalar>::infinity()); - VERIFY((numext::isnan)(numext::zeta(Scalar(0.9), Scalar(1.2345)))); // The second scalar does not matter - - // Check the polygamma against scipy.special.polygamma examples - VERIFY_IS_APPROX(numext::polygamma(Scalar(1), Scalar(2)), RealScalar(0.644934066848)); - VERIFY_IS_APPROX(numext::polygamma(Scalar(1), Scalar(3)), RealScalar(0.394934066848)); - VERIFY_IS_APPROX(numext::polygamma(Scalar(1), Scalar(25.5)), RealScalar(0.0399946696496)); - VERIFY((numext::isnan)(numext::polygamma(Scalar(1.5), Scalar(1.2345)))); // The second scalar does not matter - - // Check the polygamma function over a larger range of values - VERIFY_IS_APPROX(numext::polygamma(Scalar(17), Scalar(4.7)), RealScalar(293.334565435)); - VERIFY_IS_APPROX(numext::polygamma(Scalar(31), Scalar(11.8)), RealScalar(0.445487887616)); - VERIFY_IS_APPROX(numext::polygamma(Scalar(28), Scalar(17.7)), RealScalar(-2.47810300902e-07)); - VERIFY_IS_APPROX(numext::polygamma(Scalar(8), Scalar(30.2)), RealScalar(-8.29668781082e-09)); - /* The following tests only pass for doubles because floats cannot handle the large values of - the gamma function. - VERIFY_IS_APPROX(numext::polygamma(Scalar(42), Scalar(15.8)), RealScalar(-0.434562276666)); - VERIFY_IS_APPROX(numext::polygamma(Scalar(147), Scalar(54.1)), RealScalar(0.567742190178)); - VERIFY_IS_APPROX(numext::polygamma(Scalar(170), Scalar(64)), RealScalar(-0.0108615497927)); - */ - - { - // Test various propreties of igamma & igammac. These are normalized - // gamma integrals where - // igammac(a, x) = Gamma(a, x) / Gamma(a) - // igamma(a, x) = gamma(a, x) / Gamma(a) - // where Gamma and gamma are considered the standard unnormalized - // upper and lower incomplete gamma functions, respectively. - ArrayType a = m1.abs() + 2; - ArrayType x = m2.abs() + 2; - ArrayType zero = ArrayType::Zero(rows, cols); - ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0)); - ArrayType a_m1 = a - one; - ArrayType Gamma_a_x = Eigen::igammac(a, x) * a.lgamma().exp(); - ArrayType Gamma_a_m1_x = Eigen::igammac(a_m1, x) * a_m1.lgamma().exp(); - ArrayType gamma_a_x = Eigen::igamma(a, x) * a.lgamma().exp(); - ArrayType gamma_a_m1_x = Eigen::igamma(a_m1, x) * a_m1.lgamma().exp(); - - // Gamma(a, 0) == Gamma(a) - VERIFY_IS_APPROX(Eigen::igammac(a, zero), one); - - // Gamma(a, x) + gamma(a, x) == Gamma(a) - VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp()); - - // Gamma(a, x) == (a - 1) * Gamma(a-1, x) + x^(a-1) * exp(-x) - VERIFY_IS_APPROX(Gamma_a_x, (a - 1) * Gamma_a_m1_x + x.pow(a-1) * (-x).exp()); - - // gamma(a, x) == (a - 1) * gamma(a-1, x) - x^(a-1) * exp(-x) - VERIFY_IS_APPROX(gamma_a_x, (a - 1) * gamma_a_m1_x - x.pow(a-1) * (-x).exp()); - } - - // Check exact values of igamma and igammac against a third party calculation. - Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; - Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; - - // location i*6+j corresponds to a_s[i], x_s[j]. - Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); - Scalar igamma_s[][6] = {{0.0, nan, nan, nan, nan, nan}, - {0.0, 0.6321205588285578, 0.7768698398515702, - 0.9816843611112658, 9.999500016666262e-05, 1.0}, - {0.0, 0.4275932955291202, 0.608374823728911, - 0.9539882943107686, 7.522076445089201e-07, 1.0}, - {0.0, 0.01898815687615381, 0.06564245437845008, - 0.5665298796332909, 4.166333347221828e-18, 1.0}, - {0.0, 0.9999780593618628, 0.9999899967080838, - 0.9999996219837988, 0.9991370418689945, 1.0}, - {0.0, 0.0, 0.0, 0.0, 0.0, 0.5042041932513908}}; - Scalar igammac_s[][6] = {{nan, nan, nan, nan, nan, nan}, - {1.0, 0.36787944117144233, 0.22313016014842982, - 0.018315638888734182, 0.9999000049998333, 0.0}, - {1.0, 0.5724067044708798, 0.3916251762710878, - 0.04601170568923136, 0.9999992477923555, 0.0}, - {1.0, 0.9810118431238462, 0.9343575456215499, - 0.4334701203667089, 1.0, 0.0}, - {1.0, 2.1940638138146658e-05, 1.0003291916285e-05, - 3.7801620118431334e-07, 0.0008629581310054535, - 0.0}, - {1.0, 1.0, 1.0, 1.0, 1.0, 0.49579580674813944}}; - for (int i = 0; i < 6; ++i) { - for (int j = 0; j < 6; ++j) { - if ((std::isnan)(igamma_s[i][j])) { - VERIFY((std::isnan)(numext::igamma(a_s[i], x_s[j]))); - } else { - VERIFY_IS_APPROX(numext::igamma(a_s[i], x_s[j]), igamma_s[i][j]); - } - - if ((std::isnan)(igammac_s[i][j])) { - VERIFY((std::isnan)(numext::igammac(a_s[i], x_s[j]))); - } else { - VERIFY_IS_APPROX(numext::igammac(a_s[i], x_s[j]), igammac_s[i][j]); - } - } - } - } -#endif // EIGEN_HAS_C99_MATH - // check inplace transpose m3 = m1; m3.transposeInPlace(); @@ -525,7 +408,7 @@ template<typename ArrayType> void array_complex(const ArrayType& m) // scalar by array division Scalar s1 = internal::random<Scalar>(); - const RealScalar tiny = sqrt(std::numeric_limits<RealScalar>::epsilon()); + const RealScalar tiny = std::sqrt(std::numeric_limits<RealScalar>::epsilon()); s1 += Scalar(tiny); m1 += ArrayType::Constant(rows,cols,Scalar(tiny)); VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse()); @@ -605,7 +488,7 @@ void test_array() VERIFY((internal::is_same< internal::global_math_functions_filtering_base<int>::type, int >::value)); VERIFY((internal::is_same< internal::global_math_functions_filtering_base<float>::type, float >::value)); VERIFY((internal::is_same< internal::global_math_functions_filtering_base<Array2i>::type, ArrayBase<Array2i> >::value)); - typedef CwiseUnaryOp<internal::scalar_multiple_op<double>, ArrayXd > Xpr; + typedef CwiseUnaryOp<internal::scalar_abs_op<double>, ArrayXd > Xpr; VERIFY((internal::is_same< internal::global_math_functions_filtering_base<Xpr>::type, ArrayBase<Xpr> >::value)); diff --git a/test/array_for_matrix.cpp b/test/array_for_matrix.cpp index db5f3b34a..97e03be83 100644 --- a/test/array_for_matrix.cpp +++ b/test/array_for_matrix.cpp @@ -45,7 +45,7 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m) VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum().sum() - m1.sum(), m1.squaredNorm()); VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum() + m2.colwise().sum() - (m1+m2).colwise().sum(), (m1+m2).squaredNorm()); VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum() - m2.rowwise().sum() - (m1-m2).rowwise().sum(), (m1-m2).squaredNorm()); - VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar>())); + VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar,Scalar>())); // vector-wise ops m3 = m1; @@ -144,9 +144,21 @@ template<typename MatrixType> void comparisons(const MatrixType& m) template<typename VectorType> void lpNorm(const VectorType& v) { using std::sqrt; + typedef typename VectorType::RealScalar RealScalar; VectorType u = VectorType::Random(v.size()); - VERIFY_IS_APPROX(u.template lpNorm<Infinity>(), u.cwiseAbs().maxCoeff()); + if(v.size()==0) + { + VERIFY_IS_APPROX(u.template lpNorm<Infinity>(), RealScalar(0)); + VERIFY_IS_APPROX(u.template lpNorm<1>(), RealScalar(0)); + VERIFY_IS_APPROX(u.template lpNorm<2>(), RealScalar(0)); + VERIFY_IS_APPROX(u.template lpNorm<5>(), RealScalar(0)); + } + else + { + VERIFY_IS_APPROX(u.template lpNorm<Infinity>(), u.cwiseAbs().maxCoeff()); + } + VERIFY_IS_APPROX(u.template lpNorm<1>(), u.cwiseAbs().sum()); VERIFY_IS_APPROX(u.template lpNorm<2>(), sqrt(u.array().abs().square().sum())); VERIFY_IS_APPROX(numext::pow(u.template lpNorm<5>(), typename VectorType::RealScalar(5)), u.array().abs().pow(5).sum()); @@ -255,6 +267,8 @@ void test_array_for_matrix() CALL_SUBTEST_5( lpNorm(VectorXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); CALL_SUBTEST_4( lpNorm(VectorXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); } + CALL_SUBTEST_5( lpNorm(VectorXf(0)) ); + CALL_SUBTEST_4( lpNorm(VectorXcf(0)) ); for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_4( resize(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); CALL_SUBTEST_5( resize(MatrixXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); diff --git a/test/array_of_string.cpp b/test/array_of_string.cpp new file mode 100644 index 000000000..e23b7c59e --- /dev/null +++ b/test/array_of_string.cpp @@ -0,0 +1,32 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" + +void test_array_of_string() +{ + typedef Array<std::string,1,Dynamic> ArrayXs; + ArrayXs a1(3), a2(3), a3(3), a3ref(3); + a1 << "one", "two", "three"; + a2 << "1", "2", "3"; + a3ref << "one (1)", "two (2)", "three (3)"; + std::stringstream s1; + s1 << a1; + VERIFY_IS_EQUAL(s1.str(), std::string(" one two three")); + a3 = a1 + std::string(" (") + a2 + std::string(")"); + VERIFY((a3==a3ref).all()); + + a3 = a1; + a3 += std::string(" (") + a2 + std::string(")"); + VERIFY((a3==a3ref).all()); + + a1.swap(a3); + VERIFY((a1==a3ref).all()); + VERIFY((a3!=a3ref).all()); +} diff --git a/test/array_reverse.cpp b/test/array_reverse.cpp index a5c0d37f9..c9d9f90c3 100644 --- a/test/array_reverse.cpp +++ b/test/array_reverse.cpp @@ -117,13 +117,11 @@ template<typename MatrixType> void reverse(const MatrixType& m) m2.colwise().reverseInPlace(); VERIFY_IS_APPROX(m2,m1.colwise().reverse().eval()); - /* m1.colwise().reverse()(r, c) = x; VERIFY_IS_APPROX(x, m1(rows - 1 - r, c)); m1.rowwise().reverse()(r, c) = x; VERIFY_IS_APPROX(x, m1(r, cols - 1 - c)); - */ } void test_array_reverse() diff --git a/test/boostmultiprec.cpp b/test/boostmultiprec.cpp new file mode 100644 index 000000000..e06e9bdaf --- /dev/null +++ b/test/boostmultiprec.cpp @@ -0,0 +1,201 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include <sstream> + +#ifdef EIGEN_TEST_MAX_SIZE +#undef EIGEN_TEST_MAX_SIZE +#endif + +#define EIGEN_TEST_MAX_SIZE 50 + +#ifdef EIGEN_TEST_PART_1 +#include "cholesky.cpp" +#endif + +#ifdef EIGEN_TEST_PART_2 +#include "lu.cpp" +#endif + +#ifdef EIGEN_TEST_PART_3 +#include "qr.cpp" +#endif + +#ifdef EIGEN_TEST_PART_4 +#include "qr_colpivoting.cpp" +#endif + +#ifdef EIGEN_TEST_PART_5 +#include "qr_fullpivoting.cpp" +#endif + +#ifdef EIGEN_TEST_PART_6 +#include "eigensolver_selfadjoint.cpp" +#endif + +#ifdef EIGEN_TEST_PART_7 +#include "eigensolver_generic.cpp" +#endif + +#ifdef EIGEN_TEST_PART_8 +#include "eigensolver_generalized_real.cpp" +#endif + +#ifdef EIGEN_TEST_PART_9 +#include "jacobisvd.cpp" +#endif + +#ifdef EIGEN_TEST_PART_10 +#include "bdcsvd.cpp" +#endif + +#include <Eigen/Dense> + +#undef min +#undef max +#undef isnan +#undef isinf +#undef isfinite + +#include <boost/multiprecision/cpp_dec_float.hpp> +#include <boost/multiprecision/number.hpp> +#include <boost/math/special_functions.hpp> +#include <boost/math/complex.hpp> + +namespace mp = boost::multiprecision; +typedef mp::number<mp::cpp_dec_float<100>, mp::et_on> Real; + +namespace Eigen { + template<> struct NumTraits<Real> : GenericNumTraits<Real> { + static inline Real dummy_precision() { return 1e-50; } + }; + + template<typename T1,typename T2,typename T3,typename T4,typename T5> + struct NumTraits<boost::multiprecision::detail::expression<T1,T2,T3,T4,T5> > : NumTraits<Real> {}; + + template<> + Real test_precision<Real>() { return 1e-50; } + + // needed in C++93 mode where number does not support explicit cast. + namespace internal { + template<typename NewType> + struct cast_impl<Real,NewType> { + static inline NewType run(const Real& x) { + return x.template convert_to<NewType>(); + } + }; + + template<> + struct cast_impl<Real,std::complex<Real> > { + static inline std::complex<Real> run(const Real& x) { + return std::complex<Real>(x); + } + }; + } +} + +namespace boost { +namespace multiprecision { + // to make ADL works as expected: + using boost::math::isfinite; + using boost::math::isnan; + using boost::math::isinf; + using boost::math::copysign; + using boost::math::hypot; + + // The following is needed for std::complex<Real>: + Real fabs(const Real& a) { return abs EIGEN_NOT_A_MACRO (a); } + Real fmax(const Real& a, const Real& b) { using std::max; return max(a,b); } + + // some specialization for the unit tests: + inline bool test_isMuchSmallerThan(const Real& a, const Real& b) { + return internal::isMuchSmallerThan(a, b, test_precision<Real>()); + } + + inline bool test_isApprox(const Real& a, const Real& b) { + return internal::isApprox(a, b, test_precision<Real>()); + } + + inline bool test_isApproxOrLessThan(const Real& a, const Real& b) { + return internal::isApproxOrLessThan(a, b, test_precision<Real>()); + } + + Real get_test_precision(const Real&) { + return test_precision<Real>(); + } + + Real test_relative_error(const Real &a, const Real &b) { + using Eigen::numext::abs2; + return sqrt(abs2<Real>(a-b)/Eigen::numext::mini<Real>(abs2(a),abs2(b))); + } +} +} + +namespace Eigen { + +} + +void test_boostmultiprec() +{ + typedef Matrix<Real,Dynamic,Dynamic> Mat; + typedef Matrix<std::complex<Real>,Dynamic,Dynamic> MatC; + + std::cout << "NumTraits<Real>::epsilon() = " << NumTraits<Real>::epsilon() << std::endl; + std::cout << "NumTraits<Real>::dummy_precision() = " << NumTraits<Real>::dummy_precision() << std::endl; + std::cout << "NumTraits<Real>::lowest() = " << NumTraits<Real>::lowest() << std::endl; + std::cout << "NumTraits<Real>::highest() = " << NumTraits<Real>::highest() << std::endl; + std::cout << "NumTraits<Real>::digits10() = " << NumTraits<Real>::digits10() << std::endl; + + // chekc stream output + { + Mat A(10,10); + A.setRandom(); + std::stringstream ss; + ss << A; + } + { + MatC A(10,10); + A.setRandom(); + std::stringstream ss; + ss << A; + } + + for(int i = 0; i < g_repeat; i++) { + int s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); + + CALL_SUBTEST_1( cholesky(Mat(s,s)) ); + + CALL_SUBTEST_2( lu_non_invertible<Mat>() ); + CALL_SUBTEST_2( lu_invertible<Mat>() ); + CALL_SUBTEST_2( lu_non_invertible<MatC>() ); + CALL_SUBTEST_2( lu_invertible<MatC>() ); + + CALL_SUBTEST_3( qr(Mat(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) ); + CALL_SUBTEST_3( qr_invertible<Mat>() ); + + CALL_SUBTEST_4( qr<Mat>() ); + CALL_SUBTEST_4( cod<Mat>() ); + CALL_SUBTEST_4( qr_invertible<Mat>() ); + + CALL_SUBTEST_5( qr<Mat>() ); + CALL_SUBTEST_5( qr_invertible<Mat>() ); + + CALL_SUBTEST_6( selfadjointeigensolver(Mat(s,s)) ); + + CALL_SUBTEST_7( eigensolver(Mat(s,s)) ); + + CALL_SUBTEST_8( generalized_eigensolver_real(Mat(s,s)) ); + + TEST_SET_BUT_UNUSED_VARIABLE(s) + } + + 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))) )); + 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))) )); +} + diff --git a/test/bug1213.cpp b/test/bug1213.cpp new file mode 100644 index 000000000..581760c1a --- /dev/null +++ b/test/bug1213.cpp @@ -0,0 +1,13 @@ + +// This anonymous enum is essential to trigger the linking issue +enum { + Foo +}; + +#include "bug1213.h" + +bool bug1213_1(const Eigen::Vector3f& x) +{ + return bug1213_2(x); +} + diff --git a/test/bug1213.h b/test/bug1213.h new file mode 100644 index 000000000..040e5a470 --- /dev/null +++ b/test/bug1213.h @@ -0,0 +1,8 @@ + +#include <Eigen/Core> + +template<typename T, int dim> +bool bug1213_2(const Eigen::Matrix<T,dim,1>& x); + +bool bug1213_1(const Eigen::Vector3f& x); + diff --git a/test/bug1213_main.cpp b/test/bug1213_main.cpp new file mode 100644 index 000000000..4802c0003 --- /dev/null +++ b/test/bug1213_main.cpp @@ -0,0 +1,18 @@ + +// This is a regression unit regarding a weird linking issue with gcc. + +#include "bug1213.h" + +int main() +{ + return 0; +} + + +template<typename T, int dim> +bool bug1213_2(const Eigen::Matrix<T,dim,1>& ) +{ + return true; +} + +template bool bug1213_2<float,3>(const Eigen::Vector3f&); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index b7abc230b..8ad5ac639 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -154,6 +154,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m) SquareMatrixType symmLo = symm.template triangularView<Lower>(); LDLT<SquareMatrixType,Lower> ldltlo(symmLo); + VERIFY(ldltlo.info()==Success); VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); vecX = ldltlo.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); @@ -170,6 +171,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m) LDLT<SquareMatrixType,Upper> ldltup(symmUp); + VERIFY(ldltup.info()==Success); VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); vecX = ldltup.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); @@ -243,11 +245,13 @@ template<typename MatrixType> void cholesky(const MatrixType& m) // check matrices with a wide spectrum if(rows>=3) { + using std::pow; + using std::sqrt; RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8); Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,rows); Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(rows); for(Index k=0; k<rows; ++k) - d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); + d(k) = d(k)*pow(RealScalar(10),internal::random<RealScalar>(-s,s)); SquareMatrixType A = a * d.asDiagonal() * a.adjoint(); // Make sure a solution exists: vecX.setRandom(); @@ -263,7 +267,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m) } else { - RealScalar large_tol = std::sqrt(test_precision<RealScalar>()); + RealScalar large_tol = sqrt(test_precision<RealScalar>()); VERIFY((A * vecX).isApprox(vecB, large_tol)); ++g_test_level; @@ -329,6 +333,7 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m) RealMatrixType symmLo = symm.template triangularView<Lower>(); LDLT<RealMatrixType,Lower> ldltlo(symmLo); + VERIFY(ldltlo.info()==Success); VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); vecX = ldltlo.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); @@ -365,35 +370,90 @@ template<typename MatrixType> void cholesky_definiteness(const MatrixType& m) { mat << 1, 0, 0, -1; ldlt.compute(mat); + VERIFY(ldlt.info()==Success); VERIFY(!ldlt.isNegative()); VERIFY(!ldlt.isPositive()); } { mat << 1, 2, 2, 1; ldlt.compute(mat); + VERIFY(ldlt.info()==Success); VERIFY(!ldlt.isNegative()); VERIFY(!ldlt.isPositive()); } { mat << 0, 0, 0, 0; ldlt.compute(mat); + VERIFY(ldlt.info()==Success); VERIFY(ldlt.isNegative()); VERIFY(ldlt.isPositive()); } { mat << 0, 0, 0, 1; ldlt.compute(mat); + VERIFY(ldlt.info()==Success); VERIFY(!ldlt.isNegative()); VERIFY(ldlt.isPositive()); } { mat << -1, 0, 0, 0; ldlt.compute(mat); + VERIFY(ldlt.info()==Success); VERIFY(ldlt.isNegative()); VERIFY(!ldlt.isPositive()); } } +template<typename> +void cholesky_faillure_cases() +{ + MatrixXd mat; + LDLT<MatrixXd> ldlt; + + { + mat.resize(2,2); + mat << 0, 1, 1, 0; + ldlt.compute(mat); + VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); + VERIFY(ldlt.info()==NumericalIssue); + } +#if (!EIGEN_ARCH_i386) || defined(EIGEN_VECTORIZE_SSE2) + { + mat.resize(3,3); + mat << -1, -3, 3, + -3, -8.9999999999999999999, 1, + 3, 1, 0; + ldlt.compute(mat); + VERIFY(ldlt.info()==NumericalIssue); + VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); + } +#endif + { + mat.resize(3,3); + mat << 1, 2, 3, + 2, 4, 1, + 3, 1, 0; + ldlt.compute(mat); + VERIFY(ldlt.info()==NumericalIssue); + VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); + } + + { + mat.resize(8,8); + mat << 0.1, 0, -0.1, 0, 0, 0, 1, 0, + 0, 4.24667, 0, 2.00333, 0, 0, 0, 0, + -0.1, 0, 0.2, 0, -0.1, 0, 0, 0, + 0, 2.00333, 0, 8.49333, 0, 2.00333, 0, 0, + 0, 0, -0.1, 0, 0.1, 0, 0, 1, + 0, 0, 0, 2.00333, 0, 4.24667, 0, 0, + 1, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 0, 0; + ldlt.compute(mat); + VERIFY(ldlt.info()==NumericalIssue); + VERIFY_IS_NOT_APPROX(mat,ldlt.reconstructedMatrix()); + } +} + template<typename MatrixType> void cholesky_verify_assert() { MatrixType tmp; @@ -443,5 +503,7 @@ void test_cholesky() CALL_SUBTEST_9( LLT<MatrixXf>(10) ); CALL_SUBTEST_9( LDLT<MatrixXf>(10) ); + CALL_SUBTEST_2( cholesky_faillure_cases<void>() ); + TEST_SET_BUT_UNUSED_VARIABLE(nb_temporaries) } diff --git a/test/commainitializer.cpp b/test/commainitializer.cpp index 99102b966..9844adbd2 100644 --- a/test/commainitializer.cpp +++ b/test/commainitializer.cpp @@ -9,6 +9,62 @@ #include "main.h" + +template<int M1, int M2, int N1, int N2> +void test_blocks() +{ + Matrix<int, M1+M2, N1+N2> m_fixed; + MatrixXi m_dynamic(M1+M2, N1+N2); + + Matrix<int, M1, N1> mat11; mat11.setRandom(); + Matrix<int, M1, N2> mat12; mat12.setRandom(); + Matrix<int, M2, N1> mat21; mat21.setRandom(); + Matrix<int, M2, N2> mat22; mat22.setRandom(); + + MatrixXi matx11 = mat11, matx12 = mat12, matx21 = mat21, matx22 = mat22; + + { + VERIFY_IS_EQUAL((m_fixed << mat11, mat12, mat21, matx22).finished(), (m_dynamic << mat11, matx12, mat21, matx22).finished()); + VERIFY_IS_EQUAL((m_fixed.template topLeftCorner<M1,N1>()), mat11); + VERIFY_IS_EQUAL((m_fixed.template topRightCorner<M1,N2>()), mat12); + VERIFY_IS_EQUAL((m_fixed.template bottomLeftCorner<M2,N1>()), mat21); + VERIFY_IS_EQUAL((m_fixed.template bottomRightCorner<M2,N2>()), mat22); + VERIFY_IS_EQUAL((m_fixed << mat12, mat11, matx21, mat22).finished(), (m_dynamic << mat12, matx11, matx21, mat22).finished()); + } + + if(N1 > 0) + { + VERIFY_RAISES_ASSERT((m_fixed << mat11, mat12, mat11, mat21, mat22)); + VERIFY_RAISES_ASSERT((m_fixed << mat11, mat12, mat21, mat21, mat22)); + } + else + { + // allow insertion of zero-column blocks: + VERIFY_IS_EQUAL((m_fixed << mat11, mat12, mat11, mat11, mat21, mat21, mat22).finished(), (m_dynamic << mat12, mat22).finished()); + } + if(M1 != M2) + { + VERIFY_RAISES_ASSERT((m_fixed << mat11, mat21, mat12, mat22)); + } +} + + +template<int N> +struct test_block_recursion +{ + static void run() + { + test_blocks<(N>>6)&3, (N>>4)&3, (N>>2)&3, N & 3>(); + test_block_recursion<N-1>::run(); + } +}; + +template<> +struct test_block_recursion<-1> +{ + static void run() { } +}; + void test_commainitializer() { Matrix3d m3; @@ -43,4 +99,8 @@ void test_commainitializer() 4, 5, 6, vec[2].transpose(); VERIFY_IS_APPROX(m3, ref); + + + // recursively test all block-sizes from 0 to 3: + test_block_recursion<(1<<8) - 1>(); } diff --git a/test/cuda_basic.cu b/test/cuda_basic.cu index b36ed888d..cb2e4167a 100644 --- a/test/cuda_basic.cu +++ b/test/cuda_basic.cu @@ -1,4 +1,11 @@ - +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2015-2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. // workaround issue between gcc >= 4.7 and cuda 5.5 #if (defined __GNUC__) && (__GNUC__>4 || __GNUC_MINOR__>=7) @@ -12,10 +19,15 @@ #define EIGEN_DEFAULT_DENSE_INDEX_TYPE int #include <math_constants.h> +#include <cuda.h> +#if defined __CUDACC_VER__ && __CUDACC_VER__ >= 70500 +#include <cuda_fp16.h> +#endif #include "main.h" #include "cuda_common.h" -#include <Eigen/Eigenvalues> +// Check that dense modules can be properly parsed by nvcc +#include <Eigen/Dense> // struct Foo{ // EIGEN_DEVICE_FUNC diff --git a/test/dynalloc.cpp b/test/dynalloc.cpp index 5f587007c..f1cc70bee 100644 --- a/test/dynalloc.cpp +++ b/test/dynalloc.cpp @@ -22,7 +22,7 @@ void check_handmade_aligned_malloc() for(int i = 1; i < 1000; i++) { char *p = (char*)internal::handmade_aligned_malloc(i); - VERIFY(size_t(p)%ALIGNMENT==0); + VERIFY(internal::UIntPtr(p)%ALIGNMENT==0); // if the buffer is wrongly allocated this will give a bad write --> check with valgrind for(int j = 0; j < i; j++) p[j]=0; internal::handmade_aligned_free(p); @@ -34,7 +34,7 @@ void check_aligned_malloc() for(int i = ALIGNMENT; i < 1000; i++) { char *p = (char*)internal::aligned_malloc(i); - VERIFY(size_t(p)%ALIGNMENT==0); + VERIFY(internal::UIntPtr(p)%ALIGNMENT==0); // if the buffer is wrongly allocated this will give a bad write --> check with valgrind for(int j = 0; j < i; j++) p[j]=0; internal::aligned_free(p); @@ -46,7 +46,7 @@ void check_aligned_new() for(int i = ALIGNMENT; i < 1000; i++) { float *p = internal::aligned_new<float>(i); - VERIFY(size_t(p)%ALIGNMENT==0); + VERIFY(internal::UIntPtr(p)%ALIGNMENT==0); // if the buffer is wrongly allocated this will give a bad write --> check with valgrind for(int j = 0; j < i; j++) p[j]=0; internal::aligned_delete(p,i); @@ -58,7 +58,7 @@ void check_aligned_stack_alloc() for(int i = ALIGNMENT; i < 400; i++) { ei_declare_aligned_stack_constructed_variable(float,p,i,0); - VERIFY(size_t(p)%ALIGNMENT==0); + VERIFY(internal::UIntPtr(p)%ALIGNMENT==0); // if the buffer is wrongly allocated this will give a bad write --> check with valgrind for(int j = 0; j < i; j++) p[j]=0; } @@ -88,7 +88,7 @@ template<typename T> void check_dynaligned() { T* obj = new T; VERIFY(T::NeedsToAlign==1); - VERIFY(size_t(obj)%ALIGNMENT==0); + VERIFY(internal::UIntPtr(obj)%ALIGNMENT==0); delete obj; } } @@ -148,15 +148,15 @@ void test_dynalloc() } { - MyStruct foo0; VERIFY(size_t(foo0.avec.data())%ALIGNMENT==0); - MyClassA fooA; VERIFY(size_t(fooA.avec.data())%ALIGNMENT==0); + MyStruct foo0; VERIFY(internal::UIntPtr(foo0.avec.data())%ALIGNMENT==0); + MyClassA fooA; VERIFY(internal::UIntPtr(fooA.avec.data())%ALIGNMENT==0); } // dynamic allocation, single object for (int i=0; i<g_repeat*100; ++i) { - MyStruct *foo0 = new MyStruct(); VERIFY(size_t(foo0->avec.data())%ALIGNMENT==0); - MyClassA *fooA = new MyClassA(); VERIFY(size_t(fooA->avec.data())%ALIGNMENT==0); + MyStruct *foo0 = new MyStruct(); VERIFY(internal::UIntPtr(foo0->avec.data())%ALIGNMENT==0); + MyClassA *fooA = new MyClassA(); VERIFY(internal::UIntPtr(fooA->avec.data())%ALIGNMENT==0); delete foo0; delete fooA; } @@ -165,8 +165,8 @@ void test_dynalloc() const int N = 10; for (int i=0; i<g_repeat*100; ++i) { - MyStruct *foo0 = new MyStruct[N]; VERIFY(size_t(foo0->avec.data())%ALIGNMENT==0); - MyClassA *fooA = new MyClassA[N]; VERIFY(size_t(fooA->avec.data())%ALIGNMENT==0); + MyStruct *foo0 = new MyStruct[N]; VERIFY(internal::UIntPtr(foo0->avec.data())%ALIGNMENT==0); + MyClassA *fooA = new MyClassA[N]; VERIFY(internal::UIntPtr(fooA->avec.data())%ALIGNMENT==0); delete[] foo0; delete[] fooA; } diff --git a/test/eigensolver_generalized_real.cpp b/test/eigensolver_generalized_real.cpp index a46a2e50e..9c0838ba4 100644 --- a/test/eigensolver_generalized_real.cpp +++ b/test/eigensolver_generalized_real.cpp @@ -1,15 +1,17 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#define EIGEN_RUNTIME_NO_MALLOC #include "main.h" #include <limits> #include <Eigen/Eigenvalues> +#include <Eigen/LU> template<typename MatrixType> void generalized_eigensolver_real(const MatrixType& m) { @@ -21,6 +23,7 @@ template<typename MatrixType> void generalized_eigensolver_real(const MatrixType Index cols = m.cols(); typedef typename MatrixType::Scalar Scalar; + typedef std::complex<Scalar> ComplexScalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; MatrixType a = MatrixType::Random(rows,cols); @@ -31,14 +34,41 @@ template<typename MatrixType> void generalized_eigensolver_real(const MatrixType MatrixType spdB = b.adjoint() * b + b1.adjoint() * b1; // lets compare to GeneralizedSelfAdjointEigenSolver - GeneralizedSelfAdjointEigenSolver<MatrixType> symmEig(spdA, spdB); - GeneralizedEigenSolver<MatrixType> eig(spdA, spdB); + { + GeneralizedSelfAdjointEigenSolver<MatrixType> symmEig(spdA, spdB); + GeneralizedEigenSolver<MatrixType> eig(spdA, spdB); + + VERIFY_IS_EQUAL(eig.eigenvalues().imag().cwiseAbs().maxCoeff(), 0); + + VectorType realEigenvalues = eig.eigenvalues().real(); + std::sort(realEigenvalues.data(), realEigenvalues.data()+realEigenvalues.size()); + VERIFY_IS_APPROX(realEigenvalues, symmEig.eigenvalues()); - VERIFY_IS_EQUAL(eig.eigenvalues().imag().cwiseAbs().maxCoeff(), 0); + // check eigenvectors + typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType D = eig.eigenvalues().asDiagonal(); + typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType V = eig.eigenvectors(); + VERIFY_IS_APPROX(spdA*V, spdB*V*D); + } - VectorType realEigenvalues = eig.eigenvalues().real(); - std::sort(realEigenvalues.data(), realEigenvalues.data()+realEigenvalues.size()); - VERIFY_IS_APPROX(realEigenvalues, symmEig.eigenvalues()); + // non symmetric case: + { + GeneralizedEigenSolver<MatrixType> eig(rows); + // TODO enable full-prealocation of required memory, this probably requires an in-place mode for HessenbergDecomposition + //Eigen::internal::set_is_malloc_allowed(false); + eig.compute(a,b); + //Eigen::internal::set_is_malloc_allowed(true); + for(Index k=0; k<cols; ++k) + { + Matrix<ComplexScalar,Dynamic,Dynamic> tmp = (eig.betas()(k)*a).template cast<ComplexScalar>() - eig.alphas()(k)*b; + if(tmp.size()>1 && tmp.norm()>(std::numeric_limits<Scalar>::min)()) + tmp /= tmp.norm(); + VERIFY_IS_MUCH_SMALLER_THAN( std::abs(tmp.determinant()), Scalar(1) ); + } + // check eigenvectors + typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType D = eig.eigenvalues().asDiagonal(); + typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType V = eig.eigenvectors(); + VERIFY_IS_APPROX(a*V, b*V*D); + } // regression test for bug 1098 { @@ -57,7 +87,7 @@ void test_eigensolver_generalized_real() s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(s,s)) ); - // some trivial but implementation-wise tricky cases + // some trivial but implementation-wise special cases CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(1,1)) ); CALL_SUBTEST_2( generalized_eigensolver_real(MatrixXd(2,2)) ); CALL_SUBTEST_3( generalized_eigensolver_real(Matrix<double,1,1>()) ); diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index 566546310..e18fbf687 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -127,16 +127,29 @@ void test_eigensolver_generic() } ); - // regression test for bug 793 #ifdef EIGEN_TEST_PART_2 { - MatrixXd a(3,3); - a << 0, 0, 1, - 1, 1, 1, - 1, 1e+200, 1; - Eigen::EigenSolver<MatrixXd> eig(a); - VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()); - VERIFY_IS_APPROX(a * eig.eigenvectors(), eig.eigenvectors() * eig.eigenvalues().asDiagonal()); + // regression test for bug 793 + MatrixXd a(3,3); + a << 0, 0, 1, + 1, 1, 1, + 1, 1e+200, 1; + Eigen::EigenSolver<MatrixXd> eig(a); + double scale = 1e-200; // scale to avoid overflow during the comparisons + VERIFY_IS_APPROX(a * eig.pseudoEigenvectors()*scale, eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()*scale); + VERIFY_IS_APPROX(a * eig.eigenvectors()*scale, eig.eigenvectors() * eig.eigenvalues().asDiagonal()*scale); + } + { + // check a case where all eigenvalues are null. + MatrixXd a(2,2); + a << 1, 1, + -1, -1; + Eigen::EigenSolver<MatrixXd> eig(a); + VERIFY_IS_APPROX(eig.pseudoEigenvectors().squaredNorm(), 2.); + VERIFY_IS_APPROX((a * eig.pseudoEigenvectors()).norm()+1., 1.); + VERIFY_IS_APPROX((eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()).norm()+1., 1.); + VERIFY_IS_APPROX((a * eig.eigenvectors()).norm()+1., 1.); + VERIFY_IS_APPROX((eig.eigenvectors() * eig.eigenvalues().asDiagonal()).norm()+1., 1.); } #endif diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index f909761a1..4ed126116 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -12,18 +12,29 @@ #include "svd_fill.h" #include <limits> #include <Eigen/Eigenvalues> +#include <Eigen/SparseCore> template<typename MatrixType> void selfadjointeigensolver_essential_check(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - RealScalar eival_eps = (std::min)(test_precision<RealScalar>(), NumTraits<Scalar>::dummy_precision()*20000); + RealScalar eival_eps = numext::mini<RealScalar>(test_precision<RealScalar>(), NumTraits<Scalar>::dummy_precision()*20000); SelfAdjointEigenSolver<MatrixType> eiSymm(m); VERIFY_IS_EQUAL(eiSymm.info(), Success); - VERIFY_IS_APPROX(m.template selfadjointView<Lower>() * eiSymm.eigenvectors(), - eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal()); + + RealScalar scaling = m.cwiseAbs().maxCoeff(); + + if(scaling<(std::numeric_limits<RealScalar>::min)()) + { + VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)()); + } + else + { + VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiSymm.eigenvectors())/scaling, + (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal())/scaling); + } VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues()); VERIFY_IS_UNITARY(eiSymm.eigenvectors()); @@ -32,7 +43,6 @@ template<typename MatrixType> void selfadjointeigensolver_essential_check(const SelfAdjointEigenSolver<MatrixType> eiDirect; eiDirect.computeDirect(m); VERIFY_IS_EQUAL(eiDirect.info(), Success); - VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiDirect.eigenvalues()); if(! eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps) ) { std::cerr << "reference eigenvalues: " << eiSymm.eigenvalues().transpose() << "\n" @@ -40,10 +50,18 @@ template<typename MatrixType> void selfadjointeigensolver_essential_check(const << "diff: " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).transpose() << "\n" << "error (eps): " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).norm() / eiSymm.eigenvalues().norm() << " (" << eival_eps << ")\n"; } - VERIFY(eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps)); - VERIFY_IS_APPROX(m.template selfadjointView<Lower>() * eiDirect.eigenvectors(), - eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal()); - VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiDirect.eigenvalues()); + if(scaling<(std::numeric_limits<RealScalar>::min)()) + { + VERIFY(eiDirect.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)()); + } + else + { + VERIFY_IS_APPROX(eiSymm.eigenvalues()/scaling, eiDirect.eigenvalues()/scaling); + VERIFY_IS_APPROX((m.template selfadjointView<Lower>() * eiDirect.eigenvectors())/scaling, + (eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal())/scaling); + VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues()/scaling, eiDirect.eigenvalues()/scaling); + } + VERIFY_IS_UNITARY(eiDirect.eigenvectors()); } } @@ -164,6 +182,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) } } +template<int> void bug_854() { Matrix3d m; @@ -173,6 +192,7 @@ void bug_854() selfadjointeigensolver_essential_check(m); } +template<int> void bug_1014() { Matrix3d m; @@ -182,6 +202,26 @@ void bug_1014() selfadjointeigensolver_essential_check(m); } +template<int> +void bug_1225() +{ + Matrix3d m1, m2; + m1.setRandom(); + m1 = m1*m1.transpose(); + m2 = m1.triangularView<Upper>(); + SelfAdjointEigenSolver<Matrix3d> eig1(m1); + SelfAdjointEigenSolver<Matrix3d> eig2(m2.selfadjointView<Upper>()); + VERIFY_IS_APPROX(eig1.eigenvalues(), eig2.eigenvalues()); +} + +template<int> +void bug_1204() +{ + SparseMatrix<double> A(2,2); + A.setIdentity(); + SelfAdjointEigenSolver<Eigen::SparseMatrix<double> > eig(A); +} + void test_eigensolver_selfadjoint() { int s = 0; @@ -210,8 +250,10 @@ void test_eigensolver_selfadjoint() CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) ); } - CALL_SUBTEST_13( bug_854() ); - CALL_SUBTEST_13( bug_1014() ); + CALL_SUBTEST_13( bug_854<0>() ); + CALL_SUBTEST_13( bug_1014<0>() ); + CALL_SUBTEST_13( bug_1204<0>() ); + CALL_SUBTEST_13( bug_1225<0>() ); // Test problem size constructors s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4); diff --git a/test/evaluators.cpp b/test/evaluators.cpp index 876dffe22..aed5a05a7 100644 --- a/test/evaluators.cpp +++ b/test/evaluators.cpp @@ -21,7 +21,7 @@ namespace Eigen { EIGEN_STRONG_INLINE DstXprType& copy_using_evaluator(const EigenBase<DstXprType> &dst, const SrcXprType &src) { - call_assignment(dst.const_cast_derived(), src.derived(), internal::assign_op<typename DstXprType::Scalar>()); + call_assignment(dst.const_cast_derived(), src.derived(), internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>()); return dst.const_cast_derived(); } @@ -29,7 +29,7 @@ namespace Eigen { EIGEN_STRONG_INLINE const DstXprType& copy_using_evaluator(const NoAlias<DstXprType, StorageBase>& dst, const SrcXprType &src) { - call_assignment(dst, src.derived(), internal::assign_op<typename DstXprType::Scalar>()); + call_assignment(dst, src.derived(), internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>()); return dst.expression(); } @@ -45,7 +45,7 @@ namespace Eigen { dst.const_cast_derived().resizeLike(src.derived()); #endif - call_assignment(dst.const_cast_derived(), src.derived(), internal::assign_op<typename DstXprType::Scalar>()); + call_assignment(dst.const_cast_derived(), src.derived(), internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>()); return dst.const_cast_derived(); } @@ -53,28 +53,28 @@ namespace Eigen { void add_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src) { typedef typename DstXprType::Scalar Scalar; - call_assignment(const_cast<DstXprType&>(dst), src.derived(), internal::add_assign_op<Scalar>()); + call_assignment(const_cast<DstXprType&>(dst), src.derived(), internal::add_assign_op<Scalar,typename SrcXprType::Scalar>()); } template<typename DstXprType, typename SrcXprType> void subtract_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src) { typedef typename DstXprType::Scalar Scalar; - call_assignment(const_cast<DstXprType&>(dst), src.derived(), internal::sub_assign_op<Scalar>()); + call_assignment(const_cast<DstXprType&>(dst), src.derived(), internal::sub_assign_op<Scalar,typename SrcXprType::Scalar>()); } template<typename DstXprType, typename SrcXprType> void multiply_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src) { typedef typename DstXprType::Scalar Scalar; - call_assignment(dst.const_cast_derived(), src.derived(), internal::mul_assign_op<Scalar>()); + call_assignment(dst.const_cast_derived(), src.derived(), internal::mul_assign_op<Scalar,typename SrcXprType::Scalar>()); } template<typename DstXprType, typename SrcXprType> void divide_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src) { typedef typename DstXprType::Scalar Scalar; - call_assignment(dst.const_cast_derived(), src.derived(), internal::div_assign_op<Scalar>()); + call_assignment(dst.const_cast_derived(), src.derived(), internal::div_assign_op<Scalar,typename SrcXprType::Scalar>()); } template<typename DstXprType, typename SrcXprType> diff --git a/test/fastmath.cpp b/test/fastmath.cpp index efdd5b313..cc5db0746 100644 --- a/test/fastmath.cpp +++ b/test/fastmath.cpp @@ -49,7 +49,8 @@ void check_inf_nan(bool dryrun) { VERIFY( !m.allFinite() ); VERIFY( m.hasNaN() ); } - m(4) /= 0.0; + T hidden_zero = (std::numeric_limits<T>::min)()*(std::numeric_limits<T>::min)(); + m(4) /= hidden_zero; if(dryrun) { std::cout << "std::isfinite(" << m(4) << ") = "; check((std::isfinite)(m(4)),false); std::cout << " ; numext::isfinite = "; check((numext::isfinite)(m(4)), false); std::cout << "\n"; diff --git a/test/first_aligned.cpp b/test/first_aligned.cpp index bf22f6b97..ae2d4bc42 100644 --- a/test/first_aligned.cpp +++ b/test/first_aligned.cpp @@ -41,7 +41,7 @@ void test_first_aligned() test_first_aligned_helper(array_double+1, 50); test_first_aligned_helper(array_double+2, 50); - double *array_double_plus_4_bytes = (double*)(size_t(array_double)+4); + double *array_double_plus_4_bytes = (double*)(internal::UIntPtr(array_double)+4); test_none_aligned_helper(array_double_plus_4_bytes, 50); test_none_aligned_helper(array_double_plus_4_bytes+1, 50); diff --git a/test/geo_alignedbox.cpp b/test/geo_alignedbox.cpp index 2bdb4b7f2..d2339a651 100644 --- a/test/geo_alignedbox.cpp +++ b/test/geo_alignedbox.cpp @@ -48,6 +48,8 @@ template<typename BoxType> void alignedbox(const BoxType& _box) b0.extend(p0); b0.extend(p1); VERIFY(b0.contains(p0*s1+(Scalar(1)-s1)*p1)); + VERIFY(b0.contains(b0.center())); + VERIFY_IS_APPROX(b0.center(),(p0+p1)/Scalar(2)); (b2 = b0).extend(b1); VERIFY(b2.contains(b0)); diff --git a/test/geo_homogeneous.cpp b/test/geo_homogeneous.cpp index bf63c69ec..2187c7bf9 100644 --- a/test/geo_homogeneous.cpp +++ b/test/geo_homogeneous.cpp @@ -58,6 +58,8 @@ template<typename Scalar,int Size> void homogeneous(void) T2MatrixType t2 = T2MatrixType::Random(); VERIFY_IS_APPROX(t2 * (v0.homogeneous().eval()), t2 * v0.homogeneous()); VERIFY_IS_APPROX(t2 * (m0.colwise().homogeneous().eval()), t2 * m0.colwise().homogeneous()); + VERIFY_IS_APPROX(t2 * (v0.homogeneous().asDiagonal()), t2 * hv0.asDiagonal()); + VERIFY_IS_APPROX((v0.homogeneous().asDiagonal()) * t2, hv0.asDiagonal() * t2); VERIFY_IS_APPROX((v0.transpose().rowwise().homogeneous().eval()) * t2, v0.transpose().rowwise().homogeneous() * t2); @@ -109,6 +111,8 @@ template<typename Scalar,int Size> void homogeneous(void) VERIFY_IS_APPROX( (v0.transpose().homogeneous() .lazyProduct( t2 )).hnormalized(), (v0.transpose().homogeneous()*t2).hnormalized() ); VERIFY_IS_APPROX( (pts.transpose().rowwise().homogeneous() .lazyProduct( t2 )).rowwise().hnormalized(), (pts1.transpose()*t2).rowwise().hnormalized() ); + + VERIFY_IS_APPROX( (t2.template triangularView<Lower>() * v0.homogeneous()).eval(), (t2.template triangularView<Lower>()*hv0) ); } void test_geo_homogeneous() diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp index c1cc691c9..e77702bc7 100644 --- a/test/geo_hyperplane.cpp +++ b/test/geo_hyperplane.cpp @@ -97,9 +97,9 @@ template<typename Scalar> void lines() Vector u = Vector::Random(); Vector v = Vector::Random(); Scalar a = internal::random<Scalar>(); - while (abs(a-1) < 1e-4) a = internal::random<Scalar>(); - while (u.norm() < 1e-4) u = Vector::Random(); - while (v.norm() < 1e-4) v = Vector::Random(); + while (abs(a-1) < Scalar(1e-4)) a = internal::random<Scalar>(); + while (u.norm() < Scalar(1e-4)) u = Vector::Random(); + while (v.norm() < Scalar(1e-4)) v = Vector::Random(); HLine line_u = HLine::Through(center + u, center + a*u); HLine line_v = HLine::Through(center + v, center + a*v); @@ -111,14 +111,14 @@ template<typename Scalar> void lines() Vector result = line_u.intersection(line_v); // the lines should intersect at the point we called "center" - if(abs(a-1) > 1e-2 && abs(v.normalized().dot(u.normalized()))<0.9) + if(abs(a-1) > Scalar(1e-2) && abs(v.normalized().dot(u.normalized()))<Scalar(0.9)) VERIFY_IS_APPROX(result, center); // check conversions between two types of lines PLine pl(line_u); // gcc 3.3 will commit suicide if we don't name this variable HLine line_u2(pl); CoeffsType converted_coeffs = line_u2.coeffs(); - if(line_u2.normal().dot(line_u.normal())<0.) + if(line_u2.normal().dot(line_u.normal())<Scalar(0)) converted_coeffs = -line_u2.coeffs(); VERIFY(line_u.coeffs().isApprox(converted_coeffs)); } diff --git a/test/geo_quaternion.cpp b/test/geo_quaternion.cpp index 761bb52b4..96889e722 100644 --- a/test/geo_quaternion.cpp +++ b/test/geo_quaternion.cpp @@ -30,8 +30,8 @@ template<typename QuatType> void check_slerp(const QuatType& q0, const QuatType& Scalar largeEps = test_precision<Scalar>(); Scalar theta_tot = AA(q1*q0.inverse()).angle(); - if(theta_tot>EIGEN_PI) - theta_tot = Scalar(2.*EIGEN_PI)-theta_tot; + if(theta_tot>Scalar(EIGEN_PI)) + theta_tot = Scalar(2.)*Scalar(EIGEN_PI)-theta_tot; for(Scalar t=0; t<=Scalar(1.001); t+=Scalar(0.1)) { QuatType q = q0.slerp(t,q1); @@ -50,13 +50,12 @@ template<typename Scalar, int Options> void quaternion(void) using std::abs; typedef Matrix<Scalar,3,1> Vector3; typedef Matrix<Scalar,3,3> Matrix3; - typedef Matrix<Scalar,4,1> Vector4; typedef Quaternion<Scalar,Options> Quaternionx; typedef AngleAxis<Scalar> AngleAxisx; Scalar largeEps = test_precision<Scalar>(); if (internal::is_same<Scalar,float>::value) - largeEps = 1e-3f; + largeEps = Scalar(1e-3); Scalar eps = internal::random<Scalar>() * Scalar(1e-2); @@ -115,8 +114,8 @@ template<typename Scalar, int Options> void quaternion(void) // Do not execute the test if the rotation angle is almost zero, or // the rotation axis and v1 are almost parallel. if (abs(aa.angle()) > 5*test_precision<Scalar>() - && (aa.axis() - v1.normalized()).norm() < 1.99 - && (aa.axis() + v1.normalized()).norm() < 1.99) + && (aa.axis() - v1.normalized()).norm() < Scalar(1.99) + && (aa.axis() + v1.normalized()).norm() < Scalar(1.99)) { VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1); } @@ -157,8 +156,8 @@ template<typename Scalar, int Options> void quaternion(void) Quaternionx *q = new Quaternionx; delete q; - q1 = AngleAxisx(a, v0.normalized()); - q2 = AngleAxisx(b, v1.normalized()); + q1 = Quaternionx::UnitRandom(); + q2 = Quaternionx::UnitRandom(); check_slerp(q1,q2); q1 = AngleAxisx(b, v1.normalized()); @@ -169,7 +168,7 @@ template<typename Scalar, int Options> void quaternion(void) q2 = AngleAxisx(-b, -v1.normalized()); check_slerp(q1,q2); - q1.coeffs() = Vector4::Random().normalized(); + q1 = Quaternionx::UnitRandom(); q2.coeffs() = -q1.coeffs(); check_slerp(q1,q2); } diff --git a/test/geo_transformations.cpp b/test/geo_transformations.cpp index 51f90036d..278e527c2 100644 --- a/test/geo_transformations.cpp +++ b/test/geo_transformations.cpp @@ -18,6 +18,11 @@ Matrix<T,2,1> angleToVec(T a) return Matrix<T,2,1>(std::cos(a), std::sin(a)); } +// This permits to workaround a bug in clang/llvm code generation. +template<typename T> +EIGEN_DONT_INLINE +void dont_over_optimize(T& x) { volatile typename T::Scalar tmp = x(0); x(0) = tmp; } + template<typename Scalar, int Mode, int Options> void non_projective_only() { /* this test covers the following files: @@ -224,12 +229,13 @@ template<typename Scalar, int Mode, int Options> void transformations() do { v3 = Vector3::Random(); + dont_over_optimize(v3); } while (v3.cwiseAbs().minCoeff()<NumTraits<Scalar>::epsilon()); Translation3 tv3(v3); Transform3 t5(tv3); t4 = tv3; VERIFY_IS_APPROX(t5.matrix(), t4.matrix()); - t4.translate(-v3); + t4.translate((-v3).eval()); VERIFY_IS_APPROX(t4.matrix(), MatrixType::Identity()); t4 *= tv3; VERIFY_IS_APPROX(t5.matrix(), t4.matrix()); @@ -328,6 +334,9 @@ template<typename Scalar, int Mode, int Options> void transformations() t0.scale(v0); t1 *= AlignedScaling3(v0); VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); + t1 = AlignedScaling3(v0) * (Translation3(v0) * Transform3(q1)); + t1 = t1 * v0.asDiagonal(); + VERIFY_IS_APPROX(t0.matrix(), t1.matrix()); // transformation * translation t0.translate(v0); t1 = t1 * Translation3(v0); @@ -466,7 +475,7 @@ template<typename Scalar, int Mode, int Options> void transformations() Scalar a2 = R0.slerp(Scalar(k+1)/Scalar(path_steps), R1).angle(); l += std::abs(a2-a1); } - VERIFY(l<=EIGEN_PI*(Scalar(1)+NumTraits<Scalar>::epsilon()*Scalar(path_steps/2))); + VERIFY(l<=Scalar(EIGEN_PI)*(Scalar(1)+NumTraits<Scalar>::epsilon()*Scalar(path_steps/2))); // check basic features { @@ -476,6 +485,79 @@ template<typename Scalar, int Mode, int Options> void transformations() Rotation2D<Scalar> r2(r1); // copy ctor VERIFY_IS_APPROX(r2.angle(),s0); } + + { + Transform3 t32(Matrix4::Random()), t33, t34; + t34 = t33 = t32; + t32.scale(v0); + t33*=AlignedScaling3(v0); + VERIFY_IS_APPROX(t32.matrix(), t33.matrix()); + t33 = t34 * AlignedScaling3(v0); + VERIFY_IS_APPROX(t32.matrix(), t33.matrix()); + } + +} + +template<typename A1, typename A2, typename P, typename Q, typename V, typename H> +void transform_associativity_left(const A1& a1, const A2& a2, const P& p, const Q& q, const V& v, const H& h) +{ + VERIFY_IS_APPROX( q*(a1*v), (q*a1)*v ); + VERIFY_IS_APPROX( q*(a2*v), (q*a2)*v ); + VERIFY_IS_APPROX( q*(p*h).hnormalized(), ((q*p)*h).hnormalized() ); +} + +template<typename A1, typename A2, typename P, typename Q, typename V, typename H> +void transform_associativity2(const A1& a1, const A2& a2, const P& p, const Q& q, const V& v, const H& h) +{ + VERIFY_IS_APPROX( a1*(q*v), (a1*q)*v ); + VERIFY_IS_APPROX( a2*(q*v), (a2*q)*v ); + VERIFY_IS_APPROX( p *(q*v).homogeneous(), (p *q)*v.homogeneous() ); + + transform_associativity_left(a1, a2,p, q, v, h); +} + +template<typename Scalar, int Dim, int Options,typename RotationType> +void transform_associativity(const RotationType& R) +{ + typedef Matrix<Scalar,Dim,1> VectorType; + typedef Matrix<Scalar,Dim+1,1> HVectorType; + typedef Matrix<Scalar,Dim,Dim> LinearType; + typedef Matrix<Scalar,Dim+1,Dim+1> MatrixType; + typedef Transform<Scalar,Dim,AffineCompact,Options> AffineCompactType; + typedef Transform<Scalar,Dim,Affine,Options> AffineType; + typedef Transform<Scalar,Dim,Projective,Options> ProjectiveType; + typedef DiagonalMatrix<Scalar,Dim> ScalingType; + typedef Translation<Scalar,Dim> TranslationType; + + AffineCompactType A1c; A1c.matrix().setRandom(); + AffineCompactType A2c; A2c.matrix().setRandom(); + AffineType A1(A1c); + AffineType A2(A2c); + ProjectiveType P1; P1.matrix().setRandom(); + VectorType v1 = VectorType::Random(); + VectorType v2 = VectorType::Random(); + HVectorType h1 = HVectorType::Random(); + Scalar s1 = internal::random<Scalar>(); + LinearType L = LinearType::Random(); + MatrixType M = MatrixType::Random(); + + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, A2, v2, h1) ); + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, A2c, v2, h1) ); + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, v1.asDiagonal(), v2, h1) ); + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, ScalingType(v1), v2, h1) ); + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, Scaling(v1), v2, h1) ); + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, Scaling(s1), v2, h1) ); + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, TranslationType(v1), v2, h1) ); + CALL_SUBTEST( transform_associativity_left(A1c, A1, P1, L, v2, h1) ); + CALL_SUBTEST( transform_associativity2(A1c, A1, P1, R, v2, h1) ); + + VERIFY_IS_APPROX( A1*(M*h1), (A1*M)*h1 ); + VERIFY_IS_APPROX( A1c*(M*h1), (A1c*M)*h1 ); + VERIFY_IS_APPROX( P1*(M*h1), (P1*M)*h1 ); + + VERIFY_IS_APPROX( M*(A1*h1), (M*A1)*h1 ); + VERIFY_IS_APPROX( M*(A1c*h1), (M*A1c)*h1 ); + VERIFY_IS_APPROX( M*(P1*h1), ((M*P1)*h1) ); } template<typename Scalar> void transform_alignment() @@ -556,5 +638,8 @@ void test_geo_transformations() CALL_SUBTEST_7(( transform_products<double,3,RowMajor|AutoAlign>() )); CALL_SUBTEST_7(( transform_products<float,2,AutoAlign>() )); + + CALL_SUBTEST_8(( transform_associativity<double,2,ColMajor>(Rotation2D<double>(internal::random<double>()*double(EIGEN_PI))) )); + CALL_SUBTEST_8(( transform_associativity<double,3,ColMajor>(Quaterniond::UnitRandom()) )); } } diff --git a/test/half_float.cpp b/test/half_float.cpp new file mode 100644 index 000000000..f8d438e2f --- /dev/null +++ b/test/half_float.cpp @@ -0,0 +1,252 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include <sstream> + +#include "main.h" + +#include <Eigen/src/Core/arch/CUDA/Half.h> + +// Make sure it's possible to forward declare Eigen::half +namespace Eigen { +struct half; +} + +using Eigen::half; + +void test_conversion() +{ + using Eigen::half_impl::__half; + + // Conversion from float. + VERIFY_IS_EQUAL(half(1.0f).x, 0x3c00); + VERIFY_IS_EQUAL(half(0.5f).x, 0x3800); + VERIFY_IS_EQUAL(half(0.33333f).x, 0x3555); + VERIFY_IS_EQUAL(half(0.0f).x, 0x0000); + VERIFY_IS_EQUAL(half(-0.0f).x, 0x8000); + VERIFY_IS_EQUAL(half(65504.0f).x, 0x7bff); + VERIFY_IS_EQUAL(half(65536.0f).x, 0x7c00); // Becomes infinity. + + // Denormals. + VERIFY_IS_EQUAL(half(-5.96046e-08f).x, 0x8001); + VERIFY_IS_EQUAL(half(5.96046e-08f).x, 0x0001); + VERIFY_IS_EQUAL(half(1.19209e-07f).x, 0x0002); + + // Verify round-to-nearest-even behavior. + float val1 = float(half(__half(0x3c00))); + float val2 = float(half(__half(0x3c01))); + float val3 = float(half(__half(0x3c02))); + VERIFY_IS_EQUAL(half(0.5f * (val1 + val2)).x, 0x3c00); + VERIFY_IS_EQUAL(half(0.5f * (val2 + val3)).x, 0x3c02); + + // Conversion from int. + VERIFY_IS_EQUAL(half(-1).x, 0xbc00); + VERIFY_IS_EQUAL(half(0).x, 0x0000); + VERIFY_IS_EQUAL(half(1).x, 0x3c00); + VERIFY_IS_EQUAL(half(2).x, 0x4000); + VERIFY_IS_EQUAL(half(3).x, 0x4200); + + // Conversion from bool. + VERIFY_IS_EQUAL(half(false).x, 0x0000); + VERIFY_IS_EQUAL(half(true).x, 0x3c00); + + // Conversion to float. + VERIFY_IS_EQUAL(float(half(__half(0x0000))), 0.0f); + VERIFY_IS_EQUAL(float(half(__half(0x3c00))), 1.0f); + + // Denormals. + VERIFY_IS_APPROX(float(half(__half(0x8001))), -5.96046e-08f); + VERIFY_IS_APPROX(float(half(__half(0x0001))), 5.96046e-08f); + VERIFY_IS_APPROX(float(half(__half(0x0002))), 1.19209e-07f); + + // NaNs and infinities. + VERIFY(!(numext::isinf)(float(half(65504.0f)))); // Largest finite number. + VERIFY(!(numext::isnan)(float(half(0.0f)))); + VERIFY((numext::isinf)(float(half(__half(0xfc00))))); + VERIFY((numext::isnan)(float(half(__half(0xfc01))))); + VERIFY((numext::isinf)(float(half(__half(0x7c00))))); + VERIFY((numext::isnan)(float(half(__half(0x7c01))))); + +#if !EIGEN_COMP_MSVC + // Visual Studio errors out on divisions by 0 + VERIFY((numext::isnan)(float(half(0.0 / 0.0)))); + VERIFY((numext::isinf)(float(half(1.0 / 0.0)))); + VERIFY((numext::isinf)(float(half(-1.0 / 0.0)))); +#endif + + // Exactly same checks as above, just directly on the half representation. + VERIFY(!(numext::isinf)(half(__half(0x7bff)))); + VERIFY(!(numext::isnan)(half(__half(0x0000)))); + VERIFY((numext::isinf)(half(__half(0xfc00)))); + VERIFY((numext::isnan)(half(__half(0xfc01)))); + VERIFY((numext::isinf)(half(__half(0x7c00)))); + VERIFY((numext::isnan)(half(__half(0x7c01)))); + +#if !EIGEN_COMP_MSVC + // Visual Studio errors out on divisions by 0 + VERIFY((numext::isnan)(half(0.0 / 0.0))); + VERIFY((numext::isinf)(half(1.0 / 0.0))); + VERIFY((numext::isinf)(half(-1.0 / 0.0))); +#endif +} + +void test_numtraits() +{ + std::cout << "epsilon = " << NumTraits<half>::epsilon() << std::endl; + std::cout << "highest = " << NumTraits<half>::highest() << std::endl; + std::cout << "lowest = " << NumTraits<half>::lowest() << std::endl; + std::cout << "inifinty = " << NumTraits<half>::infinity() << std::endl; + std::cout << "nan = " << NumTraits<half>::quiet_NaN() << std::endl; + +} + +void test_arithmetic() +{ + VERIFY_IS_EQUAL(float(half(2) + half(2)), 4); + VERIFY_IS_EQUAL(float(half(2) + half(-2)), 0); + VERIFY_IS_APPROX(float(half(0.33333f) + half(0.66667f)), 1.0f); + VERIFY_IS_EQUAL(float(half(2.0f) * half(-5.5f)), -11.0f); + VERIFY_IS_APPROX(float(half(1.0f) / half(3.0f)), 0.33333f); + VERIFY_IS_EQUAL(float(-half(4096.0f)), -4096.0f); + VERIFY_IS_EQUAL(float(-half(-4096.0f)), 4096.0f); +} + +void test_comparison() +{ + VERIFY(half(1.0f) > half(0.5f)); + VERIFY(half(0.5f) < half(1.0f)); + VERIFY(!(half(1.0f) < half(0.5f))); + VERIFY(!(half(0.5f) > half(1.0f))); + + VERIFY(!(half(4.0f) > half(4.0f))); + VERIFY(!(half(4.0f) < half(4.0f))); + + VERIFY(!(half(0.0f) < half(-0.0f))); + VERIFY(!(half(-0.0f) < half(0.0f))); + VERIFY(!(half(0.0f) > half(-0.0f))); + VERIFY(!(half(-0.0f) > half(0.0f))); + + VERIFY(half(0.2f) > half(-1.0f)); + VERIFY(half(-1.0f) < half(0.2f)); + VERIFY(half(-16.0f) < half(-15.0f)); + + VERIFY(half(1.0f) == half(1.0f)); + VERIFY(half(1.0f) != half(2.0f)); + + // Comparisons with NaNs and infinities. +#if !EIGEN_COMP_MSVC + // Visual Studio errors out on divisions by 0 + VERIFY(!(half(0.0 / 0.0) == half(0.0 / 0.0))); + VERIFY(half(0.0 / 0.0) != half(0.0 / 0.0)); + + VERIFY(!(half(1.0) == half(0.0 / 0.0))); + VERIFY(!(half(1.0) < half(0.0 / 0.0))); + VERIFY(!(half(1.0) > half(0.0 / 0.0))); + VERIFY(half(1.0) != half(0.0 / 0.0)); + + VERIFY(half(1.0) < half(1.0 / 0.0)); + VERIFY(half(1.0) > half(-1.0 / 0.0)); +#endif +} + +void test_basic_functions() +{ + VERIFY_IS_EQUAL(float(numext::abs(half(3.5f))), 3.5f); + VERIFY_IS_EQUAL(float(abs(half(3.5f))), 3.5f); + VERIFY_IS_EQUAL(float(numext::abs(half(-3.5f))), 3.5f); + VERIFY_IS_EQUAL(float(abs(half(-3.5f))), 3.5f); + + VERIFY_IS_EQUAL(float(numext::floor(half(3.5f))), 3.0f); + VERIFY_IS_EQUAL(float(floor(half(3.5f))), 3.0f); + VERIFY_IS_EQUAL(float(numext::floor(half(-3.5f))), -4.0f); + VERIFY_IS_EQUAL(float(floor(half(-3.5f))), -4.0f); + + VERIFY_IS_EQUAL(float(numext::ceil(half(3.5f))), 4.0f); + VERIFY_IS_EQUAL(float(ceil(half(3.5f))), 4.0f); + VERIFY_IS_EQUAL(float(numext::ceil(half(-3.5f))), -3.0f); + VERIFY_IS_EQUAL(float(ceil(half(-3.5f))), -3.0f); + + VERIFY_IS_APPROX(float(numext::sqrt(half(0.0f))), 0.0f); + VERIFY_IS_APPROX(float(sqrt(half(0.0f))), 0.0f); + VERIFY_IS_APPROX(float(numext::sqrt(half(4.0f))), 2.0f); + VERIFY_IS_APPROX(float(sqrt(half(4.0f))), 2.0f); + + VERIFY_IS_APPROX(float(numext::pow(half(0.0f), half(1.0f))), 0.0f); + VERIFY_IS_APPROX(float(pow(half(0.0f), half(1.0f))), 0.0f); + VERIFY_IS_APPROX(float(numext::pow(half(2.0f), half(2.0f))), 4.0f); + VERIFY_IS_APPROX(float(pow(half(2.0f), half(2.0f))), 4.0f); + + VERIFY_IS_EQUAL(float(numext::exp(half(0.0f))), 1.0f); + VERIFY_IS_EQUAL(float(exp(half(0.0f))), 1.0f); + VERIFY_IS_APPROX(float(numext::exp(half(EIGEN_PI))), 20.f + float(EIGEN_PI)); + VERIFY_IS_APPROX(float(exp(half(EIGEN_PI))), 20.f + float(EIGEN_PI)); + + VERIFY_IS_EQUAL(float(numext::log(half(1.0f))), 0.0f); + VERIFY_IS_EQUAL(float(log(half(1.0f))), 0.0f); + VERIFY_IS_APPROX(float(numext::log(half(10.0f))), 2.30273f); + VERIFY_IS_APPROX(float(log(half(10.0f))), 2.30273f); + + VERIFY_IS_EQUAL(float(numext::log1p(half(0.0f))), 0.0f); + VERIFY_IS_EQUAL(float(log1p(half(0.0f))), 0.0f); + VERIFY_IS_APPROX(float(numext::log1p(half(10.0f))), 2.3978953f); + VERIFY_IS_APPROX(float(log1p(half(10.0f))), 2.3978953f); +} + +void test_trigonometric_functions() +{ + VERIFY_IS_APPROX(numext::cos(half(0.0f)), half(cosf(0.0f))); + VERIFY_IS_APPROX(cos(half(0.0f)), half(cosf(0.0f))); + VERIFY_IS_APPROX(numext::cos(half(EIGEN_PI)), half(cosf(EIGEN_PI))); + //VERIFY_IS_APPROX(numext::cos(half(EIGEN_PI/2)), half(cosf(EIGEN_PI/2))); + //VERIFY_IS_APPROX(numext::cos(half(3*EIGEN_PI/2)), half(cosf(3*EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::cos(half(3.5f)), half(cosf(3.5f))); + + VERIFY_IS_APPROX(numext::sin(half(0.0f)), half(sinf(0.0f))); + VERIFY_IS_APPROX(sin(half(0.0f)), half(sinf(0.0f))); + // VERIFY_IS_APPROX(numext::sin(half(EIGEN_PI)), half(sinf(EIGEN_PI))); + VERIFY_IS_APPROX(numext::sin(half(EIGEN_PI/2)), half(sinf(EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::sin(half(3*EIGEN_PI/2)), half(sinf(3*EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::sin(half(3.5f)), half(sinf(3.5f))); + + VERIFY_IS_APPROX(numext::tan(half(0.0f)), half(tanf(0.0f))); + VERIFY_IS_APPROX(tan(half(0.0f)), half(tanf(0.0f))); + // VERIFY_IS_APPROX(numext::tan(half(EIGEN_PI)), half(tanf(EIGEN_PI))); + // VERIFY_IS_APPROX(numext::tan(half(EIGEN_PI/2)), half(tanf(EIGEN_PI/2))); + //VERIFY_IS_APPROX(numext::tan(half(3*EIGEN_PI/2)), half(tanf(3*EIGEN_PI/2))); + VERIFY_IS_APPROX(numext::tan(half(3.5f)), half(tanf(3.5f))); +} + +void test_array() +{ + typedef Array<half,1,Dynamic> ArrayXh; + Index size = internal::random<Index>(1,10); + Index i = internal::random<Index>(0,size-1); + ArrayXh a1 = ArrayXh::Random(size), a2 = ArrayXh::Random(size); + VERIFY_IS_APPROX( a1+a1, half(2)*a1 ); + VERIFY( (a1.abs() >= half(0)).all() ); + VERIFY_IS_APPROX( (a1*a1).sqrt(), a1.abs() ); + + VERIFY( ((a1.min)(a2) <= (a1.max)(a2)).all() ); + a1(i) = half(-10.); + VERIFY_IS_EQUAL( a1.minCoeff(), half(-10.) ); + a1(i) = half(10.); + VERIFY_IS_EQUAL( a1.maxCoeff(), half(10.) ); + + std::stringstream ss; + ss << a1; +} + +void test_half_float() +{ + CALL_SUBTEST(test_conversion()); + CALL_SUBTEST(test_numtraits()); + CALL_SUBTEST(test_arithmetic()); + CALL_SUBTEST(test_comparison()); + CALL_SUBTEST(test_basic_functions()); + CALL_SUBTEST(test_trigonometric_functions()); + CALL_SUBTEST(test_array()); +} diff --git a/test/inplace_decomposition.cpp b/test/inplace_decomposition.cpp new file mode 100644 index 000000000..92d0d91b6 --- /dev/null +++ b/test/inplace_decomposition.cpp @@ -0,0 +1,110 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "main.h" +#include <Eigen/LU> +#include <Eigen/Cholesky> +#include <Eigen/QR> + +// This file test inplace decomposition through Ref<>, as supported by Cholesky, LU, and QR decompositions. + +template<typename DecType,typename MatrixType> void inplace(bool square = false, bool SPD = false) +{ + typedef typename MatrixType::Scalar Scalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> RhsType; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> ResType; + + Index rows = MatrixType::RowsAtCompileTime==Dynamic ? internal::random<Index>(2,EIGEN_TEST_MAX_SIZE/2) : Index(MatrixType::RowsAtCompileTime); + Index cols = MatrixType::ColsAtCompileTime==Dynamic ? (square?rows:internal::random<Index>(2,rows)) : Index(MatrixType::ColsAtCompileTime); + + MatrixType A = MatrixType::Random(rows,cols); + RhsType b = RhsType::Random(rows); + ResType x(cols); + + if(SPD) + { + assert(square); + A.topRows(cols) = A.topRows(cols).adjoint() * A.topRows(cols); + A.diagonal().array() += 1e-3; + } + + MatrixType A0 = A; + MatrixType A1 = A; + + DecType dec(A); + + // Check that the content of A has been modified + VERIFY_IS_NOT_APPROX( A, A0 ); + + // Check that the decomposition is correct: + if(rows==cols) + { + VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b ); + } + else + { + VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); + } + + // Check that modifying A breaks the current dec: + A.setRandom(); + if(rows==cols) + { + VERIFY_IS_NOT_APPROX( A0 * (x = dec.solve(b)), b ); + } + else + { + VERIFY_IS_NOT_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); + } + + // Check that calling compute(A1) does not modify A1: + A = A0; + dec.compute(A1); + VERIFY_IS_EQUAL(A0,A1); + VERIFY_IS_NOT_APPROX( A, A0 ); + if(rows==cols) + { + VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b ); + } + else + { + VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); + } +} + + +void test_inplace_decomposition() +{ + EIGEN_UNUSED typedef Matrix<double,4,3> Matrix43d; + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1(( inplace<LLT<Ref<MatrixXd> >, MatrixXd>(true,true) )); + CALL_SUBTEST_1(( inplace<LLT<Ref<Matrix4d> >, Matrix4d>(true,true) )); + + CALL_SUBTEST_2(( inplace<LDLT<Ref<MatrixXd> >, MatrixXd>(true,true) )); + CALL_SUBTEST_2(( inplace<LDLT<Ref<Matrix4d> >, Matrix4d>(true,true) )); + + CALL_SUBTEST_3(( inplace<PartialPivLU<Ref<MatrixXd> >, MatrixXd>(true,false) )); + CALL_SUBTEST_3(( inplace<PartialPivLU<Ref<Matrix4d> >, Matrix4d>(true,false) )); + + CALL_SUBTEST_4(( inplace<FullPivLU<Ref<MatrixXd> >, MatrixXd>(true,false) )); + CALL_SUBTEST_4(( inplace<FullPivLU<Ref<Matrix4d> >, Matrix4d>(true,false) )); + + CALL_SUBTEST_5(( inplace<HouseholderQR<Ref<MatrixXd> >, MatrixXd>(false,false) )); + CALL_SUBTEST_5(( inplace<HouseholderQR<Ref<Matrix43d> >, Matrix43d>(false,false) )); + + CALL_SUBTEST_6(( inplace<ColPivHouseholderQR<Ref<MatrixXd> >, MatrixXd>(false,false) )); + CALL_SUBTEST_6(( inplace<ColPivHouseholderQR<Ref<Matrix43d> >, Matrix43d>(false,false) )); + + CALL_SUBTEST_7(( inplace<FullPivHouseholderQR<Ref<MatrixXd> >, MatrixXd>(false,false) )); + CALL_SUBTEST_7(( inplace<FullPivHouseholderQR<Ref<Matrix43d> >, Matrix43d>(false,false) )); + + CALL_SUBTEST_8(( inplace<CompleteOrthogonalDecomposition<Ref<MatrixXd> >, MatrixXd>(false,false) )); + CALL_SUBTEST_8(( inplace<CompleteOrthogonalDecomposition<Ref<Matrix43d> >, Matrix43d>(false,false) )); + } +} diff --git a/test/integer_types.cpp b/test/integer_types.cpp index 950f8e9be..a21f73a81 100644 --- a/test/integer_types.cpp +++ b/test/integer_types.cpp @@ -158,4 +158,12 @@ void test_integer_types() CALL_SUBTEST_8( integer_type_tests(Matrix<unsigned long long, Dynamic, 5>(1, 5)) ); } +#ifdef EIGEN_TEST_PART_9 + VERIFY_IS_EQUAL(internal::scalar_div_cost<int>::value, 8); + VERIFY_IS_EQUAL(internal::scalar_div_cost<unsigned int>::value, 8); + if(sizeof(long)>sizeof(int)) { + VERIFY(internal::scalar_div_cost<long>::value > internal::scalar_div_cost<int>::value); + VERIFY(internal::scalar_div_cost<unsigned long>::value > internal::scalar_div_cost<int>::value); + } +#endif } diff --git a/test/is_same_dense.cpp b/test/is_same_dense.cpp index 6d7904bac..2c7838ce9 100644 --- a/test/is_same_dense.cpp +++ b/test/is_same_dense.cpp @@ -9,6 +9,8 @@ #include "main.h" +using internal::is_same_dense; + void test_is_same_dense() { typedef Matrix<double,Dynamic,Dynamic,ColMajor> ColMatrixXd; diff --git a/test/linearstructure.cpp b/test/linearstructure.cpp index 292f33969..17474af10 100644 --- a/test/linearstructure.cpp +++ b/test/linearstructure.cpp @@ -9,7 +9,7 @@ // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. static bool g_called; -#define EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN { g_called = true; } +#define EIGEN_SCALAR_BINARY_OP_PLUGIN { g_called |= (!internal::is_same<LhsScalar,RhsScalar>::value); } #include "main.h" @@ -21,6 +21,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m) */ typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; Index rows = m.rows(); Index cols = m.cols(); @@ -32,7 +33,7 @@ template<typename MatrixType> void linearStructure(const MatrixType& m) m3(rows, cols); Scalar s1 = internal::random<Scalar>(); - while (abs(s1)<1e-3) s1 = internal::random<Scalar>(); + while (abs(s1)<RealScalar(1e-3)) s1 = internal::random<Scalar>(); Index r = internal::random<Index>(0, rows-1), c = internal::random<Index>(0, cols-1); @@ -92,6 +93,22 @@ template<typename MatrixType> void real_complex(DenseIndex rows = MatrixType::Ro g_called = false; VERIFY_IS_APPROX(m1/s, m1/Scalar(s)); VERIFY(g_called && "matrix<complex> / real not properly optimized"); + + g_called = false; + VERIFY_IS_APPROX(s+m1.array(), Scalar(s)+m1.array()); + VERIFY(g_called && "real + matrix<complex> not properly optimized"); + + g_called = false; + VERIFY_IS_APPROX(m1.array()+s, m1.array()+Scalar(s)); + VERIFY(g_called && "matrix<complex> + real not properly optimized"); + + g_called = false; + VERIFY_IS_APPROX(s-m1.array(), Scalar(s)-m1.array()); + VERIFY(g_called && "real - matrix<complex> not properly optimized"); + + g_called = false; + VERIFY_IS_APPROX(m1.array()-s, m1.array()-Scalar(s)); + VERIFY(g_called && "matrix<complex> - real not properly optimized"); } void test_linearstructure() diff --git a/test/main.h b/test/main.h index b0e3b7818..74ff96a23 100644 --- a/test/main.h +++ b/test/main.h @@ -279,8 +279,8 @@ inline void verify_impl(bool condition, const char *testname, const char *file, #define VERIFY_LE(a, b) ::verify_impl(a <= b, g_test_stack.back().c_str(), __FILE__, __LINE__, EI_PP_MAKE_STRING(a <= b)) -#define VERIFY_IS_EQUAL(a, b) VERIFY(test_is_equal(a, b)) -#define VERIFY_IS_NOT_EQUAL(a, b) VERIFY(!test_is_equal(a, b)) +#define VERIFY_IS_EQUAL(a, b) VERIFY(test_is_equal(a, b, true)) +#define VERIFY_IS_NOT_EQUAL(a, b) VERIFY(test_is_equal(a, b, false)) #define VERIFY_IS_APPROX(a, b) VERIFY(verifyIsApprox(a, b)) #define VERIFY_IS_NOT_APPROX(a, b) VERIFY(!test_isApprox(a, b)) #define VERIFY_IS_MUCH_SMALLER_THAN(a, b) VERIFY(test_isMuchSmallerThan(a, b)) @@ -302,7 +302,7 @@ namespace Eigen { template<typename T> inline typename NumTraits<T>::Real test_precision() { return NumTraits<T>::dummy_precision(); } template<> inline float test_precision<float>() { return 1e-3f; } template<> inline double test_precision<double>() { return 1e-6; } -template<> inline long double test_precision<long double>() { return 1e-6; } +template<> inline long double test_precision<long double>() { return 1e-6l; } template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); } template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); } template<> inline long double test_precision<std::complex<long double> >() { return test_precision<long double>(); } @@ -452,20 +452,20 @@ T test_relative_error(const AngleAxis<T> &a, const AngleAxis<T> &b) } template<typename Type1, typename Type2> -inline bool test_isApprox(const Type1& a, const Type2& b) +inline bool test_isApprox(const Type1& a, const Type2& b, typename Type1::Scalar* = 0) // Enabled for Eigen's type only { return a.isApprox(b, test_precision<typename Type1::Scalar>()); } // get_test_precision is a small wrapper to test_precision allowing to return the scalar precision for either scalars or expressions template<typename T> -typename NumTraits<typename T::Scalar>::Real get_test_precision(const typename T::Scalar* = 0) +typename NumTraits<typename T::Scalar>::Real get_test_precision(const T&, const typename T::Scalar* = 0) { return test_precision<typename NumTraits<typename T::Scalar>::Real>(); } template<typename T> -typename NumTraits<T>::Real get_test_precision(typename internal::enable_if<internal::is_arithmetic<typename NumTraits<T>::Real>::value, T>::type* = 0) +typename NumTraits<T>::Real get_test_precision(const T&,typename internal::enable_if<internal::is_arithmetic<typename NumTraits<T>::Real>::value, T>::type* = 0) { return test_precision<typename NumTraits<T>::Real>(); } @@ -477,7 +477,7 @@ inline bool verifyIsApprox(const Type1& a, const Type2& b) bool ret = test_isApprox(a,b); if(!ret) { - std::cerr << "Difference too large wrt tolerance " << get_test_precision<Type1>() << ", relative error is: " << test_relative_error(a,b) << std::endl; + std::cerr << "Difference too large wrt tolerance " << get_test_precision(a) << ", relative error is: " << test_relative_error(a,b) << std::endl; } return ret; } @@ -517,17 +517,17 @@ inline bool test_isUnitary(const MatrixBase<Derived>& m) // Forward declaration to avoid ICC warning template<typename T, typename U> -bool test_is_equal(const T& actual, const U& expected); +bool test_is_equal(const T& actual, const U& expected, bool expect_equal=true); template<typename T, typename U> -bool test_is_equal(const T& actual, const U& expected) +bool test_is_equal(const T& actual, const U& expected, bool expect_equal) { - if (actual==expected) + if ((actual==expected) == expect_equal) return true; // false: std::cerr - << std::endl << " actual = " << actual - << std::endl << " expected = " << expected << std::endl << std::endl; + << "\n actual = " << actual + << "\n expected " << (expect_equal ? "= " : "!=") << expected << "\n\n"; return false; } @@ -736,3 +736,8 @@ int main(int argc, char *argv[]) // remark #1572: floating-point equality and inequality comparisons are unreliable #pragma warning disable 279 383 1418 1572 #endif + +#ifdef _MSC_VER + // 4503 - decorated name length exceeded, name was truncated + #pragma warning( disable : 4503) +#endif diff --git a/test/mapped_matrix.cpp b/test/mapped_matrix.cpp index 88653e887..6a84c5897 100644 --- a/test/mapped_matrix.cpp +++ b/test/mapped_matrix.cpp @@ -25,7 +25,7 @@ template<typename VectorType> void map_class_vector(const VectorType& m) Scalar* array1 = internal::aligned_new<Scalar>(size); Scalar* array2 = internal::aligned_new<Scalar>(size); Scalar* array3 = new Scalar[size+1]; - Scalar* array3unaligned = (std::size_t(array3)%EIGEN_MAX_ALIGN_BYTES) == 0 ? array3+1 : array3; + Scalar* array3unaligned = (internal::UIntPtr(array3)%EIGEN_MAX_ALIGN_BYTES) == 0 ? array3+1 : array3; Scalar array4[EIGEN_TESTMAP_MAX_SIZE]; Map<VectorType, AlignedMax>(array1, size) = VectorType::Random(size); @@ -65,7 +65,7 @@ template<typename MatrixType> void map_class_matrix(const MatrixType& m) // array3unaligned -> unaligned pointer to heap Scalar* array3 = new Scalar[size+1]; for(int i = 0; i < size+1; i++) array3[i] = Scalar(1); - Scalar* array3unaligned = size_t(array3)%EIGEN_MAX_ALIGN_BYTES == 0 ? array3+1 : array3; + Scalar* array3unaligned = internal::UIntPtr(array3)%EIGEN_MAX_ALIGN_BYTES == 0 ? array3+1 : array3; Scalar array4[256]; if(size<=256) for(int i = 0; i < size; i++) array4[i] = Scalar(1); @@ -129,7 +129,7 @@ template<typename VectorType> void map_static_methods(const VectorType& m) Scalar* array1 = internal::aligned_new<Scalar>(size); Scalar* array2 = internal::aligned_new<Scalar>(size); Scalar* array3 = new Scalar[size+1]; - Scalar* array3unaligned = size_t(array3)%EIGEN_MAX_ALIGN_BYTES == 0 ? array3+1 : array3; + Scalar* array3unaligned = internal::UIntPtr(array3)%EIGEN_MAX_ALIGN_BYTES == 0 ? array3+1 : array3; VectorType::MapAligned(array1, size) = VectorType::Random(size); VectorType::Map(array2, size) = VectorType::Map(array1, size); diff --git a/test/mapstride.cpp b/test/mapstride.cpp index ee2414248..4858f8fea 100644 --- a/test/mapstride.cpp +++ b/test/mapstride.cpp @@ -23,7 +23,7 @@ template<int Alignment,typename VectorType> void map_class_vector(const VectorTy Scalar* a_array = internal::aligned_new<Scalar>(arraysize+1); Scalar* array = a_array; if(Alignment!=Aligned) - array = (Scalar*)(ptrdiff_t(a_array) + (internal::packet_traits<Scalar>::AlignedOnScalar?sizeof(Scalar):sizeof(typename NumTraits<Scalar>::Real))); + array = (Scalar*)(internal::IntPtr(a_array) + (internal::packet_traits<Scalar>::AlignedOnScalar?sizeof(Scalar):sizeof(typename NumTraits<Scalar>::Real))); { Map<VectorType, Alignment, InnerStride<3> > map(array, size); @@ -63,14 +63,14 @@ template<int Alignment,typename MatrixType> void map_class_matrix(const MatrixTy Scalar* a_array1 = internal::aligned_new<Scalar>(arraysize+1); Scalar* array1 = a_array1; if(Alignment!=Aligned) - array1 = (Scalar*)(std::ptrdiff_t(a_array1) + (internal::packet_traits<Scalar>::AlignedOnScalar?sizeof(Scalar):sizeof(typename NumTraits<Scalar>::Real))); + array1 = (Scalar*)(internal::IntPtr(a_array1) + (internal::packet_traits<Scalar>::AlignedOnScalar?sizeof(Scalar):sizeof(typename NumTraits<Scalar>::Real))); Scalar a_array2[256]; Scalar* array2 = a_array2; if(Alignment!=Aligned) - array2 = (Scalar*)(std::ptrdiff_t(a_array2) + (internal::packet_traits<Scalar>::AlignedOnScalar?sizeof(Scalar):sizeof(typename NumTraits<Scalar>::Real))); + array2 = (Scalar*)(internal::IntPtr(a_array2) + (internal::packet_traits<Scalar>::AlignedOnScalar?sizeof(Scalar):sizeof(typename NumTraits<Scalar>::Real))); else - array2 = (Scalar*)(((std::size_t(a_array2)+EIGEN_MAX_ALIGN_BYTES-1)/EIGEN_MAX_ALIGN_BYTES)*EIGEN_MAX_ALIGN_BYTES); + array2 = (Scalar*)(((internal::UIntPtr(a_array2)+EIGEN_MAX_ALIGN_BYTES-1)/EIGEN_MAX_ALIGN_BYTES)*EIGEN_MAX_ALIGN_BYTES); Index maxsize2 = a_array2 - array2 + 256; // test no inner stride and some dynamic outer stride diff --git a/test/mixingtypes.cpp b/test/mixingtypes.cpp index 0b381ec6c..ad9c2c652 100644 --- a/test/mixingtypes.cpp +++ b/test/mixingtypes.cpp @@ -23,10 +23,18 @@ #endif +static bool g_called; +#define EIGEN_SCALAR_BINARY_OP_PLUGIN { g_called |= (!internal::is_same<LhsScalar,RhsScalar>::value); } + #include "main.h" using namespace std; +#define VERIFY_MIX_SCALAR(XPR,REF) \ + g_called = false; \ + VERIFY_IS_APPROX(XPR,REF); \ + VERIFY( g_called && #XPR" not properly optimized"); + template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) { typedef std::complex<float> CF; @@ -42,6 +50,7 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) Mat_f mf = Mat_f::Random(size,size); Mat_d md = mf.template cast<double>(); + //Mat_d rd = md; Mat_cf mcf = Mat_cf::Random(size,size); Mat_cd mcd = mcf.template cast<complex<double> >(); Mat_cd rcd = mcd; @@ -54,25 +63,59 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) complex<float> scf = internal::random<complex<float> >(); complex<double> scd = internal::random<complex<double> >(); - mf+mf; - VERIFY_RAISES_ASSERT(mf+md); -#ifndef EIGEN_HAS_STD_RESULT_OF - // this one does not even compile with C++11 - VERIFY_RAISES_ASSERT(mf+mcf); -#endif + + float epsf = std::sqrt(std::numeric_limits<float> ::min EIGEN_EMPTY ()); + double epsd = std::sqrt(std::numeric_limits<double>::min EIGEN_EMPTY ()); + + while(std::abs(sf )<epsf) sf = internal::random<float>(); + while(std::abs(sd )<epsd) sf = internal::random<double>(); + while(std::abs(scf)<epsf) scf = internal::random<CF>(); + while(std::abs(scd)<epsd) scd = internal::random<CD>(); + +// VERIFY_RAISES_ASSERT(mf+md); // does not even compile #ifdef EIGEN_DONT_VECTORIZE VERIFY_RAISES_ASSERT(vf=vd); VERIFY_RAISES_ASSERT(vf+=vd); - VERIFY_RAISES_ASSERT(mcd=md); #endif // check scalar products - VERIFY_IS_APPROX(vcf * sf , vcf * complex<float>(sf)); - VERIFY_IS_APPROX(sd * vcd, complex<double>(sd) * vcd); - VERIFY_IS_APPROX(vf * scf , vf.template cast<complex<float> >() * scf); - VERIFY_IS_APPROX(scd * vd, scd * vd.template cast<complex<double> >()); + VERIFY_MIX_SCALAR(vcf * sf , vcf * complex<float>(sf)); + VERIFY_MIX_SCALAR(sd * vcd , complex<double>(sd) * vcd); + VERIFY_MIX_SCALAR(vf * scf , vf.template cast<complex<float> >() * scf); + VERIFY_MIX_SCALAR(scd * vd , scd * vd.template cast<complex<double> >()); + + VERIFY_MIX_SCALAR(vcf * 2 , vcf * complex<float>(2)); + VERIFY_MIX_SCALAR(vcf * 2.1 , vcf * complex<float>(2.1)); + VERIFY_MIX_SCALAR(2 * vcf, vcf * complex<float>(2)); + VERIFY_MIX_SCALAR(2.1 * vcf , vcf * complex<float>(2.1)); + + // check scalar quotients + VERIFY_MIX_SCALAR(vcf / sf , vcf / complex<float>(sf)); + VERIFY_MIX_SCALAR(vf / scf , vf.template cast<complex<float> >() / scf); + VERIFY_MIX_SCALAR(vf.array() / scf, vf.template cast<complex<float> >().array() / scf); + VERIFY_MIX_SCALAR(scd / vd.array() , scd / vd.template cast<complex<double> >().array()); + + // check scalar increment + VERIFY_MIX_SCALAR(vcf.array() + sf , vcf.array() + complex<float>(sf)); + VERIFY_MIX_SCALAR(sd + vcd.array(), complex<double>(sd) + vcd.array()); + VERIFY_MIX_SCALAR(vf.array() + scf, vf.template cast<complex<float> >().array() + scf); + VERIFY_MIX_SCALAR(scd + vd.array() , scd + vd.template cast<complex<double> >().array()); + + // check scalar subtractions + VERIFY_MIX_SCALAR(vcf.array() - sf , vcf.array() - complex<float>(sf)); + VERIFY_MIX_SCALAR(sd - vcd.array(), complex<double>(sd) - vcd.array()); + VERIFY_MIX_SCALAR(vf.array() - scf, vf.template cast<complex<float> >().array() - scf); + VERIFY_MIX_SCALAR(scd - vd.array() , scd - vd.template cast<complex<double> >().array()); + + // check scalar powers + VERIFY_MIX_SCALAR( pow(vcf.array(), sf), Eigen::pow(vcf.array(), complex<float>(sf)) ); + VERIFY_MIX_SCALAR( vcf.array().pow(sf) , Eigen::pow(vcf.array(), complex<float>(sf)) ); + VERIFY_MIX_SCALAR( pow(sd, vcd.array()), Eigen::pow(complex<double>(sd), vcd.array()) ); + VERIFY_MIX_SCALAR( Eigen::pow(vf.array(), scf), Eigen::pow(vf.template cast<complex<float> >().array(), scf) ); + VERIFY_MIX_SCALAR( vf.array().pow(scf) , Eigen::pow(vf.template cast<complex<float> >().array(), scf) ); + VERIFY_MIX_SCALAR( Eigen::pow(scd, vd.array()), Eigen::pow(scd, vd.template cast<complex<double> >().array()) ); // check dot product vf.dot(vf); @@ -184,6 +227,63 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) Mat_cd((scd * mcd * md.template cast<CD>().eval()).template triangularView<Upper>())); VERIFY_IS_APPROX(Mat_cd(rcd.template triangularView<Upper>() = scd * md * mcd), Mat_cd((scd * md.template cast<CD>().eval() * mcd).template triangularView<Upper>())); + + + VERIFY_IS_APPROX( md.array() * mcd.array(), md.template cast<CD>().eval().array() * mcd.array() ); + VERIFY_IS_APPROX( mcd.array() * md.array(), mcd.array() * md.template cast<CD>().eval().array() ); + + VERIFY_IS_APPROX( md.array() + mcd.array(), md.template cast<CD>().eval().array() + mcd.array() ); + VERIFY_IS_APPROX( mcd.array() + md.array(), mcd.array() + md.template cast<CD>().eval().array() ); + + VERIFY_IS_APPROX( md.array() - mcd.array(), md.template cast<CD>().eval().array() - mcd.array() ); + VERIFY_IS_APPROX( mcd.array() - md.array(), mcd.array() - md.template cast<CD>().eval().array() ); + + if(mcd.array().abs().minCoeff()>epsd) + { + VERIFY_IS_APPROX( md.array() / mcd.array(), md.template cast<CD>().eval().array() / mcd.array() ); + } + if(md.array().abs().minCoeff()>epsd) + { + VERIFY_IS_APPROX( mcd.array() / md.array(), mcd.array() / md.template cast<CD>().eval().array() ); + } + + if(md.array().abs().minCoeff()>epsd || mcd.array().abs().minCoeff()>epsd) + { + VERIFY_IS_APPROX( md.array().pow(mcd.array()), md.template cast<CD>().eval().array().pow(mcd.array()) ); + VERIFY_IS_APPROX( mcd.array().pow(md.array()), mcd.array().pow(md.template cast<CD>().eval().array()) ); + + VERIFY_IS_APPROX( pow(md.array(),mcd.array()), md.template cast<CD>().eval().array().pow(mcd.array()) ); + VERIFY_IS_APPROX( pow(mcd.array(),md.array()), mcd.array().pow(md.template cast<CD>().eval().array()) ); + } + + rcd = mcd; + VERIFY_IS_APPROX( rcd = md, md.template cast<CD>().eval() ); + rcd = mcd; + VERIFY_IS_APPROX( rcd += md, mcd + md.template cast<CD>().eval() ); + rcd = mcd; + VERIFY_IS_APPROX( rcd -= md, mcd - md.template cast<CD>().eval() ); + rcd = mcd; + VERIFY_IS_APPROX( rcd.array() *= md.array(), mcd.array() * md.template cast<CD>().eval().array() ); + rcd = mcd; + if(md.array().abs().minCoeff()>epsd) + { + VERIFY_IS_APPROX( rcd.array() /= md.array(), mcd.array() / md.template cast<CD>().eval().array() ); + } + + rcd = mcd; + VERIFY_IS_APPROX( rcd.noalias() += md + mcd*md, mcd + (md.template cast<CD>().eval()) + mcd*(md.template cast<CD>().eval())); + + VERIFY_IS_APPROX( rcd.noalias() = md*md, ((md*md).eval().template cast<CD>()) ); + rcd = mcd; + VERIFY_IS_APPROX( rcd.noalias() += md*md, mcd + ((md*md).eval().template cast<CD>()) ); + rcd = mcd; + VERIFY_IS_APPROX( rcd.noalias() -= md*md, mcd - ((md*md).eval().template cast<CD>()) ); + + VERIFY_IS_APPROX( rcd.noalias() = mcd + md*md, mcd + ((md*md).eval().template cast<CD>()) ); + rcd = mcd; + VERIFY_IS_APPROX( rcd.noalias() += mcd + md*md, mcd + mcd + ((md*md).eval().template cast<CD>()) ); + rcd = mcd; + VERIFY_IS_APPROX( rcd.noalias() -= mcd + md*md, - ((md*md).eval().template cast<CD>()) ); } void test_mixingtypes() diff --git a/test/nesting_ops.cpp b/test/nesting_ops.cpp index 2f5025305..a419b0e44 100644 --- a/test/nesting_ops.cpp +++ b/test/nesting_ops.cpp @@ -75,8 +75,8 @@ template <typename MatrixType> void run_nesting_ops_2(const MatrixType& _m) } else { - VERIFY( verify_eval_type<1>(2*m1, 2*m1) ); - VERIFY( verify_eval_type<2>(2*m1, m1) ); + VERIFY( verify_eval_type<2>(2*m1, 2*m1) ); + VERIFY( verify_eval_type<3>(2*m1, m1) ); } VERIFY( verify_eval_type<2>(m1+m1, m1+m1) ); VERIFY( verify_eval_type<3>(m1+m1, m1) ); diff --git a/test/nullary.cpp b/test/nullary.cpp index cb87695ee..9063c6de8 100644 --- a/test/nullary.cpp +++ b/test/nullary.cpp @@ -104,13 +104,29 @@ void testVectorType(const VectorType& base) template<typename MatrixType> void testMatrixType(const MatrixType& m) { + using std::abs; const Index rows = m.rows(); const Index cols = m.cols(); + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + + Scalar s1; + do { + s1 = internal::random<Scalar>(); + } while(abs(s1)<RealScalar(1e-5) && (!NumTraits<Scalar>::IsInteger)); MatrixType A; A.setIdentity(rows, cols); VERIFY(equalsIdentity(A)); VERIFY(equalsIdentity(MatrixType::Identity(rows, cols))); + + + A = MatrixType::Constant(rows,cols,s1); + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); + VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1)(i,j), s1 ); + VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1).coeff(i,j), s1 ); + VERIFY_IS_APPROX( A(i,j), s1 ); } void test_nullary() @@ -137,4 +153,47 @@ void test_nullary() // Assignment of a RowVectorXd to a MatrixXd (regression test for bug #79). VERIFY( (MatrixXd(RowVectorXd::LinSpaced(3, 0, 1)) - RowVector3d(0, 0.5, 1)).norm() < std::numeric_limits<double>::epsilon() ); #endif + +#ifdef EIGEN_TEST_PART_10 + // check some internal logic + VERIFY(( internal::has_nullary_operator<internal::scalar_constant_op<double> >::value )); + VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<double> >::value )); + VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<double> >::value )); + VERIFY(( internal::functor_has_linear_access<internal::scalar_constant_op<double> >::ret )); + + VERIFY(( !internal::has_nullary_operator<internal::scalar_identity_op<double> >::value )); + VERIFY(( !internal::has_unary_operator<internal::scalar_identity_op<double> >::value )); + VERIFY(( internal::has_binary_operator<internal::scalar_identity_op<double> >::value )); + VERIFY(( !internal::functor_has_linear_access<internal::scalar_identity_op<double> >::ret )); + + VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<float,float,false> >::value )); + VERIFY(( internal::has_unary_operator<internal::linspaced_op<float,float,false> >::value )); + VERIFY(( !internal::has_binary_operator<internal::linspaced_op<float,float,false> >::value )); + VERIFY(( internal::functor_has_linear_access<internal::linspaced_op<float,float,false> >::ret )); + + // Regression unit test for a weird MSVC bug. + // Search "nullary_wrapper_workaround_msvc" in CoreEvaluators.h for the details. + // See also traits<Ref>::match. + { + MatrixXf A = MatrixXf::Random(3,3); + Ref<const MatrixXf> R = 2.0*A; + VERIFY_IS_APPROX(R, A+A); + + Ref<const MatrixXf> R1 = MatrixXf::Random(3,3)+A; + + VectorXi V = VectorXi::Random(3); + Ref<const VectorXi> R2 = VectorXi::LinSpaced(3,1,3)+V; + VERIFY_IS_APPROX(R2, V+Vector3i(1,2,3)); + + VERIFY(( internal::has_nullary_operator<internal::scalar_constant_op<float> >::value )); + VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<float> >::value )); + VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<float> >::value )); + VERIFY(( internal::functor_has_linear_access<internal::scalar_constant_op<float> >::ret )); + + VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<int,int,false> >::value )); + VERIFY(( internal::has_unary_operator<internal::linspaced_op<int,int,false> >::value )); + VERIFY(( !internal::has_binary_operator<internal::linspaced_op<int,int,false> >::value )); + VERIFY(( internal::functor_has_linear_access<internal::linspaced_op<int,int,false> >::ret )); + } +#endif } diff --git a/test/packetmath.cpp b/test/packetmath.cpp index c2346e1cd..20addf1ad 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -9,7 +9,11 @@ // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "main.h" +#include "unsupported/Eigen/SpecialFunctions" +#if defined __GNUC__ && __GNUC__>=6 + #pragma GCC diagnostic ignored "-Wignored-attributes" +#endif // using namespace Eigen; namespace Eigen { @@ -368,7 +372,15 @@ template<typename Scalar> void packetmath_real() VERIFY_IS_EQUAL(std::exp(-std::numeric_limits<Scalar>::denorm_min()), data2[1]); } -#ifdef EIGEN_HAS_C99_MATH + if (PacketTraits::HasTanh) { + // NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details. + data1[0] = std::numeric_limits<Scalar>::quiet_NaN(); + packet_helper<internal::packet_traits<Scalar>::HasTanh,Packet> h; + h.store(data2, internal::ptanh(h.load(data1))); + VERIFY((numext::isnan)(data2[0])); + } + +#if EIGEN_HAS_C99_MATH { data1[0] = std::numeric_limits<Scalar>::quiet_NaN(); packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h; @@ -395,11 +407,12 @@ template<typename Scalar> void packetmath_real() data2[i] = internal::random<Scalar>(0,1) * std::pow(Scalar(10), internal::random<Scalar>(-6,6)); } - if(internal::random<float>(0,1)<0.1) + if(internal::random<float>(0,1)<0.1f) data1[internal::random<int>(0, PacketSize)] = 0; CHECK_CWISE1_IF(PacketTraits::HasSqrt, std::sqrt, internal::psqrt); CHECK_CWISE1_IF(PacketTraits::HasLog, std::log, internal::plog); -#if defined(EIGEN_HAS_C99_MATH) && (__cplusplus > 199711L) +#if EIGEN_HAS_C99_MATH && (__cplusplus > 199711L) + CHECK_CWISE1_IF(PacketTraits::HasLog1p, std::log1p, internal::plog1p); CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma); CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf); CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc); @@ -432,7 +445,7 @@ template<typename Scalar> void packetmath_real() // VERIFY_IS_EQUAL(std::log(std::numeric_limits<Scalar>::denorm_min()), data2[0]); VERIFY((numext::isnan)(data2[1])); - data1[0] = -1.0f; + data1[0] = Scalar(-1.0f); h.store(data2, internal::plog(h.load(data1))); VERIFY((numext::isnan)(data2[0])); #if !EIGEN_FAST_MATH diff --git a/test/prec_inverse_4x4.cpp b/test/prec_inverse_4x4.cpp index c4ef2d4bd..eb6ad18c9 100644 --- a/test/prec_inverse_4x4.cpp +++ b/test/prec_inverse_4x4.cpp @@ -53,14 +53,29 @@ template<typename MatrixType> void inverse_general_4x4(int repeat) // FIXME that 1.25 used to be 1.2 until we tested gcc 4.1 on 30 June 2010 and got 1.21. VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.25)); VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 20.0)); + + { + int s = 5;//internal::random<int>(4,10); + int i = 0;//internal::random<int>(0,s-4); + int j = 0;//internal::random<int>(0,s-4); + Matrix<Scalar,5,5> mat(s,s); + mat.setRandom(); + MatrixType submat = mat.template block<4,4>(i,j); + MatrixType mat_inv = mat.template block<4,4>(i,j).inverse(); + VERIFY_IS_APPROX(mat_inv, submat.inverse()); + mat.template block<4,4>(i,j) = submat.inverse(); + VERIFY_IS_APPROX(mat_inv, (mat.template block<4,4>(i,j))); + } } void test_prec_inverse_4x4() { CALL_SUBTEST_1((inverse_permutation_4x4<Matrix4f>())); CALL_SUBTEST_1(( inverse_general_4x4<Matrix4f>(200000 * g_repeat) )); + CALL_SUBTEST_1(( inverse_general_4x4<Matrix<float,4,4,RowMajor> >(200000 * g_repeat) )); CALL_SUBTEST_2((inverse_permutation_4x4<Matrix<double,4,4,RowMajor> >())); + CALL_SUBTEST_2(( inverse_general_4x4<Matrix<double,4,4,ColMajor> >(200000 * g_repeat) )); CALL_SUBTEST_2(( inverse_general_4x4<Matrix<double,4,4,RowMajor> >(200000 * g_repeat) )); CALL_SUBTEST_3((inverse_permutation_4x4<Matrix4cf>())); diff --git a/test/product.h b/test/product.h index 27976a4ae..3b6511270 100644 --- a/test/product.h +++ b/test/product.h @@ -119,6 +119,14 @@ template<typename MatrixType> void product(const MatrixType& m) res.noalias() -= square + m1 * m2.transpose(); VERIFY_IS_APPROX(res, square + m1 * m2.transpose()); + // test d ?= a-b*c rules + res.noalias() = square - m1 * m2.transpose(); + VERIFY_IS_APPROX(res, square - m1 * m2.transpose()); + res.noalias() += square - m1 * m2.transpose(); + VERIFY_IS_APPROX(res, 2*(square - m1 * m2.transpose())); + res.noalias() -= square - m1 * m2.transpose(); + VERIFY_IS_APPROX(res, square - m1 * m2.transpose()); + tm1 = m1; VERIFY_IS_APPROX(tm1.transpose() * v1, m1.transpose() * v1); @@ -160,6 +168,29 @@ template<typename MatrixType> void product(const MatrixType& m) VERIFY_IS_APPROX(res2.block(0,0,1,cols).noalias() = m1.block(0,0,1,cols) * square2, (ref2.row(0) = m1.row(0) * square2)); } + // vector.block() (see bug 1283) + { + RowVectorType w1(rows); + VERIFY_IS_APPROX(square * v1.block(0,0,rows,1), square * v1); + VERIFY_IS_APPROX(w1.noalias() = square * v1.block(0,0,rows,1), square * v1); + VERIFY_IS_APPROX(w1.block(0,0,rows,1).noalias() = square * v1.block(0,0,rows,1), square * v1); + + Matrix<Scalar,1,MatrixType::ColsAtCompileTime> w2(cols); + VERIFY_IS_APPROX(vc2.block(0,0,cols,1).transpose() * square2, vc2.transpose() * square2); + VERIFY_IS_APPROX(w2.noalias() = vc2.block(0,0,cols,1).transpose() * square2, vc2.transpose() * square2); + VERIFY_IS_APPROX(w2.block(0,0,1,cols).noalias() = vc2.block(0,0,cols,1).transpose() * square2, vc2.transpose() * square2); + + vc2 = square2.block(0,0,1,cols).transpose(); + VERIFY_IS_APPROX(square2.block(0,0,1,cols) * square2, vc2.transpose() * square2); + VERIFY_IS_APPROX(w2.noalias() = square2.block(0,0,1,cols) * square2, vc2.transpose() * square2); + VERIFY_IS_APPROX(w2.block(0,0,1,cols).noalias() = square2.block(0,0,1,cols) * square2, vc2.transpose() * square2); + + vc2 = square2.block(0,0,cols,1); + VERIFY_IS_APPROX(square2.block(0,0,cols,1).transpose() * square2, vc2.transpose() * square2); + VERIFY_IS_APPROX(w2.noalias() = square2.block(0,0,cols,1).transpose() * square2, vc2.transpose() * square2); + VERIFY_IS_APPROX(w2.block(0,0,1,cols).noalias() = square2.block(0,0,cols,1).transpose() * square2, vc2.transpose() * square2); + } + // inner product { Scalar x = square2.row(c) * square2.col(c2); @@ -196,4 +227,5 @@ template<typename MatrixType> void product(const MatrixType& m) VERIFY_IS_APPROX(square * (s1*(square*square)), s1 * square * square * square); VERIFY_IS_APPROX(square * (square*square).conjugate(), square * square.conjugate() * square.conjugate()); } + } diff --git a/test/product_extra.cpp b/test/product_extra.cpp index d253fd7ed..e4990ac8c 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp @@ -256,6 +256,51 @@ Index compute_block_size() return ret; } + + +template<int> +void bug_1308() +{ + int n = 10; + MatrixXd r(n,n); + VectorXd v = VectorXd::Random(n); + r = v * RowVectorXd::Ones(n); + VERIFY_IS_APPROX(r, v.rowwise().replicate(n)); + r = VectorXd::Ones(n) * v.transpose(); + VERIFY_IS_APPROX(r, v.rowwise().replicate(n).transpose()); + + Matrix4d ones44 = Matrix4d::Ones(); + Matrix4d m44 = Matrix4d::Ones() * Matrix4d::Ones(); + VERIFY_IS_APPROX(m44,Matrix4d::Constant(4)); + VERIFY_IS_APPROX(m44.noalias()=ones44*Matrix4d::Ones(), Matrix4d::Constant(4)); + VERIFY_IS_APPROX(m44.noalias()=ones44.transpose()*Matrix4d::Ones(), Matrix4d::Constant(4)); + VERIFY_IS_APPROX(m44.noalias()=Matrix4d::Ones()*ones44, Matrix4d::Constant(4)); + VERIFY_IS_APPROX(m44.noalias()=Matrix4d::Ones()*ones44.transpose(), Matrix4d::Constant(4)); + + typedef Matrix<double,4,4,RowMajor> RMatrix4d; + RMatrix4d r44 = Matrix4d::Ones() * Matrix4d::Ones(); + VERIFY_IS_APPROX(r44,Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=ones44*Matrix4d::Ones(), Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=ones44.transpose()*Matrix4d::Ones(), Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=Matrix4d::Ones()*ones44, Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=Matrix4d::Ones()*ones44.transpose(), Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=ones44*RMatrix4d::Ones(), Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=ones44.transpose()*RMatrix4d::Ones(), Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=RMatrix4d::Ones()*ones44, Matrix4d::Constant(4)); + VERIFY_IS_APPROX(r44.noalias()=RMatrix4d::Ones()*ones44.transpose(), Matrix4d::Constant(4)); + +// RowVector4d r4; + m44.setOnes(); + r44.setZero(); + VERIFY_IS_APPROX(r44.noalias() += m44.row(0).transpose() * RowVector4d::Ones(), ones44); + r44.setZero(); + VERIFY_IS_APPROX(r44.noalias() += m44.col(0) * RowVector4d::Ones(), ones44); + r44.setZero(); + VERIFY_IS_APPROX(r44.noalias() += Vector4d::Ones() * m44.row(0), ones44); + r44.setZero(); + VERIFY_IS_APPROX(r44.noalias() += Vector4d::Ones() * m44.col(0).transpose(), ones44); +} + void test_product_extra() { for(int i = 0; i < g_repeat; i++) { @@ -268,8 +313,10 @@ void test_product_extra() } CALL_SUBTEST_5( bug_127<0>() ); CALL_SUBTEST_5( bug_817<0>() ); + CALL_SUBTEST_5( bug_1308<0>() ); CALL_SUBTEST_6( unaligned_objects<0>() ); CALL_SUBTEST_7( compute_block_size<float>() ); CALL_SUBTEST_7( compute_block_size<double>() ); CALL_SUBTEST_7( compute_block_size<std::complex<double> >() ); + } diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index 5a3f3a01a..2bb19a681 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -56,6 +56,9 @@ template<typename MatrixType> void product_notemporary(const MatrixType& m) VERIFY_EVALUATION_COUNT( m3.noalias() = m3 + m1 * m2.transpose(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() += m3 + m1 * m2.transpose(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() -= m3 + m1 * m2.transpose(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m3 - m1 * m2.transpose(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() += m3 - m1 * m2.transpose(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= m3 - m1 * m2.transpose(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * m2.adjoint(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * (m1*s3+m2*s2).adjoint(), 1); diff --git a/test/product_small.cpp b/test/product_small.cpp index c35db6f65..fdfdd9f6c 100644 --- a/test/product_small.cpp +++ b/test/product_small.cpp @@ -12,6 +12,7 @@ #include <Eigen/LU> // regression test for bug 447 +template<int> void product1x1() { Matrix<float,1,3> matAstatic; @@ -177,15 +178,66 @@ void test_lazy_l3() CALL_SUBTEST(( test_lazy_all_layout<T,4,-1,-1>(4,cols,depth) )); } +template<typename T,int N,int M,int K> +void test_linear_but_not_vectorizable() +{ + // Check tricky cases for which the result of the product is a vector and thus must exhibit the LinearBit flag, + // but is not vectorizable along the linear dimension. + Index n = N==Dynamic ? internal::random<Index>(1,32) : N; + Index m = M==Dynamic ? internal::random<Index>(1,32) : M; + Index k = K==Dynamic ? internal::random<Index>(1,32) : K; + + { + Matrix<T,N,M+1> A; A.setRandom(n,m+1); + Matrix<T,M*2,K> B; B.setRandom(m*2,k); + Matrix<T,1,K> C; + Matrix<T,1,K> R; + + C.noalias() = A.template topLeftCorner<1,M>() * (B.template topRows<M>()+B.template bottomRows<M>()); + R.noalias() = A.template topLeftCorner<1,M>() * (B.template topRows<M>()+B.template bottomRows<M>()).eval(); + VERIFY_IS_APPROX(C,R); + } + + { + Matrix<T,M+1,N,RowMajor> A; A.setRandom(m+1,n); + Matrix<T,K,M*2,RowMajor> B; B.setRandom(k,m*2); + Matrix<T,K,1> C; + Matrix<T,K,1> R; + + C.noalias() = (B.template leftCols<M>()+B.template rightCols<M>()) * A.template topLeftCorner<M,1>(); + R.noalias() = (B.template leftCols<M>()+B.template rightCols<M>()).eval() * A.template topLeftCorner<M,1>(); + VERIFY_IS_APPROX(C,R); + } +} + +template<int Rows> +void bug_1311() +{ + Matrix< double, Rows, 2 > A; A.setRandom(); + Vector2d b = Vector2d::Random() ; + Matrix<double,Rows,1> res; + res.noalias() = 1. * (A * b); + VERIFY_IS_APPROX(res, A*b); + res.noalias() = 1.*A * b; + VERIFY_IS_APPROX(res, A*b); + res.noalias() = (1.*A).lazyProduct(b); + VERIFY_IS_APPROX(res, A*b); + res.noalias() = (1.*A).lazyProduct(1.*b); + VERIFY_IS_APPROX(res, A*b); + res.noalias() = (A).lazyProduct(1.*b); + VERIFY_IS_APPROX(res, A*b); +} + void test_product_small() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1( product(Matrix<float, 3, 2>()) ); - CALL_SUBTEST_2( product(Matrix<int, 3, 5>()) ); + CALL_SUBTEST_2( product(Matrix<int, 3, 17>()) ); + CALL_SUBTEST_8( product(Matrix<double, 3, 17>()) ); CALL_SUBTEST_3( product(Matrix3d()) ); CALL_SUBTEST_4( product(Matrix4d()) ); CALL_SUBTEST_5( product(Matrix4f()) ); - CALL_SUBTEST_6( product1x1() ); + CALL_SUBTEST_6( product1x1<0>() ); CALL_SUBTEST_11( test_lazy_l1<float>() ); CALL_SUBTEST_12( test_lazy_l2<float>() ); @@ -202,6 +254,13 @@ void test_product_small() CALL_SUBTEST_41( test_lazy_l1<std::complex<double> >() ); CALL_SUBTEST_42( test_lazy_l2<std::complex<double> >() ); CALL_SUBTEST_43( test_lazy_l3<std::complex<double> >() ); + + CALL_SUBTEST_7(( test_linear_but_not_vectorizable<float,2,1,Dynamic>() )); + CALL_SUBTEST_7(( test_linear_but_not_vectorizable<float,3,1,Dynamic>() )); + CALL_SUBTEST_7(( test_linear_but_not_vectorizable<float,2,1,16>() )); + + CALL_SUBTEST_6( bug_1311<3>() ); + CALL_SUBTEST_6( bug_1311<5>() ); } #ifdef EIGEN_TEST_PART_6 diff --git a/test/qr.cpp b/test/qr.cpp index 98738777f..dfcc1e8f9 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -86,7 +86,7 @@ template<typename MatrixType> void qr_invertible() VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); // This test is tricky if the determinant becomes too small. // Since we generate random numbers with magnitude rrange [0,1], the average determinant is 0.5^size - VERIFY_IS_MUCH_SMALLER_THAN( abs(absdet-qr.absDeterminant()), (max)(RealScalar(pow(0.5,size)),(max)(abs(absdet),abs(qr.absDeterminant()))) ); + VERIFY_IS_MUCH_SMALLER_THAN( abs(absdet-qr.absDeterminant()), numext::maxi(RealScalar(pow(0.5,size)),numext::maxi<RealScalar>(abs(absdet),abs(qr.absDeterminant()))) ); } diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 46c54b74f..057bb014c 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -93,6 +93,7 @@ void cod_fixedsize() { template<typename MatrixType> void qr() { + using std::sqrt; typedef typename MatrixType::Index Index; Index rows = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE), cols2 = internal::random<Index>(2,EIGEN_TEST_MAX_SIZE); @@ -120,14 +121,14 @@ template<typename MatrixType> void qr() // Verify that the absolute value of the diagonal elements in R are // non-increasing until they reach the singularity threshold. RealScalar threshold = - std::sqrt(RealScalar(rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon(); + sqrt(RealScalar(rows)) * numext::abs(r(0, 0)) * NumTraits<Scalar>::epsilon(); for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) { - RealScalar x = (std::abs)(r(i, i)); - RealScalar y = (std::abs)(r(i + 1, i + 1)); + RealScalar x = numext::abs(r(i, i)); + RealScalar y = numext::abs(r(i + 1, i + 1)); if (x < threshold && y < threshold) continue; if (!test_isApproxOrLessThan(y, x)) { for (Index j = 0; j < (std::min)(rows, cols); ++j) { - std::cout << "i = " << j << ", |r_ii| = " << (std::abs)(r(j, j)) << std::endl; + std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; } std::cout << "Failure at i=" << i << ", rank=" << rank << ", threshold=" << threshold << std::endl; @@ -144,6 +145,8 @@ template<typename MatrixType> void qr() template<typename MatrixType, int Cols2> void qr_fixedsize() { + using std::sqrt; + using std::abs; enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -169,14 +172,14 @@ template<typename MatrixType, int Cols2> void qr_fixedsize() // Verify that the absolute value of the diagonal elements in R are // non-increasing until they reache the singularity threshold. RealScalar threshold = - std::sqrt(RealScalar(Rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon(); + sqrt(RealScalar(Rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon(); for (Index i = 0; i < (std::min)(int(Rows), int(Cols)) - 1; ++i) { - RealScalar x = (std::abs)(r(i, i)); - RealScalar y = (std::abs)(r(i + 1, i + 1)); + RealScalar x = numext::abs(r(i, i)); + RealScalar y = numext::abs(r(i + 1, i + 1)); if (x < threshold && y < threshold) continue; if (!test_isApproxOrLessThan(y, x)) { for (Index j = 0; j < (std::min)(int(Rows), int(Cols)); ++j) { - std::cout << "i = " << j << ", |r_ii| = " << (std::abs)(r(j, j)) << std::endl; + std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; } std::cout << "Failure at i=" << i << ", rank=" << rank << ", threshold=" << threshold << std::endl; @@ -194,6 +197,8 @@ template<typename MatrixType, int Cols2> void qr_fixedsize() // page 3 for more detail. template<typename MatrixType> void qr_kahan_matrix() { + using std::sqrt; + using std::abs; typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -204,23 +209,25 @@ template<typename MatrixType> void qr_kahan_matrix() m1.setZero(rows,cols); RealScalar s = std::pow(NumTraits<RealScalar>::epsilon(), 1.0 / rows); RealScalar c = std::sqrt(1 - s*s); + RealScalar pow_s_i(1.0); // pow(s,i) for (Index i = 0; i < rows; ++i) { - m1(i, i) = pow(s, i); - m1.row(i).tail(rows - i - 1) = -pow(s, i) * c * MatrixType::Ones(1, rows - i - 1); + m1(i, i) = pow_s_i; + m1.row(i).tail(rows - i - 1) = -pow_s_i * c * MatrixType::Ones(1, rows - i - 1); + pow_s_i *= s; } m1 = (m1 + m1.transpose()).eval(); ColPivHouseholderQR<MatrixType> qr(m1); MatrixType r = qr.matrixQR().template triangularView<Upper>(); RealScalar threshold = - std::sqrt(RealScalar(rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon(); + std::sqrt(RealScalar(rows)) * numext::abs(r(0, 0)) * NumTraits<Scalar>::epsilon(); for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) { - RealScalar x = (std::abs)(r(i, i)); - RealScalar y = (std::abs)(r(i + 1, i + 1)); + RealScalar x = numext::abs(r(i, i)); + RealScalar y = numext::abs(r(i + 1, i + 1)); if (x < threshold && y < threshold) continue; if (!test_isApproxOrLessThan(y, x)) { for (Index j = 0; j < (std::min)(rows, cols); ++j) { - std::cout << "i = " << j << ", |r_ii| = " << (std::abs)(r(j, j)) << std::endl; + std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; } std::cout << "Failure at i=" << i << ", rank=" << qr.rank() << ", threshold=" << threshold << std::endl; diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index d82e123d0..05a705887 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -15,8 +15,12 @@ template<typename MatrixType> void qr() { typedef typename MatrixType::Index Index; - Index rows = internal::random<Index>(20,200), cols = internal::random<int>(20,200), cols2 = internal::random<int>(20,200); - Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1); + Index max_size = EIGEN_TEST_MAX_SIZE; + Index min_size = numext::maxi(1,EIGEN_TEST_MAX_SIZE/10); + Index rows = internal::random<Index>(min_size,max_size), + cols = internal::random<Index>(min_size,max_size), + cols2 = internal::random<Index>(min_size,max_size), + rank = internal::random<Index>(1, (std::min)(rows, cols)-1); typedef typename MatrixType::Scalar Scalar; typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; @@ -59,7 +63,9 @@ template<typename MatrixType> void qr_invertible() typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef typename MatrixType::Scalar Scalar; - int size = internal::random<int>(10,50); + Index max_size = numext::mini(50,EIGEN_TEST_MAX_SIZE); + Index min_size = numext::maxi(1,EIGEN_TEST_MAX_SIZE/10); + Index size = internal::random<Index>(min_size,max_size); MatrixType m1(size, size), m2(size, size), m3(size, size); m1 = MatrixType::Random(size,size); diff --git a/test/rand.cpp b/test/rand.cpp index eeec34191..51cf01773 100644 --- a/test/rand.cpp +++ b/test/rand.cpp @@ -9,6 +9,8 @@ #include "main.h" +typedef long long int64; + template<typename Scalar> Scalar check_in_range(Scalar x, Scalar y) { Scalar r = internal::random<Scalar>(x,y); @@ -35,31 +37,49 @@ template<typename Scalar> void check_all_in_range(Scalar x, Scalar y) VERIFY( (mask>0).all() ); } +template<typename Scalar> void check_histogram(Scalar x, Scalar y, int bins) +{ + Array<int,1,Dynamic> hist(bins); + hist.fill(0); + int f = 100000; + int n = bins*f; + int64 range = int64(y)-int64(x); + int divisor = int((range+1)/bins); + assert(((range+1)%bins)==0); + for(int k=0; k<n; ++k) + { + Scalar r = check_in_range(x,y); + hist( int((int64(r)-int64(x))/divisor) )++; + } + VERIFY( (((hist.cast<double>()/double(f))-1.0).abs()<0.02).all() ); +} + void test_rand() { long long_ref = NumTraits<long>::highest()/10; signed char char_offset = (std::min)(g_repeat,64); signed char short_offset = (std::min)(g_repeat,16000); - - for(int i = 0; i < g_repeat*10; i++) { + + for(int i = 0; i < g_repeat*10000; i++) { CALL_SUBTEST(check_in_range<float>(10,11)); CALL_SUBTEST(check_in_range<float>(1.24234523,1.24234523)); CALL_SUBTEST(check_in_range<float>(-1,1)); CALL_SUBTEST(check_in_range<float>(-1432.2352,-1432.2352)); - + CALL_SUBTEST(check_in_range<double>(10,11)); CALL_SUBTEST(check_in_range<double>(1.24234523,1.24234523)); CALL_SUBTEST(check_in_range<double>(-1,1)); CALL_SUBTEST(check_in_range<double>(-1432.2352,-1432.2352)); - + CALL_SUBTEST(check_in_range<int>(0,-1)); CALL_SUBTEST(check_in_range<short>(0,-1)); CALL_SUBTEST(check_in_range<long>(0,-1)); CALL_SUBTEST(check_in_range<int>(-673456,673456)); + CALL_SUBTEST(check_in_range<int>(-RAND_MAX+10,RAND_MAX-10)); CALL_SUBTEST(check_in_range<short>(-24345,24345)); CALL_SUBTEST(check_in_range<long>(-long_ref,long_ref)); } - + CALL_SUBTEST(check_all_in_range<signed char>(11,11)); CALL_SUBTEST(check_all_in_range<signed char>(11,11+char_offset)); CALL_SUBTEST(check_all_in_range<signed char>(-5,5)); @@ -67,25 +87,32 @@ void test_rand() CALL_SUBTEST(check_all_in_range<signed char>(-126,-126+char_offset)); CALL_SUBTEST(check_all_in_range<signed char>(126-char_offset,126)); CALL_SUBTEST(check_all_in_range<signed char>(-126,126)); - + CALL_SUBTEST(check_all_in_range<short>(11,11)); CALL_SUBTEST(check_all_in_range<short>(11,11+short_offset)); CALL_SUBTEST(check_all_in_range<short>(-5,5)); CALL_SUBTEST(check_all_in_range<short>(-11-short_offset,-11)); CALL_SUBTEST(check_all_in_range<short>(-24345,-24345+short_offset)); CALL_SUBTEST(check_all_in_range<short>(24345,24345+short_offset)); - + CALL_SUBTEST(check_all_in_range<int>(11,11)); CALL_SUBTEST(check_all_in_range<int>(11,11+g_repeat)); CALL_SUBTEST(check_all_in_range<int>(-5,5)); CALL_SUBTEST(check_all_in_range<int>(-11-g_repeat,-11)); CALL_SUBTEST(check_all_in_range<int>(-673456,-673456+g_repeat)); CALL_SUBTEST(check_all_in_range<int>(673456,673456+g_repeat)); - + CALL_SUBTEST(check_all_in_range<long>(11,11)); CALL_SUBTEST(check_all_in_range<long>(11,11+g_repeat)); CALL_SUBTEST(check_all_in_range<long>(-5,5)); CALL_SUBTEST(check_all_in_range<long>(-11-g_repeat,-11)); CALL_SUBTEST(check_all_in_range<long>(-long_ref,-long_ref+g_repeat)); CALL_SUBTEST(check_all_in_range<long>( long_ref, long_ref+g_repeat)); + + CALL_SUBTEST(check_histogram<int>(-5,5,11)); + int bins = 100; + CALL_SUBTEST(check_histogram<int>(-3333,-3333+bins*(3333/bins)-1,bins)); + bins = 1000; + CALL_SUBTEST(check_histogram<int>(-RAND_MAX+10,-RAND_MAX+10+bins*(RAND_MAX/bins)-1,bins)); + CALL_SUBTEST(check_histogram<int>(-RAND_MAX+10,-int64(RAND_MAX)+10+bins*(2*int64(RAND_MAX)/bins)-1,bins)); } diff --git a/test/real_qz.cpp b/test/real_qz.cpp index a1766c6d9..99ac31235 100644 --- a/test/real_qz.cpp +++ b/test/real_qz.cpp @@ -7,6 +7,7 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#define EIGEN_RUNTIME_NO_MALLOC #include "main.h" #include <limits> #include <Eigen/Eigenvalues> @@ -41,7 +42,11 @@ template<typename MatrixType> void real_qz(const MatrixType& m) break; } - RealQZ<MatrixType> qz(A,B); + RealQZ<MatrixType> qz(dim); + // TODO enable full-prealocation of required memory, this probably requires an in-place mode for HessenbergDecomposition + //Eigen::internal::set_is_malloc_allowed(false); + qz.compute(A,B); + //Eigen::internal::set_is_malloc_allowed(true); VERIFY_IS_EQUAL(qz.info(), Success); // check for zeros @@ -49,11 +54,20 @@ template<typename MatrixType> void real_qz(const MatrixType& m) for (Index i=0; i<A.cols(); i++) for (Index j=0; j<i; j++) { if (abs(qz.matrixT()(i,j))!=Scalar(0.0)) + { + std::cerr << "Error: T(" << i << "," << j << ") = " << qz.matrixT()(i,j) << std::endl; all_zeros = false; + } if (j<i-1 && abs(qz.matrixS()(i,j))!=Scalar(0.0)) + { + std::cerr << "Error: S(" << i << "," << j << ") = " << qz.matrixS()(i,j) << std::endl; all_zeros = false; + } if (j==i-1 && j>0 && abs(qz.matrixS()(i,j))!=Scalar(0.0) && abs(qz.matrixS()(i-1,j-1))!=Scalar(0.0)) + { + std::cerr << "Error: S(" << i << "," << j << ") = " << qz.matrixS()(i,j) << " && S(" << i-1 << "," << j-1 << ") = " << qz.matrixS()(i-1,j-1) << std::endl; all_zeros = false; + } } VERIFY_IS_EQUAL(all_zeros, true); VERIFY_IS_APPROX(qz.matrixQ()*qz.matrixS()*qz.matrixZ(), A); diff --git a/test/rvalue_types.cpp b/test/rvalue_types.cpp index 3eebfc61b..8887f1b1b 100644 --- a/test/rvalue_types.cpp +++ b/test/rvalue_types.cpp @@ -11,7 +11,9 @@ #include <Eigen/Core> -#ifdef EIGEN_HAVE_RVALUE_REFERENCES +using internal::UIntPtr; + +#if EIGEN_HAS_RVALUE_REFERENCES template <typename MatrixType> void rvalue_copyassign(const MatrixType& m) { @@ -20,11 +22,11 @@ void rvalue_copyassign(const MatrixType& m) // create a temporary which we are about to destroy by moving MatrixType tmp = m; - long src_address = reinterpret_cast<long>(tmp.data()); + UIntPtr src_address = reinterpret_cast<UIntPtr>(tmp.data()); // move the temporary to n MatrixType n = std::move(tmp); - long dst_address = reinterpret_cast<long>(n.data()); + UIntPtr dst_address = reinterpret_cast<UIntPtr>(n.data()); if (MatrixType::RowsAtCompileTime==Dynamic|| MatrixType::ColsAtCompileTime==Dynamic) { diff --git a/test/schur_real.cpp b/test/schur_real.cpp index cfe4570d4..4aede87df 100644 --- a/test/schur_real.cpp +++ b/test/schur_real.cpp @@ -82,7 +82,7 @@ template<typename MatrixType> void schur(int size = MatrixType::ColsAtCompileTim Atriangular.template triangularView<StrictlyLower>().setZero(); rs3.setMaxIterations(1).compute(Atriangular); // triangular matrices do not need any iterations VERIFY_IS_EQUAL(rs3.info(), Success); - VERIFY_IS_EQUAL(rs3.matrixT(), Atriangular); + VERIFY_IS_APPROX(rs3.matrixT(), Atriangular); // approx because of scaling... VERIFY_IS_EQUAL(rs3.matrixU(), MatrixType::Identity(size, size)); // Test computation of only T, not U diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index cb8ebaedf..7b5f3eb38 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -157,18 +157,15 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re initSparse<Scalar>(density, refM3, m3); initSparse<Scalar>(density, refM4, m4); + if(internal::random<bool>()) + m1.makeCompressed(); + VERIFY_IS_APPROX(m1*s1, refM1*s1); VERIFY_IS_APPROX(m1+m2, refM1+refM2); VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3); VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2)); VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2); - VERIFY_IS_APPROX(m1*=s1, refM1*=s1); - VERIFY_IS_APPROX(m1/=s1, refM1/=s1); - - VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); - VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); - if(SparseMatrixType::IsRowMajor) VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0))); else @@ -197,11 +194,29 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(refM4 - m3, refM4 - refM3); VERIFY_IS_APPROX(m3 - refM4, refM3 - refM4); + VERIFY_IS_APPROX(m1.sum(), refM1.sum()); + + VERIFY_IS_APPROX(m1*=s1, refM1*=s1); + VERIFY_IS_APPROX(m1/=s1, refM1/=s1); + + VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); + VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); + // test aliasing VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1)); VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval())); VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval())); VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1)); + + if(m1.isCompressed()) + { + VERIFY_IS_APPROX(m1.coeffs().sum(), m1.sum()); + m1.coeffs() += s1; + for(Index j = 0; j<m1.outerSize(); ++j) + for(typename SparseMatrixType::InnerIterator it(m1,j); it; ++it) + refM1(it.row(), it.col()) += s1; + VERIFY_IS_APPROX(m1, refM1); + } } // test transpose @@ -232,11 +247,11 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re for (Index i=0; i<m2.rows(); ++i) { float x = internal::random<float>(0,1); - if (x<0.1) + if (x<0.1f) { // do nothing } - else if (x<0.5) + else if (x<0.5f) { countFalseNonZero++; m2.insert(i,j) = Scalar(0); @@ -312,6 +327,17 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3); VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3); } + + Index i = internal::random<Index>(0,rows-1); + Index j = internal::random<Index>(0,cols-1); + m2.coeffRef(i,j) = 123; + if(internal::random<bool>()) + m2.makeCompressed(); + Map<SparseMatrixType> mapMat2(rows, cols, m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr()); + VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(123)); + VERIFY_IS_EQUAL(mapMat2.coeff(i,j),Scalar(123)); + mapMat2.coeffRef(i,j) = -123; + VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(-123)); } // test triangularView @@ -372,6 +398,12 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re SparseMatrixType m2(rows, rows); initSparse<Scalar>(density, refMat2, m2); VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval()); + + // sparse view on expressions: + VERIFY_IS_APPROX((s1*m2).eval(), (s1*refMat2).sparseView().eval()); + VERIFY_IS_APPROX((m2+m2).eval(), (refMat2+refMat2).sparseView().eval()); + VERIFY_IS_APPROX((m2*m2).eval(), (refMat2.lazyProduct(refMat2)).sparseView().eval()); + VERIFY_IS_APPROX((m2*m2).eval(), (refMat2*refMat2).sparseView().eval()); } // test diagonal @@ -546,7 +578,7 @@ void test_sparse_basic() CALL_SUBTEST_4((big_sparse_triplet<SparseMatrix<double, ColMajor, long int> >(10000, 10000, 0.125))); // Regression test for bug 1105 -#ifdef EIGEN_TEST_PART_6 +#ifdef EIGEN_TEST_PART_7 { int n = Eigen::internal::random<int>(200,600); SparseMatrix<std::complex<double>,0, long> mat(n, n); diff --git a/test/sparse_block.cpp b/test/sparse_block.cpp index 8a6e0687c..49a5f135e 100644 --- a/test/sparse_block.cpp +++ b/test/sparse_block.cpp @@ -17,6 +17,7 @@ template<typename SparseMatrixType> void sparse_block(const SparseMatrixType& re const Index outer = ref.outerSize(); typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::StorageIndex StorageIndex; double density = (std::max)(8./(rows*cols), 0.01); typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; @@ -123,7 +124,7 @@ template<typename SparseMatrixType> void sparse_block(const SparseMatrixType& re m3.reserve(VectorXi::Constant(outer,int(inner/2))); for(Index j=0; j<outer; ++j) for(Index k=0; k<(std::min)(j,inner); ++k) - m3.insertByOuterInner(j,k) = k+1; + m3.insertByOuterInner(j,k) = internal::convert_index<StorageIndex>(k+1); for(Index j=0; j<(std::min)(outer, inner); ++j) { VERIFY(j==numext::real(m3.innerVector(j).nonZeros())); @@ -150,7 +151,7 @@ template<typename SparseMatrixType> void sparse_block(const SparseMatrixType& re DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); SparseMatrixType m2(rows, cols); initSparse<Scalar>(density, refMat2, m2); - if(internal::random<float>(0,1)>0.5) m2.makeCompressed(); + if(internal::random<float>(0,1)>0.5f) m2.makeCompressed(); Index j0 = internal::random<Index>(0,outer-2); Index j1 = internal::random<Index>(0,outer-2); Index n0 = internal::random<Index>(1,outer-(std::max)(j0,j1)); diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 7ec5270e8..c7c93373d 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -7,8 +7,26 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +static long int nb_temporaries; + +inline void on_temporary_creation() { + // here's a great place to set a breakpoint when debugging failures in this test! + nb_temporaries++; +} + +#define EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN { on_temporary_creation(); } + #include "sparse.h" +#define VERIFY_EVALUATION_COUNT(XPR,N) {\ + nb_temporaries = 0; \ + CALL_SUBTEST( XPR ); \ + if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \ + VERIFY( (#XPR) && nb_temporaries==N ); \ + } + + + template<typename SparseMatrixType> void sparse_product() { typedef typename SparseMatrixType::StorageIndex StorageIndex; @@ -76,6 +94,24 @@ template<typename SparseMatrixType> void sparse_product() VERIFY_IS_APPROX(m4=(m2t.transpose()*m3t.transpose()).pruned(0), refMat4=refMat2t.transpose()*refMat3t.transpose()); VERIFY_IS_APPROX(m4=(m2*m3t.transpose()).pruned(0), refMat4=refMat2*refMat3t.transpose()); + // make sure the right product implementation is called: + if((!SparseMatrixType::IsRowMajor) && m2.rows()<=m3.cols()) + { + VERIFY_EVALUATION_COUNT(m4 = m2*m3, 3); // 1 temp for the result + 2 for transposing and get a sorted result. + VERIFY_EVALUATION_COUNT(m4 = (m2*m3).pruned(0), 1); + VERIFY_EVALUATION_COUNT(m4 = (m2*m3).eval().pruned(0), 4); + } + + // and that pruning is effective: + { + DenseMatrix Ad(2,2); + Ad << -1, 1, 1, 1; + SparseMatrixType As(Ad.sparseView()), B(2,2); + VERIFY_IS_EQUAL( (As*As.transpose()).eval().nonZeros(), 4); + VERIFY_IS_EQUAL( (Ad*Ad.transpose()).eval().sparseView().eval().nonZeros(), 2); + VERIFY_IS_EQUAL( (As*As.transpose()).pruned(1e-6).eval().nonZeros(), 2); + } + // dense ?= sparse * sparse VERIFY_IS_APPROX(dm4 =m2*m3, refMat4 =refMat2*refMat3); VERIFY_IS_APPROX(dm4+=m2*m3, refMat4+=refMat2*refMat3); @@ -245,7 +281,7 @@ template<typename SparseMatrixType> void sparse_product() for (int k=0; k<mS.outerSize(); ++k) for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it) if (it.index() == k) - it.valueRef() *= 0.5; + it.valueRef() *= Scalar(0.5); VERIFY_IS_APPROX(refS.adjoint(), refS); VERIFY_IS_APPROX(mS.adjoint(), mS); @@ -256,6 +292,10 @@ template<typename SparseMatrixType> void sparse_product() VERIFY_IS_APPROX(x=mUp.template selfadjointView<Upper>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mLo.template selfadjointView<Lower>()*b, refX=refS*b); VERIFY_IS_APPROX(x=mS.template selfadjointView<Upper|Lower>()*b, refX=refS*b); + + VERIFY_IS_APPROX(x.noalias()+=mUp.template selfadjointView<Upper>()*b, refX+=refS*b); + VERIFY_IS_APPROX(x.noalias()-=mLo.template selfadjointView<Lower>()*b, refX-=refS*b); + VERIFY_IS_APPROX(x.noalias()+=mS.template selfadjointView<Upper|Lower>()*b, refX+=refS*b); // sparse selfadjointView with sparse matrices SparseMatrixType mSres(rows,rows); diff --git a/test/sparse_ref.cpp b/test/sparse_ref.cpp index f4aefbb48..5e9607234 100644 --- a/test/sparse_ref.cpp +++ b/test/sparse_ref.cpp @@ -87,8 +87,8 @@ void call_ref() VERIFY_EVALUATION_COUNT( call_ref_3(B, B), 1); VERIFY_EVALUATION_COUNT( call_ref_2(B.transpose(), B.transpose()), 0); VERIFY_EVALUATION_COUNT( call_ref_3(B.transpose(), B.transpose()), 0); - VERIFY_EVALUATION_COUNT( call_ref_2(A*A, AA), 1); - VERIFY_EVALUATION_COUNT( call_ref_3(A*A, AA), 1); + VERIFY_EVALUATION_COUNT( call_ref_2(A*A, AA), 3); + VERIFY_EVALUATION_COUNT( call_ref_3(A*A, AA), 3); VERIFY(!C.isCompressed()); VERIFY_EVALUATION_COUNT( call_ref_3(C, C), 1); diff --git a/test/sparse_solver.h b/test/sparse_solver.h index b67653496..fd6199f3e 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -11,6 +11,33 @@ #include <Eigen/SparseCore> #include <sstream> +template<typename Solver, typename Rhs, typename Guess,typename Result> +void solve_with_guess(IterativeSolverBase<Solver>& solver, const MatrixBase<Rhs>& b, const Guess& g, Result &x) { + if(internal::random<bool>()) + { + // With a temporary through evaluator<SolveWithGuess> + x = solver.derived().solveWithGuess(b,g) + Result::Zero(x.rows(), x.cols()); + } + else + { + // direct evaluation within x through Assignment<Result,SolveWithGuess> + x = solver.derived().solveWithGuess(b.derived(),g); + } +} + +template<typename Solver, typename Rhs, typename Guess,typename Result> +void solve_with_guess(SparseSolverBase<Solver>& solver, const MatrixBase<Rhs>& b, const Guess& , Result& x) { + if(internal::random<bool>()) + x = solver.derived().solve(b) + Result::Zero(x.rows(), x.cols()); + else + x = solver.derived().solve(b); +} + +template<typename Solver, typename Rhs, typename Guess,typename Result> +void solve_with_guess(SparseSolverBase<Solver>& solver, const SparseMatrixBase<Rhs>& b, const Guess& , Result& x) { + x = solver.derived().solve(b); +} + template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs> void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db) { @@ -37,6 +64,12 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, } VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); VERIFY(x.isApprox(refX,test_precision<Scalar>())); + + x.setZero(); + solve_with_guess(solver, b, x, x); + VERIFY(solver.info() == Success && "solving failed when using analyzePattern/factorize API"); + VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x.isApprox(refX,test_precision<Scalar>())); x.setZero(); // test the analyze/factorize API diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp index d95f301d5..b3e1dda25 100644 --- a/test/sparse_vector.cpp +++ b/test/sparse_vector.cpp @@ -12,7 +12,7 @@ template<typename Scalar,typename StorageIndex> void sparse_vector(int rows, int cols) { double densityMat = (std::max)(8./(rows*cols), 0.01); - double densityVec = (std::max)(8./float(rows), 0.1); + double densityVec = (std::max)(8./(rows), 0.1); typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; typedef Matrix<Scalar,Dynamic,1> DenseVector; typedef SparseVector<Scalar,0,StorageIndex> SparseVectorType; diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp index 50d1fcdf2..e8605fd21 100644 --- a/test/sparseqr.cpp +++ b/test/sparseqr.cpp @@ -54,7 +54,7 @@ template<typename Scalar> void test_sparseqr_scalar() b = dA * DenseVector::Random(A.cols()); solver.compute(A); - if(internal::random<float>(0,1)>0.5) + if(internal::random<float>(0,1)>0.5f) solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change. if (solver.info() != Success) { diff --git a/test/stdvector.cpp b/test/stdvector.cpp index 6e173c678..50cb3341d 100644 --- a/test/stdvector.cpp +++ b/test/stdvector.cpp @@ -34,7 +34,7 @@ void check_stdvector_matrix(const MatrixType& m) VERIFY_IS_APPROX(v[21], y); v.push_back(x); VERIFY_IS_APPROX(v[22], x); - VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(MatrixType)); + VERIFY((internal::UIntPtr)&(v[22]) == (internal::UIntPtr)&(v[21]) + sizeof(MatrixType)); // do a lot of push_back such that the vector gets internally resized // (with memory reallocation) @@ -69,7 +69,7 @@ void check_stdvector_transform(const TransformType&) VERIFY_IS_APPROX(v[21], y); v.push_back(x); VERIFY_IS_APPROX(v[22], x); - VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(TransformType)); + VERIFY((internal::UIntPtr)&(v[22]) == (internal::UIntPtr)&(v[21]) + sizeof(TransformType)); // do a lot of push_back such that the vector gets internally resized // (with memory reallocation) @@ -104,7 +104,7 @@ void check_stdvector_quaternion(const QuaternionType&) VERIFY_IS_APPROX(v[21], y); v.push_back(x); VERIFY_IS_APPROX(v[22], x); - VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(QuaternionType)); + VERIFY((internal::UIntPtr)&(v[22]) == (internal::UIntPtr)&(v[21]) + sizeof(QuaternionType)); // do a lot of push_back such that the vector gets internally resized // (with memory reallocation) diff --git a/test/stdvector_overload.cpp b/test/stdvector_overload.cpp index 736ff0ee7..959665954 100644 --- a/test/stdvector_overload.cpp +++ b/test/stdvector_overload.cpp @@ -48,7 +48,7 @@ void check_stdvector_matrix(const MatrixType& m) VERIFY_IS_APPROX(v[21], y); v.push_back(x); VERIFY_IS_APPROX(v[22], x); - VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(MatrixType)); + VERIFY((internal::UIntPtr)&(v[22]) == (internal::UIntPtr)&(v[21]) + sizeof(MatrixType)); // do a lot of push_back such that the vector gets internally resized // (with memory reallocation) @@ -83,7 +83,7 @@ void check_stdvector_transform(const TransformType&) VERIFY_IS_APPROX(v[21], y); v.push_back(x); VERIFY_IS_APPROX(v[22], x); - VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(TransformType)); + VERIFY((internal::UIntPtr)&(v[22]) == (internal::UIntPtr)&(v[21]) + sizeof(TransformType)); // do a lot of push_back such that the vector gets internally resized // (with memory reallocation) @@ -118,7 +118,7 @@ void check_stdvector_quaternion(const QuaternionType&) VERIFY_IS_APPROX(v[21], y); v.push_back(x); VERIFY_IS_APPROX(v[22], x); - VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(QuaternionType)); + VERIFY((internal::UIntPtr)&(v[22]) == (internal::UIntPtr)&(v[21]) + sizeof(QuaternionType)); // do a lot of push_back such that the vector gets internally resized // (with memory reallocation) diff --git a/test/svd_common.h b/test/svd_common.h index d8611b541..605d5dfef 100644 --- a/test/svd_common.h +++ b/test/svd_common.h @@ -42,9 +42,14 @@ void svd_check_full(const MatrixType& m, const SvdType& svd) MatrixUType u = svd.matrixU(); MatrixVType v = svd.matrixV(); RealScalar scaling = m.cwiseAbs().maxCoeff(); - if(scaling<=(std::numeric_limits<RealScalar>::min)()) - scaling = RealScalar(1); - VERIFY_IS_APPROX(m/scaling, u * (sigma/scaling) * v.adjoint()); + if(scaling<(std::numeric_limits<RealScalar>::min)()) + { + VERIFY(sigma.cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)()); + } + else + { + VERIFY_IS_APPROX(m/scaling, u * (sigma/scaling) * v.adjoint()); + } VERIFY_IS_UNITARY(u); VERIFY_IS_UNITARY(v); } @@ -141,14 +146,14 @@ void svd_least_square(const MatrixType& m, unsigned int computationOptions) using std::abs; SolutionType y(x); - y.row(k) = (1.+2*NumTraits<RealScalar>::epsilon())*x.row(k); + y.row(k) = (RealScalar(1)+2*NumTraits<RealScalar>::epsilon())*x.row(k); RealScalar residual_y = (m*y-rhs).norm(); VERIFY( test_isMuchSmallerThan(abs(residual_y-residual), rhs_norm) || residual < residual_y ); if(internal::is_same<RealScalar,float>::value) ++g_test_level; VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); if(internal::is_same<RealScalar,float>::value) --g_test_level; - y.row(k) = (1.-2*NumTraits<RealScalar>::epsilon())*x.row(k); + y.row(k) = (RealScalar(1)-2*NumTraits<RealScalar>::epsilon())*x.row(k); residual_y = (m*y-rhs).norm(); VERIFY( test_isMuchSmallerThan(abs(residual_y-residual), rhs_norm) || residual < residual_y ); if(internal::is_same<RealScalar,float>::value) ++g_test_level; @@ -336,7 +341,7 @@ void svd_underoverflow() M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3)); svd.compute(M,ComputeFullU|ComputeFullV); CALL_SUBTEST( svd_check_full(M,svd) ); - + id(k)++; if(id(k)>=value_set.size()) { @@ -344,7 +349,7 @@ void svd_underoverflow() id.head(k).setZero(); k=0; } - + } while((id<int(value_set.size())).all()); #if defined __INTEL_COMPILER diff --git a/test/svd_fill.h b/test/svd_fill.h index 1bbe645ee..3877c0c7e 100644 --- a/test/svd_fill.h +++ b/test/svd_fill.h @@ -7,9 +7,20 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +template<typename T> +Array<T,4,1> four_denorms(); + +template<> +Array4f four_denorms() { return Array4f(5.60844e-39f, -5.60844e-39f, 4.94e-44f, -4.94e-44f); } +template<> +Array4d four_denorms() { return Array4d(5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324); } +template<typename T> +Array<T,4,1> four_denorms() { return four_denorms<double>().cast<T>(); } + template<typename MatrixType> void svd_fill_random(MatrixType &m, int Option = 0) { + using std::pow; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; @@ -18,7 +29,7 @@ void svd_fill_random(MatrixType &m, int Option = 0) s = internal::random<RealScalar>(1,s); Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); for(Index k=0; k<diagSize; ++k) - d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); + d(k) = d(k)*pow(RealScalar(10),internal::random<RealScalar>(-s,s)); bool dup = internal::random<int>(0,10) < 3; bool unit_uv = internal::random<int>(0,10) < (dup?7:3); // if we duplicate some diagonal entries, then increase the chance to preserve them using unitary U and V factors @@ -53,8 +64,9 @@ void svd_fill_random(MatrixType &m, int Option = 0) VT.setRandom(); } - Matrix<Scalar,Dynamic,1> samples(7); - samples << 0, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -1./NumTraits<RealScalar>::highest(), 1./NumTraits<RealScalar>::highest(); + Matrix<Scalar,Dynamic,1> samples(9); + samples << 0, four_denorms<RealScalar>(), + -RealScalar(1)/NumTraits<RealScalar>::highest(), RealScalar(1)/NumTraits<RealScalar>::highest(), (std::numeric_limits<RealScalar>::min)(), pow((std::numeric_limits<RealScalar>::min)(),0.8); if(Option==Symmetric) { diff --git a/test/triangular.cpp b/test/triangular.cpp index 936c2aef3..b96856486 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -65,7 +65,7 @@ template<typename MatrixType> void triangular_square(const MatrixType& m) m1 = MatrixType::Random(rows, cols); for (int i=0; i<rows; ++i) - while (numext::abs2(m1(i,i))<1e-1) m1(i,i) = internal::random<Scalar>(); + while (numext::abs2(m1(i,i))<RealScalar(1e-1)) m1(i,i) = internal::random<Scalar>(); Transpose<MatrixType> trm4(m4); // test back and forward subsitution with a vector as the rhs @@ -78,7 +78,7 @@ template<typename MatrixType> void triangular_square(const MatrixType& m) m3 = m1.template triangularView<Lower>(); VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps)); - // test back and forward subsitution with a matrix as the rhs + // test back and forward substitution with a matrix as the rhs m3 = m1.template triangularView<Upper>(); VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps)); m3 = m1.template triangularView<Lower>(); @@ -121,6 +121,14 @@ template<typename MatrixType> void triangular_square(const MatrixType& m) VERIFY_IS_APPROX(m1.template triangularView<Upper>() * m5, m3*m5); VERIFY_IS_APPROX(m6*m1.template triangularView<Upper>(), m6*m3); + m1up = m1.template triangularView<Upper>(); + VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up); + VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up); + VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint()); + VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint()); + + VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().diagonal(), m1.diagonal()); + } diff --git a/test/unalignedassert.cpp b/test/unalignedassert.cpp index e2f03ffca..731a08977 100644 --- a/test/unalignedassert.cpp +++ b/test/unalignedassert.cpp @@ -94,7 +94,7 @@ template<typename T> void construct_at_boundary(int boundary) { char buf[sizeof(T)+256]; - size_t _buf = reinterpret_cast<size_t>(buf); + size_t _buf = reinterpret_cast<internal::UIntPtr>(buf); _buf += (EIGEN_MAX_ALIGN_BYTES - (_buf % EIGEN_MAX_ALIGN_BYTES)); // make 16/32/...-byte aligned _buf += boundary; // make exact boundary-aligned T *x = ::new(reinterpret_cast<void*>(_buf)) T; diff --git a/test/vectorization_logic.cpp b/test/vectorization_logic.cpp index ee446c3c1..83c1439ad 100644 --- a/test/vectorization_logic.cpp +++ b/test/vectorization_logic.cpp @@ -7,6 +7,14 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#ifdef EIGEN_TEST_PART_1 +#define EIGEN_UNALIGNED_VECTORIZE 1 +#endif + +#ifdef EIGEN_TEST_PART_2 +#define EIGEN_UNALIGNED_VECTORIZE 0 +#endif + #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR #undef EIGEN_DEFAULT_TO_ROW_MAJOR #endif @@ -21,7 +29,7 @@ using internal::demangle_unrolling; template<typename Dst, typename Src> bool test_assign(const Dst&, const Src&, int traversal, int unrolling) { - typedef internal::copy_using_evaluator_traits<internal::evaluator<Dst>,internal::evaluator<Src>, internal::assign_op<typename Dst::Scalar> > traits; + typedef internal::copy_using_evaluator_traits<internal::evaluator<Dst>,internal::evaluator<Src>, internal::assign_op<typename Dst::Scalar,typename Src::Scalar> > traits; bool res = traits::Traversal==traversal; if(unrolling==InnerUnrolling+CompleteUnrolling) res = res && (int(traits::Unrolling)==InnerUnrolling || int(traits::Unrolling)==CompleteUnrolling); @@ -45,7 +53,7 @@ bool test_assign(const Dst&, const Src&, int traversal, int unrolling) template<typename Dst, typename Src> bool test_assign(int traversal, int unrolling) { - typedef internal::copy_using_evaluator_traits<internal::evaluator<Dst>,internal::evaluator<Src>, internal::assign_op<typename Dst::Scalar> > traits; + typedef internal::copy_using_evaluator_traits<internal::evaluator<Dst>,internal::evaluator<Src>, internal::assign_op<typename Dst::Scalar,typename Src::Scalar> > traits; bool res = traits::Traversal==traversal && traits::Unrolling==unrolling; if(!res) { @@ -65,7 +73,8 @@ bool test_assign(int traversal, int unrolling) template<typename Xpr> bool test_redux(const Xpr&, int traversal, int unrolling) { - typedef internal::redux_traits<internal::scalar_sum_op<typename Xpr::Scalar>,internal::redux_evaluator<Xpr> > traits; + typedef typename Xpr::Scalar Scalar; + typedef internal::redux_traits<internal::scalar_sum_op<Scalar,Scalar>,internal::redux_evaluator<Xpr> > traits; bool res = traits::Traversal==traversal && traits::Unrolling==unrolling; if(!res) @@ -144,10 +153,16 @@ struct vectorization_logic InnerVectorizedTraversal,InnerUnrolling)); VERIFY(test_assign(Matrix44u(),Matrix44()+Matrix44(), - LinearTraversal,NoUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? InnerVectorizedTraversal : LinearTraversal, + EIGEN_UNALIGNED_VECTORIZE ? InnerUnrolling : NoUnrolling)); + + VERIFY(test_assign(Matrix1(),Matrix1()+Matrix1(), + (Matrix1::InnerSizeAtCompileTime % PacketSize)==0 ? InnerVectorizedTraversal : LinearVectorizedTraversal, + CompleteUnrolling)); VERIFY(test_assign(Matrix1u(),Matrix1()+Matrix1(), - LinearTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? ((Matrix1::InnerSizeAtCompileTime % PacketSize)==0 ? InnerVectorizedTraversal : LinearVectorizedTraversal) + : LinearTraversal, CompleteUnrolling)); VERIFY(test_assign(Matrix44c().col(1),Matrix44c().col(2)+Matrix44c().col(3), InnerVectorizedTraversal,CompleteUnrolling)); @@ -158,19 +173,29 @@ struct vectorization_logic if(PacketSize>1) { typedef Matrix<Scalar,3,3,ColMajor> Matrix33c; + typedef Matrix<Scalar,3,1,ColMajor> Vector3; VERIFY(test_assign(Matrix33c().row(2),Matrix33c().row(1)+Matrix33c().row(1), LinearTraversal,CompleteUnrolling)); + VERIFY(test_assign(Vector3(),Vector3()+Vector3(), + EIGEN_UNALIGNED_VECTORIZE ? (HalfPacketSize==1 ? InnerVectorizedTraversal : LinearVectorizedTraversal) : (HalfPacketSize==1 ? InnerVectorizedTraversal : LinearTraversal), CompleteUnrolling)); VERIFY(test_assign(Matrix33c().col(0),Matrix33c().col(1)+Matrix33c().col(1), - LinearTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? (HalfPacketSize==1 ? InnerVectorizedTraversal : LinearVectorizedTraversal) : (HalfPacketSize==1 ? SliceVectorizedTraversal : LinearTraversal), + ((!EIGEN_UNALIGNED_VECTORIZE) && HalfPacketSize==1) ? NoUnrolling : CompleteUnrolling)); VERIFY(test_assign(Matrix3(),Matrix3().cwiseProduct(Matrix3()), LinearVectorizedTraversal,CompleteUnrolling)); VERIFY(test_assign(Matrix<Scalar,17,17>(),Matrix<Scalar,17,17>()+Matrix<Scalar,17,17>(), - HalfPacketSize==1 ? InnerVectorizedTraversal : LinearTraversal,NoUnrolling)); + HalfPacketSize==1 ? InnerVectorizedTraversal : + EIGEN_UNALIGNED_VECTORIZE ? LinearVectorizedTraversal : + LinearTraversal, + NoUnrolling)); + + VERIFY(test_assign(Matrix11(), Matrix11()+Matrix11(),InnerVectorizedTraversal,CompleteUnrolling)); + VERIFY(test_assign(Matrix11(),Matrix<Scalar,17,17>().template block<PacketSize,PacketSize>(2,3)+Matrix<Scalar,17,17>().template block<PacketSize,PacketSize>(8,4), - DefaultTraversal,PacketSize>4?InnerUnrolling:CompleteUnrolling)); + (EIGEN_UNALIGNED_VECTORIZE) ? InnerVectorizedTraversal : DefaultTraversal, CompleteUnrolling|InnerUnrolling)); VERIFY(test_assign(Vector1(),Matrix11()*Vector1(), InnerVectorizedTraversal,CompleteUnrolling)); @@ -208,7 +233,7 @@ struct vectorization_logic VERIFY((test_assign< Map<Matrix<Scalar,EIGEN_PLAIN_ENUM_MAX(2,PacketSize),EIGEN_PLAIN_ENUM_MAX(2,PacketSize)>, AlignedMax, InnerStride<3*PacketSize> >, Matrix<Scalar,EIGEN_PLAIN_ENUM_MAX(2,PacketSize),EIGEN_PLAIN_ENUM_MAX(2,PacketSize)> - >(DefaultTraversal,CompleteUnrolling))); + >(DefaultTraversal,PacketSize>=8?InnerUnrolling:CompleteUnrolling))); VERIFY((test_assign(Matrix11(), Matrix<Scalar,PacketSize,EIGEN_PLAIN_ENUM_MIN(2,PacketSize)>()*Matrix<Scalar,EIGEN_PLAIN_ENUM_MIN(2,PacketSize),PacketSize>(), InnerVectorizedTraversal, CompleteUnrolling))); @@ -270,6 +295,12 @@ struct vectorization_logic_half InnerVectorizedTraversal,CompleteUnrolling)); VERIFY(test_assign(Vector1(),Vector1()+Vector1(), InnerVectorizedTraversal,CompleteUnrolling)); + VERIFY(test_assign(Vector1(),Vector1().template segment<PacketSize>(0).derived(), + EIGEN_UNALIGNED_VECTORIZE ? InnerVectorizedTraversal : LinearVectorizedTraversal,CompleteUnrolling)); + VERIFY(test_assign(Vector1(),Scalar(2.1)*Vector1()-Vector1(), + InnerVectorizedTraversal,CompleteUnrolling)); + VERIFY(test_assign(Vector1(),(Scalar(2.1)*Vector1().template segment<PacketSize>(0)-Vector1().template segment<PacketSize>(0)).derived(), + EIGEN_UNALIGNED_VECTORIZE ? InnerVectorizedTraversal : LinearVectorizedTraversal,CompleteUnrolling)); VERIFY(test_assign(Vector1(),Vector1().cwiseProduct(Vector1()), InnerVectorizedTraversal,CompleteUnrolling)); VERIFY(test_assign(Vector1(),Vector1().template cast<Scalar>(), @@ -287,10 +318,11 @@ struct vectorization_logic_half InnerVectorizedTraversal,InnerUnrolling)); VERIFY(test_assign(Matrix57u(),Matrix57()+Matrix57(), - LinearTraversal,NoUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? InnerVectorizedTraversal : LinearTraversal, + EIGEN_UNALIGNED_VECTORIZE ? InnerUnrolling : NoUnrolling)); VERIFY(test_assign(Matrix1u(),Matrix1()+Matrix1(), - LinearTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? ((Matrix1::InnerSizeAtCompileTime % PacketSize)==0 ? InnerVectorizedTraversal : LinearVectorizedTraversal) : LinearTraversal,CompleteUnrolling)); if(PacketSize>1) { @@ -298,16 +330,17 @@ struct vectorization_logic_half VERIFY(test_assign(Matrix33c().row(2),Matrix33c().row(1)+Matrix33c().row(1), LinearTraversal,CompleteUnrolling)); VERIFY(test_assign(Matrix33c().col(0),Matrix33c().col(1)+Matrix33c().col(1), - LinearTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? (PacketSize==1 ? InnerVectorizedTraversal : LinearVectorizedTraversal) : LinearTraversal,CompleteUnrolling)); VERIFY(test_assign(Matrix3(),Matrix3().cwiseQuotient(Matrix3()), PacketTraits::HasDiv ? LinearVectorizedTraversal : LinearTraversal,CompleteUnrolling)); VERIFY(test_assign(Matrix<Scalar,17,17>(),Matrix<Scalar,17,17>()+Matrix<Scalar,17,17>(), - LinearTraversal,NoUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? (PacketSize==1 ? InnerVectorizedTraversal : LinearVectorizedTraversal) : LinearTraversal, + NoUnrolling)); VERIFY(test_assign(Matrix11(),Matrix<Scalar,17,17>().template block<PacketSize,PacketSize>(2,3)+Matrix<Scalar,17,17>().template block<PacketSize,PacketSize>(8,4), - DefaultTraversal,PacketSize>4?InnerUnrolling:CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? InnerVectorizedTraversal : DefaultTraversal,PacketSize>4?InnerUnrolling:CompleteUnrolling)); VERIFY(test_assign(Vector1(),Matrix11()*Vector1(), InnerVectorizedTraversal,CompleteUnrolling)); @@ -337,7 +370,7 @@ struct vectorization_logic_half >(DefaultTraversal,CompleteUnrolling))); VERIFY((test_assign(Matrix57(), Matrix<Scalar,5*PacketSize,3>()*Matrix<Scalar,3,7>(), - InnerVectorizedTraversal, CompleteUnrolling))); + InnerVectorizedTraversal, InnerUnrolling|CompleteUnrolling))); #endif } }; @@ -367,19 +400,19 @@ void test_vectorization_logic() if(internal::packet_traits<float>::Vectorizable) { VERIFY(test_assign(Matrix<float,3,3>(),Matrix<float,3,3>()+Matrix<float,3,3>(), - LinearTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? LinearVectorizedTraversal : LinearTraversal,CompleteUnrolling)); VERIFY(test_redux(Matrix<float,5,2>(), - DefaultTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? LinearVectorizedTraversal : DefaultTraversal,CompleteUnrolling)); } if(internal::packet_traits<double>::Vectorizable) { VERIFY(test_assign(Matrix<double,3,3>(),Matrix<double,3,3>()+Matrix<double,3,3>(), - LinearTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? LinearVectorizedTraversal : LinearTraversal,CompleteUnrolling)); VERIFY(test_redux(Matrix<double,7,3>(), - DefaultTraversal,CompleteUnrolling)); + EIGEN_UNALIGNED_VECTORIZE ? LinearVectorizedTraversal : DefaultTraversal,CompleteUnrolling)); } #endif // EIGEN_VECTORIZE diff --git a/test/vectorwiseop.cpp b/test/vectorwiseop.cpp index 3cc198772..739eacaf3 100644 --- a/test/vectorwiseop.cpp +++ b/test/vectorwiseop.cpp @@ -233,10 +233,10 @@ template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m) Matrix<Scalar,1,MatrixType::RowsAtCompileTime> tmp(rows); VERIFY_EVALUATION_COUNT( tmp = (m1 * m1.transpose()).colwise().sum(), (MatrixType::RowsAtCompileTime==Dynamic ? 1 : 0)); - m2 = m1.rowwise() - (m1.colwise().sum()/m1.rows()).eval(); - m1 = m1.rowwise() - (m1.colwise().sum()/m1.rows()); + m2 = m1.rowwise() - (m1.colwise().sum()/RealScalar(m1.rows())).eval(); + m1 = m1.rowwise() - (m1.colwise().sum()/RealScalar(m1.rows())); VERIFY_IS_APPROX( m1, m2 ); - VERIFY_EVALUATION_COUNT( m2 = (m1.rowwise() - m1.colwise().sum()/m1.rows()), (MatrixType::RowsAtCompileTime==Dynamic && MatrixType::ColsAtCompileTime!=1 ? 1 : 0) ); + VERIFY_EVALUATION_COUNT( m2 = (m1.rowwise() - m1.colwise().sum()/RealScalar(m1.rows())), (MatrixType::RowsAtCompileTime==Dynamic && MatrixType::ColsAtCompileTime!=1 ? 1 : 0) ); } void test_vectorwiseop() |