diff options
author | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-10-04 14:22:56 -0700 |
---|---|---|
committer | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-10-04 14:22:56 -0700 |
commit | 3ed67cb0bb4af65fbf243df598604a8c7630bf7d (patch) | |
tree | 02c61529c1a3edab6c9894f271100a7488cd7cdc /test | |
parent | 6af5ac7e2749bdea7a31323855ef3b4333b91c3e (diff) |
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.
Diffstat (limited to 'test')
-rw-r--r-- | test/packetmath.cpp | 31 |
1 files changed, 14 insertions, 17 deletions
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<typename Scalar> void packetmath() internal::pstore(data2+3*PacketSize, A3); VERIFY(areApprox(ref, data2, 4*PacketSize) && "internal::pbroadcast4"); } - + { for (int i=0; i<PacketSize*2; ++i) ref[i] = data1[i/PacketSize]; @@ -203,9 +203,9 @@ template<typename Scalar> 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<Packet>(data1))) && "internal::pfirst"); - + if(PacketSize>1) { for(int offset=0;offset<4;++offset) @@ -315,7 +315,7 @@ template<typename Scalar> 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<size; ++i) { data1[i] = internal::random<Scalar>(-1,1); @@ -440,12 +440,9 @@ template<typename Scalar> 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<typename Scalar> 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<Scalar,Dynamic,1>::Map(data1, PacketTraits::size*4).setRandom(); ref[0] = data1[0]; @@ -478,7 +475,7 @@ template<typename Scalar> void packetmath_notcomplex() for (int i=0; i<PacketSize; ++i) ref[0] = (std::max)(ref[0],data1[i]); VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max"); - + for (int i=0; i<PacketSize; ++i) ref[i] = data1[0]+Scalar(i); internal::pstore(data2, internal::plset<Packet>(data1[0])); @@ -490,12 +487,12 @@ template<typename Scalar,bool ConjLhs,bool ConjRhs> void test_conj_helper(Scalar typedef internal::packet_traits<Scalar> PacketTraits; typedef typename PacketTraits::type Packet; const int PacketSize = PacketTraits::size; - + internal::conj_if<ConjLhs> cj0; internal::conj_if<ConjRhs> cj1; internal::conj_helper<Scalar,Scalar,ConjLhs,ConjRhs> cj; internal::conj_helper<Packet,Packet,ConjLhs,ConjRhs> pcj; - + for(int i=0;i<PacketSize;++i) { ref[i] = cj0(data1[i]) * cj1(data2[i]); @@ -503,7 +500,7 @@ template<typename Scalar,bool ConjLhs,bool ConjRhs> void test_conj_helper(Scalar } internal::pstore(pval,pcj.pmul(internal::pload<Packet>(data1),internal::pload<Packet>(data2))); VERIFY(areApprox(ref, pval, PacketSize) && "conj_helper pmul"); - + for(int i=0;i<PacketSize;++i) { Scalar tmp = ref[i]; @@ -531,12 +528,12 @@ template<typename Scalar> void packetmath_complex() data1[i] = internal::random<Scalar>() * Scalar(1e2); data2[i] = internal::random<Scalar>() * Scalar(1e2); } - + test_conj_helper<Scalar,false,false> (data1,data2,ref,pval); test_conj_helper<Scalar,false,true> (data1,data2,ref,pval); test_conj_helper<Scalar,true,false> (data1,data2,ref,pval); test_conj_helper<Scalar,true,true> (data1,data2,ref,pval); - + { for(int i=0;i<PacketSize;++i) ref[i] = Scalar(std::imag(data1[i]),std::real(data1[i])); @@ -556,9 +553,9 @@ template<typename Scalar> void packetmath_scatter_gather() for (int i=0; i<PacketSize; ++i) { data1[i] = internal::random<Scalar>()/RealScalar(PacketSize); } - + int stride = internal::random<int>(1,20); - + EIGEN_ALIGN_MAX Scalar buffer[PacketSize*20]; memset(buffer, 0, 20*sizeof(Packet)); Packet packet = internal::pload<Packet>(data1); @@ -594,7 +591,7 @@ void test_packetmath() CALL_SUBTEST_1( packetmath_notcomplex<float>() ); CALL_SUBTEST_2( packetmath_notcomplex<double>() ); CALL_SUBTEST_3( packetmath_notcomplex<int>() ); - + CALL_SUBTEST_1( packetmath_real<float>() ); CALL_SUBTEST_2( packetmath_real<double>() ); |