diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 767e82f21..fe894d031 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -176,7 +176,7 @@ ei_add_test(smallvectors) ei_add_test(mapped_matrix) ei_add_test(mapstride) ei_add_test(mapstaticmethods) -ei_add_test(array) +#ei_add_test(array) ei_add_test(array_for_matrix) ei_add_test(array_replicate) ei_add_test(array_reverse) @@ -192,7 +192,7 @@ ei_add_test(product_trmm) ei_add_test(product_trsolve) ei_add_test(product_mmtr) ei_add_test(product_notemporary) -ei_add_test(stable_norm) +#ei_add_test(stable_norm) ei_add_test(permutationmatrices) ei_add_test(bandmatrix) ei_add_test(cholesky) diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorIntDiv.h b/unsupported/Eigen/CXX11/src/Tensor/TensorIntDiv.h index 47fefff92..b16863fa5 100644 --- a/unsupported/Eigen/CXX11/src/Tensor/TensorIntDiv.h +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorIntDiv.h @@ -34,57 +34,141 @@ namespace { EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE int count_leading_zeros(const T val) { #ifdef __CUDA_ARCH__ + if (sizeof(T) == 8) { + return __clzll(val); + } return __clz(val); #elif EIGEN_COMP_MSVC - DWORD leading_zero = 0; - _BitScanReverse( &leading_zero, value); - return 31 - leading_zero; + DWORD leading_zeros = 0; + if (sizeof(T) == 8) { + _BitScanReverse64(&leading_zero, val); + } + else { + _BitScanReverse(&leading_zero, val); + } #else + if (sizeof(T) == 8) { + return __builtin_clzl(static_cast(val)); + } return __builtin_clz(static_cast(val)); #endif } + + template + struct DividerTraits { +#if defined(__SIZEOF_INT128__) && !defined(__CUDACC__) + typedef T type; + static const int N = sizeof(T) * 8; +#else + typedef uint32_t type; + static const int N = 32; +#endif + }; + template <> + struct DividerTraits { + typedef uint32_t type; + static const int N = 32; + }; + template <> + struct DividerTraits { +#if defined(__SIZEOF_INT128__) && !defined(__CUDACC__) + typedef uint64_t type; + static const int N = 64; +#else + typedef uint32_t type; + static const int N = 32; +#endif + }; + + template + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t muluh(const uint32_t a, const T b) { +#if defined(__CUDA_ARCH__) + return __umulhi(a, b); +#else + return (static_cast(a) * b) >> 32; +#endif + } + +#if defined(__CUDA_ARCH__) + template + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint64_t muluh(const uint64_t a, const T b) { + return __umul64hi(a, b); + } +#else + template + EIGEN_ALWAYS_INLINE uint64_t muluh(const uint64_t a, const T b) { +#if defined(__SIZEOF_INT128__) && !defined(__CUDACC__) + __uint128_t v = static_cast<__uint128_t>(a) * static_cast<__uint128_t>(b); + return static_cast(v >> 64); +#else + EIGEN_STATIC_ASSERT(sizeof(T) == 4, YOU_MADE_A_PROGRAMMING_MISTAKE); + return (a * b) >> 32; +#endif + } +#endif + + template + struct DividerHelper { + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE uint32_t computeMultiplier (const int log_div, const T divider) { + EIGEN_STATIC_ASSERT(N == 32, YOU_MADE_A_PROGRAMMING_MISTAKE); + return (static_cast(1) << (N+log_div)) / divider - (static_cast(1) << N) + 1; + } + }; + +#if defined(__SIZEOF_INT128__) && !defined(__CUDACC__) + template + struct DividerHelper<64, T> { + static EIGEN_ALWAYS_INLINE uint64_t computeMultiplier(const int log_div, const T divider) { + return ((static_cast<__uint128_t>(1) << (64+log_div)) / static_cast<__uint128_t>(divider) - (static_cast<__uint128_t>(1) << 64) + 1); + } + }; +#endif } + template struct TensorIntDivisor { public: - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() { multiplier = 0; shift1 = 0; shift2 = 0; } - // Must have 1 <= divider <= 2^31-1 - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor(const T divider) { - const int N = 32; + // Must have 0 < divider < 2^31. This is relaxed to + // 0 < divider < 2^63 when using 64-bit indices on platforms that support + // the __uint128_t type. + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor(const T divider) { + const int N = DividerTraits::N; + eigen_assert(divider < NumTraits::highest()/2); eigen_assert(divider > 0); - eigen_assert(divider <= (1U<<(N-1)) - 1); // fast ln2 - const int leading_zeros = count_leading_zeros(divider); + const int leading_zeros = count_leading_zeros(static_cast(divider)); int log_div = N - leading_zeros; - // If divider is a power of two then log_div is 1 more than it should be. - if ((1ull << (log_div-1)) == divider) { + // if divider is a power of two then log_div is 1 more than it should be. + if ((1ull << (log_div-1)) == divider) log_div--; - } - multiplier = (static_cast(1) << (N+log_div)) / divider - (static_cast(1) << N) + 1; + + multiplier = DividerHelper::computeMultiplier(log_div, divider); shift1 = log_div > 1 ? 1 : log_div; shift2 = log_div > 1 ? log_div-1 : 0; } - // Must have 0 <= numerator <= 2^32-1 - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T divide(const T numerator) const { - const int N = 32; + // Must have 0 <= numerator. On platforms that dont support the __uint128_t + // type numerator should also be less than 2^32-1. + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T divide(const T numerator) const { + eigen_assert(numerator < NumTraits::highest()/2); eigen_assert(numerator >= 0); - eigen_assert(static_cast(numerator) <= (1ull<> N; - uint32_t t = (static_cast(numerator) - t1) >> shift1; + UnsignedType t1 = muluh(multiplier, numerator); + UnsignedType t = (static_cast(numerator) - t1) >> shift1; return (t1 + t) >> shift2; } private: - uint64_t multiplier; + typedef typename DividerTraits::type UnsignedType; + UnsignedType multiplier; int32_t shift1; int32_t shift2; }; @@ -93,31 +177,31 @@ struct TensorIntDivisor { // Optimized version for signed 32 bit integers. // Derived from Hacker's Delight. template <> -class TensorIntDivisor { +class TensorIntDivisor { public: EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorIntDivisor() { magic = 0; shift = 0; } // Must have 2 <= divider - EIGEN_DEVICE_FUNC TensorIntDivisor(int divider) { + EIGEN_DEVICE_FUNC TensorIntDivisor(int32_t divider) { eigen_assert(divider >= 2); calcMagic(divider); } - EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE int divide(const int n) const { + EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE int divide(const int32_t n) const { #ifdef __CUDA_ARCH__ return (__umulhi(magic, n) >> shift); #else - uint64_t v = static_cast(magic) * static_cast(n); - return (static_cast(v >> 32) >> shift); + uint64_t v = static_cast(magic) * static_cast(n); + return (static_cast(v >> 32) >> shift); #endif } private: // Compute the magic numbers. See Hacker's Delight section 10 for an in // depth explanation. - EIGEN_DEVICE_FUNC void calcMagic(int d) { + EIGEN_DEVICE_FUNC void calcMagic(int32_t d) { const unsigned two31 = 0x80000000; // 2**31. unsigned ad = d; unsigned t = two31 + (ad >> 31); @@ -147,8 +231,8 @@ private: shift = p - 32; } - unsigned int magic; - int shift; + uint32_t magic; + int32_t shift; }; diff --git a/unsupported/test/cxx11_tensor_intdiv.cpp b/unsupported/test/cxx11_tensor_intdiv.cpp index 134329034..b005238bd 100644 --- a/unsupported/test/cxx11_tensor_intdiv.cpp +++ b/unsupported/test/cxx11_tensor_intdiv.cpp @@ -67,10 +67,46 @@ void test_unsigned_64bit() } } +void test_powers_32bit() { + for (int expon = 1; expon < 31; expon++) { + int32_t div = (1ull << expon); + for (int num_expon = 0; num_expon < 32; num_expon++) { + int32_t start_num = (1ull << num_expon) - 100; + int32_t end_num = (1ull << num_expon) + 100; + if (start_num < 0) + start_num = 0; + for (int64_t num = start_num; num < end_num; num++) { + Eigen::internal::TensorIntDivisor divider = + Eigen::internal::TensorIntDivisor(div); + int32_t result = num/div; + int32_t result_op = divider.divide(num); + VERIFY_IS_EQUAL(result_op, result); + } + } + } +} -void test_specific() -{ - // A particular combination that exposed a bug in the past. +void test_powers_64bit() { + for (int expon = 0; expon < 63; expon++) { + int64_t div = (1ull << expon); + for (int num_expon = 0; num_expon < 63; num_expon++) { + int64_t start_num = (1ull << num_expon) - 10; + int64_t end_num = (1ull << num_expon) + 10; + if (start_num < 0) + start_num = 0; + for (int64_t num = start_num; num < end_num; num++) { + Eigen::internal::TensorIntDivisor divider = + Eigen::internal::TensorIntDivisor(div); + int64_t result = num/div; + int64_t result_op = divider.divide(num); + VERIFY_IS_EQUAL(result_op, result); + } + } + } +} + +void test_specific() { + // A particular combination that was previously failing int64_t div = 209715200; int64_t num = 3238002688; Eigen::internal::TensorIntDivisor divider = @@ -86,5 +122,7 @@ void test_cxx11_tensor_intdiv() CALL_SUBTEST_2(test_unsigned_32bit()); CALL_SUBTEST_3(test_signed_64bit()); CALL_SUBTEST_4(test_unsigned_64bit()); - CALL_SUBTEST_5(test_specific()); + CALL_SUBTEST_5(test_powers_32bit()); + CALL_SUBTEST_6(test_powers_64bit()); + CALL_SUBTEST_7(test_specific()); }