Add reciprocal packet op and fast specializations for float with SSE, AVX, and AVX512.
diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h
index 7223428..4f1ff6b 100644
--- a/Eigen/src/Core/GenericPacketMath.h
+++ b/Eigen/src/Core/GenericPacketMath.h
@@ -65,6 +65,7 @@
HasCmp = 0,
HasDiv = 0,
+ HasReciprocal = 0,
HasSqrt = 0,
HasRsqrt = 0,
HasExp = 0,
@@ -816,13 +817,6 @@
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet psqrt(const Packet& a) { return numext::sqrt(a); }
-/** \internal \returns the reciprocal square-root of \a a (coeff-wise) */
-template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet prsqrt(const Packet& a) {
- typedef typename internal::unpacket_traits<Packet>::type Scalar;
- return pdiv(pset1<Packet>(Scalar(1)), psqrt(a));
-}
-
/** \internal \returns the rounded value of \a a (coeff-wise) */
template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
Packet pround(const Packet& a) { using numext::round; return round(a); }
@@ -1035,6 +1029,19 @@
return ifPacket.select[0] ? thenPacket : elsePacket;
}
+/** \internal \returns 1 / a (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC inline Packet preciprocal(const Packet& a) {
+ using Scalar = typename unpacket_traits<Packet>::type;
+ return pdiv(pset1<Packet>(Scalar(1)), a);
+}
+
+/** \internal \returns the reciprocal square-root of \a a (coeff-wise) */
+template<typename Packet> EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet prsqrt(const Packet& a) {
+ return preciprocal<Packet>(psqrt(a));
+}
+
} // end namespace internal
} // end namespace Eigen
diff --git a/Eigen/src/Core/MathFunctionsImpl.h b/Eigen/src/Core/MathFunctionsImpl.h
index 2c9bbb5..182dd37 100644
--- a/Eigen/src/Core/MathFunctionsImpl.h
+++ b/Eigen/src/Core/MathFunctionsImpl.h
@@ -17,6 +17,35 @@
namespace internal {
+/** \internal Fast reciprocal using Newton-Raphson's method.
+ We assume that the starting guess provided in approx_a_recip has at least
+ half the leading mantissa bits in the correct result, such that a single
+ Newton-Raphson step is sufficient to get within 1-2 ulps of the currect result.
+*/
+template <typename Packet, int Steps>
+struct generic_reciprocal_newton_step {
+ static_assert(Steps > 0, "Steps must be at least 1.");
+ EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Packet
+ run(const Packet& a, const Packet& approx_a_recip) {
+ using Scalar = typename unpacket_traits<Packet>::type;
+ const Packet two = pset1<Packet>(Scalar(2));
+ const Packet neg_a = pnegate(a);
+ // Refine the approximation using one Newton-Raphson step:
+ // x_{i} = x_{i-1} * (2 - a * x_{i-1})
+ const Packet x =
+ generic_reciprocal_newton_step<Packet,Steps - 1>::run(a, approx_a_recip);
+ return pmul(x, pmadd(neg_a, x, two));
+ }
+};
+
+template<typename Packet>
+struct generic_reciprocal_newton_step<Packet, 0> {
+ EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Packet
+ run(const Packet& /*unused*/, const Packet& approx_a_recip) {
+ return approx_a_recip;
+ }
+};
+
/** \internal \returns the hyperbolic tan of \a a (coeff-wise)
Doesn't do anything fancy, just a 13/6-degree rational interpolant which
is accurate up to a couple of ulps in the (approximate) range [-8, 8],
diff --git a/Eigen/src/Core/arch/AVX/MathFunctions.h b/Eigen/src/Core/arch/AVX/MathFunctions.h
index d517dff..9f61af6 100644
--- a/Eigen/src/Core/arch/AVX/MathFunctions.h
+++ b/Eigen/src/Core/arch/AVX/MathFunctions.h
@@ -166,19 +166,12 @@
return pselect<Packet8f>(not_normal_finite_mask, y_approx, y_newton);
}
-#else
-template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
-Packet8f prsqrt<Packet8f>(const Packet8f& _x) {
- EIGEN_DECLARE_CONST_Packet8f(one, 1.0f);
- return _mm256_div_ps(p8f_one, _mm256_sqrt_ps(_x));
+template<> EIGEN_STRONG_INLINE Packet8f preciprocal<Packet8f>(const Packet8f& a) {
+ return generic_reciprocal_newton_step<Packet8f, /*Steps=*/1>::run(a, _mm256_rcp_ps(a));
}
+
#endif
-template <> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
-Packet4d prsqrt<Packet4d>(const Packet4d& _x) {
- EIGEN_DECLARE_CONST_Packet4d(one, 1.0);
- return _mm256_div_pd(p4d_one, _mm256_sqrt_pd(_x));
-}
F16_PACKET_FUNCTION(Packet8f, Packet8h, psin)
F16_PACKET_FUNCTION(Packet8f, Packet8h, pcos)
@@ -190,6 +183,7 @@
F16_PACKET_FUNCTION(Packet8f, Packet8h, ptanh)
F16_PACKET_FUNCTION(Packet8f, Packet8h, psqrt)
F16_PACKET_FUNCTION(Packet8f, Packet8h, prsqrt)
+F16_PACKET_FUNCTION(Packet8f, Packet8h, preciprocal)
template <>
EIGEN_STRONG_INLINE Packet8h pfrexp(const Packet8h& a, Packet8h& exponent) {
@@ -214,6 +208,7 @@
BF16_PACKET_FUNCTION(Packet8f, Packet8bf, ptanh)
BF16_PACKET_FUNCTION(Packet8f, Packet8bf, psqrt)
BF16_PACKET_FUNCTION(Packet8f, Packet8bf, prsqrt)
+BF16_PACKET_FUNCTION(Packet8f, Packet8bf, preciprocal)
template <>
EIGEN_STRONG_INLINE Packet8bf pfrexp(const Packet8bf& a, Packet8bf& exponent) {
diff --git a/Eigen/src/Core/arch/AVX/PacketMath.h b/Eigen/src/Core/arch/AVX/PacketMath.h
index 6b20da6..bf832c9 100644
--- a/Eigen/src/Core/arch/AVX/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX/PacketMath.h
@@ -78,6 +78,7 @@
HasCmp = 1,
HasDiv = 1,
+ HasReciprocal = EIGEN_FAST_MATH,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
diff --git a/Eigen/src/Core/arch/AVX512/MathFunctions.h b/Eigen/src/Core/arch/AVX512/MathFunctions.h
index 54be6cf..6caed2d 100644
--- a/Eigen/src/Core/arch/AVX512/MathFunctions.h
+++ b/Eigen/src/Core/arch/AVX512/MathFunctions.h
@@ -253,13 +253,6 @@
// return rsqrt(+inf) = 0, rsqrt(x) = NaN if x < 0, and rsqrt(0) = +inf.
return _mm512_mask_blend_ps(not_finite_pos_mask, y_newton, y_approx);
}
-#else
-
-template <>
-EIGEN_STRONG_INLINE Packet16f prsqrt<Packet16f>(const Packet16f& x) {
- EIGEN_DECLARE_CONST_Packet16f(one, 1.0f);
- return _mm512_div_ps(p16f_one, _mm512_sqrt_ps(x));
-}
#endif
F16_PACKET_FUNCTION(Packet16f, Packet16h, prsqrt)
@@ -304,12 +297,17 @@
// return rsqrt(+inf) = 0, rsqrt(x) = NaN if x < 0, and rsqrt(0) = +inf.
return _mm512_mask_blend_pd(not_finite_pos_mask, y_newton, y_approx);
}
+
+template<> EIGEN_STRONG_INLINE Packet16f preciprocal<Packet16f>(const Packet16f& a) {
+#ifdef EIGEN_VECTORIZE_AVX512ER
+ return _mm512_rcp28_ps(a));
#else
-template <>
-EIGEN_STRONG_INLINE Packet8d prsqrt<Packet8d>(const Packet8d& x) {
- EIGEN_DECLARE_CONST_Packet8d(one, 1.0f);
- return _mm512_div_pd(p8d_one, _mm512_sqrt_pd(x));
+ return generic_reciprocal_newton_step<Packet16f, /*Steps=*/1>::run(a, _mm512_rcp14_ps(a));
+#endif
}
+
+F16_PACKET_FUNCTION(Packet16f, Packet16h, preciprocal)
+BF16_PACKET_FUNCTION(Packet16f, Packet16bf, preciprocal)
#endif
template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
diff --git a/Eigen/src/Core/arch/AVX512/PacketMath.h b/Eigen/src/Core/arch/AVX512/PacketMath.h
index ff0dd29..8a00c62 100644
--- a/Eigen/src/Core/arch/AVX512/PacketMath.h
+++ b/Eigen/src/Core/arch/AVX512/PacketMath.h
@@ -120,6 +120,7 @@
HasExp = 1,
HasSqrt = EIGEN_FAST_MATH,
HasRsqrt = EIGEN_FAST_MATH,
+ HasReciprocal = EIGEN_FAST_MATH,
HasTanh = EIGEN_FAST_MATH,
HasErf = EIGEN_FAST_MATH,
#endif
diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h
index 5a063d3..10bc107 100644
--- a/Eigen/src/Core/arch/SSE/MathFunctions.h
+++ b/Eigen/src/Core/arch/SSE/MathFunctions.h
@@ -148,21 +148,13 @@
return pselect<Packet4f>(not_normal_finite_mask, y_approx, y_newton);
}
-#else
-
-template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
-Packet4f prsqrt<Packet4f>(const Packet4f& x) {
- // Unfortunately we can't use the much faster mm_rsqrt_ps since it only provides an approximation.
- return _mm_div_ps(pset1<Packet4f>(1.0f), _mm_sqrt_ps(x));
-}
-
#endif
-template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED
-Packet2d prsqrt<Packet2d>(const Packet2d& x) {
- return _mm_div_pd(pset1<Packet2d>(1.0), _mm_sqrt_pd(x));
+template<> EIGEN_STRONG_INLINE Packet4f preciprocal<Packet4f>(const Packet4f& a) {
+ return generic_reciprocal_newton_step<Packet4f, /*Steps=*/1>::run(a, _mm_rcp_ps(a));
}
+
// Hyperbolic Tangent function.
template <>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet4f
diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h
index a843226..4de3d47 100755
--- a/Eigen/src/Core/arch/SSE/PacketMath.h
+++ b/Eigen/src/Core/arch/SSE/PacketMath.h
@@ -136,6 +136,7 @@
HasCmp = 1,
HasDiv = 1,
+ HasReciprocal = EIGEN_FAST_MATH,
HasSin = EIGEN_FAST_MATH,
HasCos = EIGEN_FAST_MATH,
HasLog = 1,
diff --git a/Eigen/src/Core/functors/UnaryFunctors.h b/Eigen/src/Core/functors/UnaryFunctors.h
index d56aae5..e7bccaf 100644
--- a/Eigen/src/Core/functors/UnaryFunctors.h
+++ b/Eigen/src/Core/functors/UnaryFunctors.h
@@ -723,13 +723,18 @@
EIGEN_DEVICE_FUNC inline Scalar operator() (const Scalar& a) const { return Scalar(1)/a; }
template<typename Packet>
EIGEN_DEVICE_FUNC inline const Packet packetOp(const Packet& a) const
- { return internal::pdiv(pset1<Packet>(Scalar(1)),a); }
+ { return internal::preciprocal(a); }
};
template <typename Scalar>
struct functor_traits<scalar_inverse_op<Scalar> > {
enum {
PacketAccess = packet_traits<Scalar>::HasDiv,
- Cost = scalar_div_cost<Scalar, PacketAccess>::value
+ // If packet_traits<Scalar>::HasReciprocal then the Estimated cost is that
+ // of computing an approximation plus a single Newton-Raphson step, which
+ // consists of 1 pmul + 1 pmadd.
+ Cost = (packet_traits<Scalar>::HasReciprocal
+ ? 4 * NumTraits<Scalar>::MulCost
+ : scalar_div_cost<Scalar, PacketAccess>::value)
};
};
@@ -1033,7 +1038,7 @@
}
};
-// TODO(rmlarsen): Enable the following on host when integer_packet is defined
+// TODO(rmlarsen): Enable the following on host when integer_packet is defined
// for the relevant packet types.
#ifdef EIGEN_GPU_CC
diff --git a/Eigen/src/LU/arch/InverseSize4.h b/Eigen/src/LU/arch/InverseSize4.h
index 220da03..69f5e79 100644
--- a/Eigen/src/LU/arch/InverseSize4.h
+++ b/Eigen/src/LU/arch/InverseSize4.h
@@ -58,10 +58,10 @@
const float* data = matrix.data();
const Index stride = matrix.innerStride();
- Packet4f L1_ = ploadt<Packet4f,MatrixAlignment>(data);
- Packet4f L2_ = ploadt<Packet4f,MatrixAlignment>(data + stride*4);
- Packet4f L3_ = ploadt<Packet4f,MatrixAlignment>(data + stride*8);
- Packet4f L4_ = ploadt<Packet4f,MatrixAlignment>(data + stride*12);
+ Packet4f L1 = ploadt<Packet4f,MatrixAlignment>(data);
+ Packet4f L2 = ploadt<Packet4f,MatrixAlignment>(data + stride*4);
+ Packet4f L3 = ploadt<Packet4f,MatrixAlignment>(data + stride*8);
+ Packet4f L4 = ploadt<Packet4f,MatrixAlignment>(data + stride*12);
// Four 2x2 sub-matrices of the input matrix
// input = [[A, B],
@@ -70,17 +70,17 @@
if (!StorageOrdersMatch)
{
- A = vec4f_unpacklo(L1_, L2_);
- B = vec4f_unpacklo(L3_, L4_);
- C = vec4f_unpackhi(L1_, L2_);
- D = vec4f_unpackhi(L3_, L4_);
+ A = vec4f_unpacklo(L1, L2);
+ B = vec4f_unpacklo(L3, L4);
+ C = vec4f_unpackhi(L1, L2);
+ D = vec4f_unpackhi(L3, L4);
}
else
{
- A = vec4f_movelh(L1_, L2_);
- B = vec4f_movehl(L2_, L1_);
- C = vec4f_movelh(L3_, L4_);
- D = vec4f_movehl(L4_, L3_);
+ A = vec4f_movelh(L1, L2);
+ B = vec4f_movehl(L2, L1);
+ C = vec4f_movelh(L3, L4);
+ D = vec4f_movehl(L4, L3);
}
Packet4f AB, DC;
@@ -120,7 +120,7 @@
Packet4f det = vec4f_duplane(psub(padd(d1, d2), d), 0);
// reciprocal of the determinant of the input matrix, rd = 1/det
- Packet4f rd = pdiv(pset1<Packet4f>(1.0f), det);
+ Packet4f rd = preciprocal(det);
// Four sub-matrices of the inverse
Packet4f iA, iB, iC, iD;
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index fe0a9f6..455ecab 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -28,6 +28,10 @@
return a / b;
}
template <typename T>
+inline T REF_RECIPROCAL(const T& a) {
+ return T(1) / a;
+}
+template <typename T>
inline T REF_ABS_DIFF(const T& a, const T& b) {
return a > b ? a - b : b - a;
}
@@ -464,9 +468,11 @@
CHECK_CWISE2_IF(PacketTraits::HasMul, REF_MUL, internal::pmul);
CHECK_CWISE2_IF(PacketTraits::HasDiv, REF_DIV, internal::pdiv);
- if (PacketTraits::HasNegate) CHECK_CWISE1(internal::negate, internal::pnegate);
+ CHECK_CWISE1_IF(PacketTraits::HasNegate, internal::negate, internal::pnegate);
+ CHECK_CWISE1_IF(PacketTraits::HasReciprocal, REF_RECIPROCAL, internal::preciprocal);
CHECK_CWISE1(numext::conj, internal::pconj);
+
for (int offset = 0; offset < 3; ++offset) {
for (int i = 0; i < PacketSize; ++i) ref[i] = data1[offset];
internal::pstore(data2, internal::pset1<Packet>(data1[offset]));
diff --git a/test/prec_inverse_4x4.cpp b/test/prec_inverse_4x4.cpp
index 86f0571..3eb061d 100644
--- a/test/prec_inverse_4x4.cpp
+++ b/test/prec_inverse_4x4.cpp
@@ -19,9 +19,7 @@
{
MatrixType m = PermutationMatrix<4>(indices);
MatrixType inv = m.inverse();
- double error = double( (m*inv-MatrixType::Identity()).norm() / NumTraits<Scalar>::epsilon() );
- EIGEN_DEBUG_VAR(error)
- VERIFY(error == 0.0);
+ VERIFY_IS_APPROX(m*inv, MatrixType::Identity());
std::next_permutation(indices.data(),indices.data()+4);
}
}
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
index e618042..c1609f1 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
@@ -601,13 +601,12 @@
ScalarType(6.79019408009981274425e-9)
};
const T eight = pset1<T>(ScalarType(8.0));
- const T one = pset1<T>(ScalarType(1));
const T neg_two = pset1<T>(ScalarType(-2));
T x, x0, x1, z;
x = psqrt(pmul(neg_two, plog(b)));
x0 = psub(x, pdiv(plog(x), x));
- z = pdiv(one, x);
+ z = preciprocal(x);
x1 = pmul(
z, pselect(
pcmp_lt(x, eight),
diff --git a/unsupported/test/cxx11_tensor_expr.cpp b/unsupported/test/cxx11_tensor_expr.cpp
index 27c2845..ddc8132 100644
--- a/unsupported/test/cxx11_tensor_expr.cpp
+++ b/unsupported/test/cxx11_tensor_expr.cpp
@@ -130,7 +130,7 @@
Tensor<float, 3, RowMajor> mat4(2,3,7);
mat4 = mat2 * 3.14f;
Tensor<float, 3> mat5(2,3,7);
- mat5 = mat1.inverse().log();
+ mat5 = (mat1 + mat1.constant(1)).inverse().log();
Tensor<float, 3, RowMajor> mat6(2,3,7);
mat6 = mat2.pow(0.5f) * 3.14f;
Tensor<float, 3> mat7(2,3,7);
@@ -150,7 +150,7 @@
for (int k = 0; k < 7; ++k) {
VERIFY_IS_APPROX(mat3(i,j,k), val + val);
VERIFY_IS_APPROX(mat4(i,j,k), val * 3.14f);
- VERIFY_IS_APPROX(mat5(i,j,k), logf(1.0f/val));
+ VERIFY_IS_APPROX(mat5(i,j,k), logf(1.0f/(val + 1)));
VERIFY_IS_APPROX(mat6(i,j,k), sqrtf(val) * 3.14f);
VERIFY_IS_APPROX(mat7(i,j,k), expf((std::max)(val, mat5(i,j,k) * 2.0f)));
VERIFY_IS_APPROX(mat8(i,j,k), expf(-val) * 3.14f);
