From 3ed67cb0bb4af65fbf243df598604a8c7630bf7d Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Tue, 4 Oct 2016 14:22:56 -0700 Subject: [PATCH 1/4] Fix a bug in the implementation of Carmack's fast sqrt algorithm in Eigen (enabled by EIGEN_FAST_MATH), which causes the vectorized parts of the computation to return -0.0 instead of NaN for negative arguments. Benchmark speed in Giga-sqrts/s Intel(R) Xeon(R) CPU E5-1650 v3 @ 3.50GHz ----------------------------------------- SSE AVX Fast=1 2.529G 4.380G Fast=0 1.944G 1.898G Fast=1 fixed 2.214G 3.739G This table illustrates the worst case in terms speed impact: It was measured by repeatedly computing the sqrt of an n=4096 float vector that fits in L1 cache. For large vectors the operation becomes memory bound and the differences between the different versions almost negligible. --- Eigen/src/Core/arch/AVX/MathFunctions.h | 24 +++++++------------ Eigen/src/Core/arch/SSE/MathFunctions.h | 19 ++++++++------- test/packetmath.cpp | 31 +++++++++++-------------- 3 files changed, 34 insertions(+), 40 deletions(-) diff --git a/Eigen/src/Core/arch/AVX/MathFunctions.h b/Eigen/src/Core/arch/AVX/MathFunctions.h index d21ec39dd..da70b6636 100644 --- a/Eigen/src/Core/arch/AVX/MathFunctions.h +++ b/Eigen/src/Core/arch/AVX/MathFunctions.h @@ -362,23 +362,17 @@ pexp(const Packet4d& _x) { template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f psqrt(const Packet8f& _x) { - _EIGEN_DECLARE_CONST_Packet8f(one_point_five, 1.5f); - _EIGEN_DECLARE_CONST_Packet8f(minus_half, -0.5f); - _EIGEN_DECLARE_CONST_Packet8f_FROM_INT(flt_min, 0x00800000); - - Packet8f neg_half = pmul(_x, p8f_minus_half); - - // select only the inverse sqrt of positive normal inputs (denormals are - // flushed to zero and cause infs as well). - Packet8f non_zero_mask = _mm256_cmp_ps(_x, p8f_flt_min, _CMP_GE_OQ); - Packet8f x = _mm256_and_ps(non_zero_mask, _mm256_rsqrt_ps(_x)); + Packet8f half = pmul(_x, pset1(.5f)); + Packet8f denormal_mask = _mm256_and_ps( + _mm256_cmpge_ps(_x, _mm256_setzero_ps()), + _mm256_cmplt_ps(_x, pset1((std::numeric_limits::min)()))); + // Compute approximate reciprocal sqrt. + Packet8f x = _mm256_rsqrt_ps(_x); // Do a single step of Newton's iteration. - x = pmul(x, pmadd(neg_half, pmul(x, x), p8f_one_point_five)); - - // Multiply the original _x by it's reciprocal square root to extract the - // square root. - return pmul(_x, x); + x = pmul(x, psub(pset1(1.5f), pmul(half, pmul(x,x)))); + // Flush results for denormals to zero. + return _mm256_andnot_ps(denormal_mask, pmul(_x,x)); } #else template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h index ac2fd8103..2c34a869d 100644 --- a/Eigen/src/Core/arch/SSE/MathFunctions.h +++ b/Eigen/src/Core/arch/SSE/MathFunctions.h @@ -32,7 +32,7 @@ Packet4f plog(const Packet4f& _x) /* the smallest non denormalized float number */ _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(min_norm_pos, 0x00800000); _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_inf, 0xff800000);//-1.f/0.f); - + /* natural logarithm computed for 4 simultaneous float return NaN for x <= 0 */ @@ -451,18 +451,21 @@ template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f psqrt(const Packet4f& _x) { Packet4f half = pmul(_x, pset1(.5f)); + Packet4f denormal_mask = _mm_and_ps( + _mm_cmpge_ps(_x, _mm_setzero_ps()), + _mm_cmplt_ps(_x, pset1((std::numeric_limits::min)()))); - /* select only the inverse sqrt of non-zero inputs */ - Packet4f non_zero_mask = _mm_cmpge_ps(_x, pset1((std::numeric_limits::min)())); - Packet4f x = _mm_and_ps(non_zero_mask, _mm_rsqrt_ps(_x)); - + // Compute approximate reciprocal sqrt. + Packet4f x = _mm_rsqrt_ps(_x); + // Do a single step of Newton's iteration. x = pmul(x, psub(pset1(1.5f), pmul(half, pmul(x,x)))); - return pmul(_x,x); + // Flush results for denormals to zero. + return _mm_andnot_ps(denormal_mask, pmul(_x,x)); } #else -template<>EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED +template<>EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f psqrt(const Packet4f& x) { return _mm_sqrt_ps(x); } #endif @@ -491,7 +494,7 @@ Packet4f prsqrt(const Packet4f& _x) { Packet4f neg_mask = _mm_cmplt_ps(_x, _mm_setzero_ps()); Packet4f zero_mask = _mm_andnot_ps(neg_mask, le_zero_mask); Packet4f infs_and_nans = _mm_or_ps(_mm_and_ps(neg_mask, p4f_nan), - _mm_and_ps(zero_mask, p4f_inf)); + _mm_and_ps(zero_mask, p4f_inf)); // Do a single step of Newton's iteration. x = pmul(x, pmadd(neg_half, pmul(x, x), p4f_one_point_five)); diff --git a/test/packetmath.cpp b/test/packetmath.cpp index 1394d9f2b..c3d3e1521 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp @@ -193,7 +193,7 @@ template void packetmath() internal::pstore(data2+3*PacketSize, A3); VERIFY(areApprox(ref, data2, 4*PacketSize) && "internal::pbroadcast4"); } - + { for (int i=0; i void packetmath() internal::pstore(data2+1*PacketSize, A1); VERIFY(areApprox(ref, data2, 2*PacketSize) && "internal::pbroadcast2"); } - + VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload(data1))) && "internal::pfirst"); - + if(PacketSize>1) { for(int offset=0;offset<4;++offset) @@ -315,7 +315,7 @@ template void packetmath_real() CHECK_CWISE1_IF(PacketTraits::HasRound, numext::round, internal::pround); CHECK_CWISE1_IF(PacketTraits::HasCeil, numext::ceil, internal::pceil); CHECK_CWISE1_IF(PacketTraits::HasFloor, numext::floor, internal::pfloor); - + for (int i=0; i(-1,1); @@ -440,12 +440,9 @@ template void packetmath_real() data1[0] = Scalar(-1.0f); h.store(data2, internal::plog(h.load(data1))); VERIFY((numext::isnan)(data2[0])); -#if !EIGEN_FAST_MATH h.store(data2, internal::psqrt(h.load(data1))); VERIFY((numext::isnan)(data2[0])); VERIFY((numext::isnan)(data2[1])); -#endif - } } @@ -459,7 +456,7 @@ template void packetmath_notcomplex() EIGEN_ALIGN_MAX Scalar data1[PacketTraits::size*4]; EIGEN_ALIGN_MAX Scalar data2[PacketTraits::size*4]; EIGEN_ALIGN_MAX Scalar ref[PacketTraits::size*4]; - + Array::Map(data1, PacketTraits::size*4).setRandom(); ref[0] = data1[0]; @@ -478,7 +475,7 @@ template void packetmath_notcomplex() for (int i=0; i(data1))) && "internal::predux_max"); - + for (int i=0; i(data1[0])); @@ -490,12 +487,12 @@ template void test_conj_helper(Scalar typedef internal::packet_traits PacketTraits; typedef typename PacketTraits::type Packet; const int PacketSize = PacketTraits::size; - + internal::conj_if cj0; internal::conj_if cj1; internal::conj_helper cj; internal::conj_helper pcj; - + for(int i=0;i void test_conj_helper(Scalar } internal::pstore(pval,pcj.pmul(internal::pload(data1),internal::pload(data2))); VERIFY(areApprox(ref, pval, PacketSize) && "conj_helper pmul"); - + for(int i=0;i void packetmath_complex() data1[i] = internal::random() * Scalar(1e2); data2[i] = internal::random() * Scalar(1e2); } - + test_conj_helper (data1,data2,ref,pval); test_conj_helper (data1,data2,ref,pval); test_conj_helper (data1,data2,ref,pval); test_conj_helper (data1,data2,ref,pval); - + { for(int i=0;i void packetmath_scatter_gather() for (int i=0; i()/RealScalar(PacketSize); } - + int stride = internal::random(1,20); - + EIGEN_ALIGN_MAX Scalar buffer[PacketSize*20]; memset(buffer, 0, 20*sizeof(Packet)); Packet packet = internal::pload(data1); @@ -594,7 +591,7 @@ void test_packetmath() CALL_SUBTEST_1( packetmath_notcomplex() ); CALL_SUBTEST_2( packetmath_notcomplex() ); CALL_SUBTEST_3( packetmath_notcomplex() ); - + CALL_SUBTEST_1( packetmath_real() ); CALL_SUBTEST_2( packetmath_real() ); From 765615609d40f485c718d60b76b35a277c36d5ce Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Tue, 4 Oct 2016 15:08:41 -0700 Subject: [PATCH 2/4] Update comment for fast sqrt. --- Eigen/src/Core/arch/SSE/MathFunctions.h | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h index 2c34a869d..7b5f948e1 100644 --- a/Eigen/src/Core/arch/SSE/MathFunctions.h +++ b/Eigen/src/Core/arch/SSE/MathFunctions.h @@ -444,9 +444,14 @@ Packet4f pcos(const Packet4f& _x) #if EIGEN_FAST_MATH -// This is based on Quake3's fast inverse square root. +// Functions for sqrt. +// The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step +// of Newton's method, at a cost of 1-2 bits of precision as opposed to the +// exact solution. It does not handle +inf, or denormalized numbers correctly. +// The main advantage of this approach is not just speed, but also the fact that +// it can be inlined and pipelined with other computations, further reducing its +// effective latency. This is similar to Quake3's fast inverse square root. // For detail see here: http://www.beyond3d.com/content/articles/8/ -// It lacks 1 (or 2 bits in some rare cases) of precision, and does not handle negative, +inf, or denormalized numbers correctly. template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f psqrt(const Packet4f& _x) { From 7f67e6dfdb635ed1e744c0794bb1f3dd4306678f Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Tue, 4 Oct 2016 15:09:11 -0700 Subject: [PATCH 3/4] Update comment for fast sqrt. --- Eigen/src/Core/arch/AVX/MathFunctions.h | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Core/arch/AVX/MathFunctions.h b/Eigen/src/Core/arch/AVX/MathFunctions.h index da70b6636..25cd992ef 100644 --- a/Eigen/src/Core/arch/AVX/MathFunctions.h +++ b/Eigen/src/Core/arch/AVX/MathFunctions.h @@ -355,9 +355,11 @@ pexp(const Packet4d& _x) { // Functions for sqrt. // The EIGEN_FAST_MATH version uses the _mm_rsqrt_ps approximation and one step // of Newton's method, at a cost of 1-2 bits of precision as opposed to the -// exact solution. The main advantage of this approach is not just speed, but -// also the fact that it can be inlined and pipelined with other computations, -// further reducing its effective latency. +// exact solution. It does not handle +inf, or denormalized numbers correctly. +// The main advantage of this approach is not just speed, but also the fact that +// it can be inlined and pipelined with other computations, further reducing its +// effective latency. This is similar to Quake3's fast inverse square root. +// For detail see here: http://www.beyond3d.com/content/articles/8/ #if EIGEN_FAST_MATH template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f From 47150af1c8e22174296d6b2077d719c5b786b1b0 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Wed, 12 Oct 2016 08:34:39 -0700 Subject: [PATCH 4/4] Fix copy-paste error: Must use _mm256_cmp_ps for AVX. --- Eigen/src/Core/arch/AVX/MathFunctions.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/arch/AVX/MathFunctions.h b/Eigen/src/Core/arch/AVX/MathFunctions.h index 25cd992ef..6af67ce2d 100644 --- a/Eigen/src/Core/arch/AVX/MathFunctions.h +++ b/Eigen/src/Core/arch/AVX/MathFunctions.h @@ -366,8 +366,9 @@ EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet8f psqrt(const Packet8f& _x) { Packet8f half = pmul(_x, pset1(.5f)); Packet8f denormal_mask = _mm256_and_ps( - _mm256_cmpge_ps(_x, _mm256_setzero_ps()), - _mm256_cmplt_ps(_x, pset1((std::numeric_limits::min)()))); + _mm256_cmp_ps(_x, pset1((std::numeric_limits::min)()), + _CMP_LT_OQ), + _mm256_cmp_ps(_x, _mm256_setzero_ps(), _CMP_GE_OQ)); // Compute approximate reciprocal sqrt. Packet8f x = _mm256_rsqrt_ps(_x);