Core: Fix fast-math trig reduction and optimization barrier spills libeigen/eigen!3389 Closes #3200 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 3efa06c..7d40a34 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h
@@ -72,7 +72,7 @@ // Check before dividing: an overflowing reciprocal can raise FE_OVERFLOW even if the two-step fallback is used. if (sqrtNorm < Accumulator(1) && invScale > Accumulator(NumTraits<Accumulator>::highest()) * sqrtNorm) return false; Accumulator localSqrtNorm = sqrtNorm; - EIGEN_OPTIMIZATION_BARRIER(localSqrtNorm) + pbarrier(localSqrtNorm); factor = invScale / localSqrtNorm; return factor >= stable_normalization_normal_min<Accumulator, Accumulator>::run(); }
diff --git a/Eigen/src/Core/GenericPacketMath.h b/Eigen/src/Core/GenericPacketMath.h index a913986..0be87b6 100644 --- a/Eigen/src/Core/GenericPacketMath.h +++ b/Eigen/src/Core/GenericPacketMath.h
@@ -321,6 +321,18 @@ } }; +// Make x opaque to the compiler without changing its value, where supported. +template <typename Packet> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void pbarrier(Packet& x) { + EIGEN_UNUSED_VARIABLE(x); + EIGEN_OPTIMIZATION_BARRIER(x) +} +// The wrapper is a class, which EIGEN_OPTIMIZATION_BARRIER passes through memory; its vector can stay in a register. +template <typename T, int unique_id> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void pbarrier(eigen_packet_wrapper<T, unique_id>& x) { + pbarrier(x.m_val); +} + /** \internal \returns a + b (coeff-wise) */ template <typename Packet> EIGEN_DEVICE_FUNC inline Packet padd(const Packet& a, const Packet& b) {
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index bbd4e46..eaa8920 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h
@@ -2229,6 +2229,12 @@ namespace internal { +// Two-product residual using FMA, including software FMA when hardware FMA is unavailable. +template <typename Scalar> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar twoprod_low_fma(const Scalar& x, const Scalar& y, const Scalar& xy) { + return numext::fma(x, y, Scalar(-xy)); +} + template <typename T> EIGEN_DEVICE_FUNC bool isfinite_impl(const std::complex<T>& x) { EIGEN_IF_CONSTEXPR ((std::is_floating_point<T>::value)) { @@ -2456,7 +2462,7 @@ const bool negative = b < RealScalar(0); const RealScalar abs_b = numext::abs(b); const RealScalar t_hi = abs_b * ln2_hi; - const RealScalar t_lo = numext::fma(abs_b, ln2_hi, -t_hi) + abs_b * ln2_lo; + const RealScalar t_lo = twoprod_low_fma(abs_b, ln2_hi, t_hi) + abs_b * ln2_lo; // |t_lo| <= ulp(t_hi) / 2, so the Taylor terms below are exact to rounding unless |t| exceeds about 2 / sqrt(eps). // The phase is within an ulp of b ln(2) while |t| < 1 / eps. RealScalar sin_lo = t_lo;
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h b/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h index 728957f..afc13de 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathComplex.h
@@ -404,7 +404,7 @@ if (predux_any(por(pcmp_lt(hi, z_abs), is_small))) { Packet z_scaled = z; // Keeps compilers from executing this branch speculatively on the common path. - EIGEN_OPTIMIZATION_BARRIER(z_scaled.v) + pbarrier(z_scaled.v); return psqrt_complex_impl(z_scaled, true); } return psqrt_complex_impl(z, false);
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h b/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h index 1e8dcf1..ab08126 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathDoubleWord.h
@@ -25,6 +25,22 @@ r = psub(x, n); } +// Knuth two-sum: x + y = s_hi + s_lo exactly, with no magnitude ordering required. +template <typename Packet> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twosum(const Packet& x, const Packet& y, Packet& s_hi, Packet& s_lo) { + s_hi = padd(x, y); + pbarrier(s_hi); + Packet y_part = psub(s_hi, x); + pbarrier(y_part); + Packet x_part = psub(s_hi, y_part); + pbarrier(x_part); + Packet x_err = psub(x, x_part); + pbarrier(x_err); + Packet y_err = psub(y, y_part); + pbarrier(y_err); + s_lo = padd(x_err, y_err); +} + // This function computes the sum {s_hi, s_lo}, such that x + y = s_hi + s_lo // holds exactly, and s_hi = fl(x+y), if |x| >= |y|. template <typename Packet> @@ -44,14 +60,14 @@ template <typename Scalar, std::enable_if_t<is_scalar<Scalar>::value, int> = 0> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar twoprod_low(const Scalar& x, const Scalar& y, const Scalar& xy) { // Error-free products require FMA even when EIGEN_SCALAR_MADD_USE_FMA disables fusion in scalar madd. - return numext::fma(x, y, Scalar(-xy)); + return twoprod_low_fma(x, y, xy); } // This function implements the extended precision product of // a pair of floating point numbers. Given {x, y}, it computes the pair // {p_hi, p_lo} such that x * y = p_hi + p_lo holds exactly and -// p_hi = fl(x * y). -template <typename Packet> +// p_hi = fl(x * y). PreventReassociation preserves the separately rounded non-FMA operations. +template <typename Packet, bool PreventReassociation = false> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) { p_hi = pmul(x, y); p_lo = twoprod_low(x, y, p_hi); @@ -79,10 +95,10 @@ template <typename Scalar, std::enable_if_t<has_fast_fma<Scalar>::value, int> = 0> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar twoprod_low(const Scalar& x, const Scalar& y, const Scalar& xy) { - return numext::fma(x, y, Scalar(-xy)); + return twoprod_low_fma(x, y, xy); } -template <typename Scalar, std::enable_if_t<has_fast_fma<Scalar>::value, int> = 0> +template <typename Scalar, bool PreventReassociation = false, std::enable_if_t<has_fast_fma<Scalar>::value, int> = 0> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Scalar& x, const Scalar& y, Scalar& p_hi, Scalar& p_lo) { p_hi = x * y; p_lo = twoprod_low(x, y, p_hi); @@ -93,49 +109,50 @@ // exactly and that half of the significand of x fits in x_hi. // This is Algorithm 3 from Jean-Michel Muller, "Elementary Functions", // 3rd edition, Birkh\"auser, 2016. -template <typename Packet> +template <typename Packet, bool PreventReassociation = false> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void veltkamp_splitting(const Packet& x, Packet& x_hi, Packet& x_lo) { using Scalar = typename unpacket_traits<Packet>::type; constexpr int shift = (NumTraits<Scalar>::digits() + 1) / 2; const Scalar shift_scale = Scalar(uint64_t(1) << shift); // Scalar constructor not necessarily constexpr. - const Packet gamma = pmul(pset1<Packet>(shift_scale + Scalar(1)), x); + Packet gamma = pmul(pset1<Packet>(shift_scale + Scalar(1)), x); + EIGEN_IF_CONSTEXPR (PreventReassociation) pbarrier(gamma); Packet rho = psub(x, gamma); + EIGEN_IF_CONSTEXPR (PreventReassociation) pbarrier(rho); x_hi = padd(rho, gamma); + EIGEN_IF_CONSTEXPR (PreventReassociation) pbarrier(x_hi); x_lo = psub(x, x_hi); } -// This function implements Dekker's algorithm for products x * y. -// Given floating point numbers {x, y} computes the pair -// {p_hi, p_lo} such that x * y = p_hi + p_lo holds exactly and -// p_hi = fl(x * y). -template <typename Packet, std::enable_if_t<!has_fast_fma<Packet>::value, int> = 0> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) { - Packet x_hi, x_lo, y_hi, y_lo; - veltkamp_splitting(x, x_hi, x_lo); - veltkamp_splitting(y, y_hi, y_lo); - - p_hi = pmul(x, y); - p_lo = pmsub(x_hi, y_hi, p_hi); - p_lo = pmadd(x_hi, y_lo, p_lo); - p_lo = pmadd(x_lo, y_hi, p_lo); - p_lo = pmadd(x_lo, y_lo, p_lo); -} - // A version of twoprod that takes x, y, and fl(x*y) as input and returns the p_lo such that -// x * y = xy + p_lo holds exactly. -template <typename Packet, std::enable_if_t<!has_fast_fma<Packet>::value, int> = 0> +// x * y = xy + p_lo holds exactly. PreventReassociation preserves the separately rounded operations. +template <typename Packet, bool PreventReassociation = false, std::enable_if_t<!has_fast_fma<Packet>::value, int> = 0> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet twoprod_low(const Packet& x, const Packet& y, const Packet& xy) { Packet x_hi, x_lo, y_hi, y_lo; - veltkamp_splitting(x, x_hi, x_lo); - veltkamp_splitting(y, y_hi, y_lo); + veltkamp_splitting<Packet, PreventReassociation>(x, x_hi, x_lo); + veltkamp_splitting<Packet, PreventReassociation>(y, y_hi, y_lo); Packet p_lo = pmsub(x_hi, y_hi, xy); + EIGEN_IF_CONSTEXPR (PreventReassociation) pbarrier(p_lo); p_lo = pmadd(x_hi, y_lo, p_lo); + EIGEN_IF_CONSTEXPR (PreventReassociation) pbarrier(p_lo); p_lo = pmadd(x_lo, y_hi, p_lo); + EIGEN_IF_CONSTEXPR (PreventReassociation) pbarrier(p_lo); p_lo = pmadd(x_lo, y_lo, p_lo); return p_lo; } +// This function implements Dekker's algorithm for products x * y. +// Given floating point numbers {x, y}, it computes {p_hi, p_lo} such that +// x * y = p_hi + p_lo exactly and p_hi = fl(x * y). +template <typename Packet, bool PreventReassociation = false, std::enable_if_t<!has_fast_fma<Packet>::value, int> = 0> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void twoprod(const Packet& x, const Packet& y, Packet& p_hi, Packet& p_lo) { + Packet hi = pmul(x, y); + EIGEN_IF_CONSTEXPR (PreventReassociation) pbarrier(hi); + const Packet lo = twoprod_low<Packet, PreventReassociation>(x, y, hi); + p_hi = hi; + p_lo = lo; +} + #endif // EIGEN_VECTORIZE_FMA // This function implements the addition of two double word numbers
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFrexpLdexp.h b/Eigen/src/Core/arch/Default/GenericPacketMathFrexpLdexp.h index 4516b05..a1a7311 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFrexpLdexp.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFrexpLdexp.h
@@ -93,9 +93,9 @@ template <typename Packet> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Packet pldexp_apply_factors(const Packet& a, const Packet& c1, const Packet& c2) { Packet out = pmul(a, c1); - EIGEN_OPTIMIZATION_BARRIER(out) + pbarrier(out); out = pmul(out, c1); - EIGEN_OPTIMIZATION_BARRIER(out) + pbarrier(out); return pmul(out, c2); }
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h index 53970a0..a08a1d4 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
@@ -654,7 +654,7 @@ if (!predux_any(pcmp_lt(pset1<Packet>(Scalar(max_exponent - 2)), pabs(n)))) return pldexp_fast(y, n); const Packet h = print(pmul(n, pset1<Packet>(Scalar(0.5)))); Packet out = pldexp_fast(y, h); - EIGEN_OPTIMIZATION_BARRIER(out) + pbarrier(out); return pldexp_fast(out, psub(n, h)); } @@ -828,7 +828,7 @@ Packet sign_a = pandnot(a, abs_a); Packet rint_a = padd(abs_a, cst_limit); // Don't compile-away addition and subtraction. - EIGEN_OPTIMIZATION_BARRIER(rint_a); + pbarrier(rint_a); rint_a = psub(rint_a, cst_limit); rint_a = por(rint_a, sign_a); // If greater than limit (or NaN), simply return a.
diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h b/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h index 0884542..a98b7c7 100644 --- a/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h +++ b/Eigen/src/Core/arch/Default/GenericPacketMathTrig.h
@@ -100,7 +100,7 @@ // Rounding trick to find nearest integer: Packet y_round = padd(y, cst_rounding_magic); - EIGEN_OPTIMIZATION_BARRIER(y_round) + pbarrier(y_round); PacketI y_int = preinterpret<PacketI>(y_round); // last 23 digits represent integer (if abs(x)<2^24) y = psub(y_round, cst_rounding_magic); // nearest integer to x * (2/pi) @@ -122,11 +122,11 @@ // and 2 ULP up to: constexpr float huge_th = (Func == TrigFunction::Sin) ? 25966.f : 18838.f; x = pmadd(y, pset1<Packet>(-1.5703125), x); // = 0xbfc90000 - EIGEN_OPTIMIZATION_BARRIER(x) + pbarrier(x); x = pmadd(y, pset1<Packet>(-0.000483989715576171875), x); // = 0xb9fdc000 - EIGEN_OPTIMIZATION_BARRIER(x) + pbarrier(x); x = pmadd(y, pset1<Packet>(1.62865035235881805419921875e-07), x); // = 0x342ee000 - EIGEN_OPTIMIZATION_BARRIER(x) + pbarrier(x); x = pmadd(y, pset1<Packet>(5.5644315544167710640977020375430583953857421875e-11), x); // = 0x2e74b9ee // For the record, the following set of coefficients maintain 2ULP up @@ -263,24 +263,21 @@ // Compute q * pio2_1 exactly as a double-word using Dekker's algorithm. Packet qp_hi, qp_lo; - twoprod(cst_pio2_1<Packet>(), q, qp_hi, qp_lo); + twoprod<Packet, true>(cst_pio2_1<Packet>(), q, qp_hi, qp_lo); // Error-free addition of x and qp_hi using Knuth's 2sum. // Returns t_hi + t_lo = x + qp_hi exactly, with t_hi = fl(x + qp_hi). - Packet t_hi = padd(x, qp_hi); - Packet v = psub(t_hi, x); - Packet t_lo = padd(psub(x, psub(t_hi, v)), psub(qp_hi, v)); + Packet t_hi, t_lo; + twosum(x, qp_hi, t_hi, t_lo); // Near a multiple of pi/2, t_hi + qp_lo cancels against q * pio2_2 too, so that product is formed exactly as well // and the two large terms are added by 2sum: only terms of order u * |q * pio2_2| are rounded. Packet qp2_hi, qp2_lo; - twoprod(cst_pio2_2<Packet>(), q, qp2_hi, qp2_lo); - Packet a_hi = padd(t_hi, qp_lo); - v = psub(a_hi, t_hi); - t_lo = padd(t_lo, padd(psub(t_hi, psub(a_hi, v)), psub(qp_lo, v))); - t_hi = padd(a_hi, qp2_hi); - v = psub(t_hi, a_hi); - t_lo = padd(t_lo, padd(psub(a_hi, psub(t_hi, v)), psub(qp2_hi, v))); + twoprod<Packet, true>(cst_pio2_2<Packet>(), q, qp2_hi, qp2_lo); + Packet a_hi, a_lo, b_lo; + twosum(t_hi, qp_lo, a_hi, a_lo); + twosum(a_hi, qp2_hi, t_hi, b_lo); + t_lo = padd(t_lo, padd(a_lo, b_lo)); // Accumulate the low part of the product and the remaining pi/2 term. t_lo = padd(t_lo, qp2_lo); @@ -319,7 +316,7 @@ if (EIGEN_PREDICT_FALSE(predux_any(pcmp_le(pset1<Packet>(small_th), x_abs)))) { // Medium path: use double-word product x * (2/pi) for precise quadrant computation. Packet prod_hi, prod_lo; - twoprod(x_abs, cst_2oPI_hi, prod_hi, prod_lo); + twoprod<Packet, true>(x_abs, cst_2oPI_hi, prod_hi, prod_lo); // Correction for 2/pi truncation: add x * lo(2/pi) prod_lo = pmadd(x_abs, cst_2oPI_lo, prod_lo); @@ -336,9 +333,26 @@ } else { // Small path: simple reduction with triple-double pi/2 split. Packet qval_noround = pmul(x_abs, cst_2oPI_hi); +#ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD q_int = pcast<Packet, PacketI>(padd(qval_noround, pset1<Packet>(0.5))); Packet q = pcast<PacketI, Packet>(q_int); s = trig_reduce_small_double(x_abs, q); +#else + // Only q_int & 3 is consumed below. In round-to-nearest mode, the bias leaves q in its low bits. + const Packet bias = pset1<Packet>(6755399441055744.0); // 1.5 * 2^52. + Packet biased_q = padd(qval_noround, bias); + pbarrier(biased_q); + q_int = preinterpret<PacketI>(biased_q); + const Packet q = psub(biased_q, bias); + // |x| < 15 implies q <= 10. The first two 49-bit parts give exact q*c products; + // the third leaves < 2^-155 error in pi/2. Generated with MPFR at 256 bits: + // c1 = trunc49(pi/2), c2 = trunc49(pi/2-c1), c3 = RN53(pi/2-c1-c2). + s = pmadd(q, pset1<Packet>(-1.5707963267948948), x_abs); + pbarrier(s); + s = pmadd(q, pset1<Packet>(-1.837589179357618e-15), s); + pbarrier(s); + s = pmadd(q, pset1<Packet>(-8.4784276603688996e-32), s); +#endif } Packet ss = pmul(s, s);
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index cff4832..9ba7a07 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h
@@ -1125,6 +1125,7 @@ // x: (SSE) Any SSE register. // (AArch64) Like w, but restricted to registers 0 to 15 inclusive. // v: (PowerPC) An Altivec vector register. +// (x86) Any EVEX-encodable SSE register, including registers 16-31. // wa:(PowerPC) A VSX register. // // "X" (uppercase) should work for all cases, though this seems to fail for @@ -1143,9 +1144,48 @@ // gcc for ppc64le, // both vectors and scalars. // -// Note that this is restricted to plain types - this will not work -// directly for std::complex<T>, Eigen::half, Eigen::bfloat16. For these, -// you will need to apply to the underlying POD type. +// GCC (including MinGW) and Clang on x86, Arm and PPC pass class, union, +// complex and array operands through memory. Elsewhere, apply the barrier to +// the underlying plain type of std::complex<T>, Eigen::half or Eigen::bfloat16. +#if (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_CLANG || EIGEN_COMP_MINGW) && \ + (EIGEN_ARCH_PPC || EIGEN_ARCH_ARM_OR_ARM64 || EIGEN_ARCH_i386_OR_x86_64) +namespace Eigen { +namespace internal { +// Macros.h precedes Meta.h and <type_traits>. Operands are classified without cv-qualifiers: __builtin_classify_type +// does not evaluate its operand, but GCC rejects a volatile one in a constant expression. +template <typename T> +struct optimization_barrier_unqualified { + using type = T; +}; +template <typename T> +struct optimization_barrier_unqualified<const T> { + using type = T; +}; +template <typename T> +struct optimization_barrier_unqualified<volatile T> { + using type = T; +}; +template <typename T> +struct optimization_barrier_unqualified<const volatile T> { + using type = T; +}; +// Aggregates go to memory. GCC rejects register and general constraints for a struct of three doubles or a double[8] +// on every target, and for an int[4] or a struct of two doubles on 32-bit Arm; on PPC, complex<float> comes back +// corrupted at -O3 even with "m" added as an alternative. A memory operand hides the object's own bytes, not the +// storage it points to. +template <typename T> +struct optimization_barrier_in_memory { + // 9 = complex, 12 = class, 13 = union. Arrays decay to pointers here and are matched below. + static constexpr int kTypeClass = __builtin_classify_type(*static_cast<T*>(nullptr)); + static constexpr bool value = kTypeClass == 9 || kTypeClass == 12 || kTypeClass == 13; +}; +template <typename T, decltype(sizeof(0)) N> +struct optimization_barrier_in_memory<T[N]> { + static constexpr bool value = true; +}; +} // namespace internal +} // namespace Eigen +#endif #if EIGEN_ARCH_PPC && EIGEN_COMP_GNUC_STRICT // These register alternatives are broken on Clang. Packet4f is loaded into a single register rather than a vector, // zeroing out some entries, and integer types generate a compile error. @@ -1158,9 +1198,8 @@ #endif namespace Eigen { namespace internal { -// Complex, class, and union operands go to memory: a class that must stay in memory is a compile error with the -// register constraints, and complex<float> comes back corrupted at -O3, even with "m" added as an alternative. The two -// constraints need separate instantiations because GCC rejects an impossible constraint even on a branch not taken. +// The two constraints need separate instantiations because GCC rejects an impossible constraint even on a branch not +// taken. template <bool InMemory> struct ppc_optimization_barrier_impl { template <typename T> @@ -1177,9 +1216,8 @@ }; template <typename T> EIGEN_ALWAYS_INLINE void ppc_optimization_barrier(T& x) { - // __builtin_classify_type does not evaluate its operand: 9 = complex, 12 = class, 13 = union. - constexpr int kTypeClass = __builtin_classify_type(*static_cast<T*>(nullptr)); - ppc_optimization_barrier_impl<kTypeClass == 9 || kTypeClass == 12 || kTypeClass == 13>::run(x); + using Operand = typename optimization_barrier_unqualified<T>::type; + ppc_optimization_barrier_impl<optimization_barrier_in_memory<Operand>::value>::run(x); } } // namespace internal } // namespace Eigen @@ -1189,12 +1227,121 @@ // Clang's PPC backend does not accept one register constraint covering all scalar and vector operands. In particular, // "wa" crashes the backend for scalar integers narrower than 64 bits. #define EIGEN_OPTIMIZATION_BARRIER(X) __asm__("" : "+m"(X)); +#elif (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_CLANG || EIGEN_COMP_MINGW) && \ + (EIGEN_ARCH_ARM_OR_ARM64 || EIGEN_ARCH_i386_OR_x86_64) +// Clang chooses memory for mixed general/vector constraints. Use a register-only constraint for float, double, +// and builtin vectors that fit in one enabled register, and memory for aggregates. Other operands retain the +// general/memory alternatives. +#if EIGEN_ARCH_ARM_OR_ARM64 && defined(__ARM_FP) +#define EIGEN_BARRIER_GENERAL_CONSTRAINT "+g,w" +#define EIGEN_BARRIER_REGISTER_CONSTRAINT "+w" +#elif EIGEN_ARCH_ARM_OR_ARM64 +// Arm without VFP or NEON: "w" will not compile, and no operand takes the register path. +#define EIGEN_BARRIER_GENERAL_CONSTRAINT "+g" +#define EIGEN_BARRIER_REGISTER_CONSTRAINT "+g" +#else +#define EIGEN_BARRIER_GENERAL_CONSTRAINT "+g,x" +#ifdef __AVX512F__ +// Unlike "x", "v" also allows registers 16-31. +#define EIGEN_BARRIER_REGISTER_CONSTRAINT "+v" +#else +#define EIGEN_BARRIER_REGISTER_CONSTRAINT "+x" +#endif +#endif +#if EIGEN_ARCH_ARM_OR_ARM64 && defined(__ARM_NEON) +#define EIGEN_BARRIER_VECTOR_BYTES 16 +#elif EIGEN_ARCH_i386_OR_x86_64 && defined(__AVX512F__) +#define EIGEN_BARRIER_VECTOR_BYTES 64 +#elif EIGEN_ARCH_i386_OR_x86_64 && defined(__AVX__) +#define EIGEN_BARRIER_VECTOR_BYTES 32 +#elif EIGEN_ARCH_i386_OR_x86_64 && defined(__SSE2__) +#define EIGEN_BARRIER_VECTOR_BYTES 16 +#else +#define EIGEN_BARRIER_VECTOR_BYTES 0 +#endif +// The SME kernels use the barrier in streaming functions with ZA state. GCC refuses to inline a helper without these +// attributes there, and Clang 23 calls it out of line, switching streaming mode and lazily saving ZA around each call. +// Clang 20 and later have __arm_agnostic; for GCC, __arm_streaming_compatible suffices. +#if EIGEN_ARCH_ARM64 && defined(__ARM_FEATURE_SME) +#define EIGEN_BARRIER_SME_ATTRIBUTES __arm_streaming_compatible +#if EIGEN_COMP_CLANG +#if !__is_identifier(__arm_agnostic) +#undef EIGEN_BARRIER_SME_ATTRIBUTES +#define EIGEN_BARRIER_SME_ATTRIBUTES __arm_streaming_compatible __arm_agnostic("sme_za_state") +#endif +#endif +#else +#define EIGEN_BARRIER_SME_ATTRIBUTES +#endif +namespace Eigen { +namespace internal { +// GCC before 14 classifies builtin vectors as -1; newer GCC and Clang use 19. nullptr_t also has category -1 and is +// excluded below. Scalars have other categories. Sizeless SVE types have no sizeof and keep the general constraint. +template <typename T, typename = void> +struct optimization_barrier_in_register { + static constexpr bool value = false; +}; +template <typename T> +struct optimization_barrier_in_register<T, decltype(void(sizeof(T)))> { + static constexpr int kTypeClass = __builtin_classify_type(*static_cast<T*>(nullptr)); + static constexpr bool value = (kTypeClass == -1 || kTypeClass == 19) && sizeof(T) <= EIGEN_BARRIER_VECTOR_BYTES; +}; +template <> +struct optimization_barrier_in_register<decltype(nullptr)> { + static constexpr bool value = false; +}; +#if (EIGEN_ARCH_ARM_OR_ARM64 && defined(__ARM_FP) && (__ARM_FP & 4)) || (EIGEN_ARCH_i386_OR_x86_64 && defined(__SSE2__)) +template <> +struct optimization_barrier_in_register<float> { + static constexpr bool value = true; +}; +#endif +#if (EIGEN_ARCH_ARM_OR_ARM64 && defined(__ARM_FP) && (__ARM_FP & 8)) || (EIGEN_ARCH_i386_OR_x86_64 && defined(__SSE2__)) +template <> +struct optimization_barrier_in_register<double> { + static constexpr bool value = true; +}; +#endif +enum class optimization_barrier_operand_kind { General, Register, Memory }; +template <optimization_barrier_operand_kind Kind> +struct optimization_barrier_impl { + template <typename T> + static EIGEN_ALWAYS_INLINE void run(T& x) EIGEN_BARRIER_SME_ATTRIBUTES { + __asm__("" : EIGEN_BARRIER_GENERAL_CONSTRAINT(x)); + } +}; +template <> +struct optimization_barrier_impl<optimization_barrier_operand_kind::Register> { + template <typename T> + static EIGEN_ALWAYS_INLINE void run(T& x) EIGEN_BARRIER_SME_ATTRIBUTES { + __asm__("" : EIGEN_BARRIER_REGISTER_CONSTRAINT(x)); + } +}; +template <> +struct optimization_barrier_impl<optimization_barrier_operand_kind::Memory> { + template <typename T> + static EIGEN_ALWAYS_INLINE void run(T& x) EIGEN_BARRIER_SME_ATTRIBUTES { + __asm__("" : "+m"(x)); + } +}; +template <typename T> +EIGEN_ALWAYS_INLINE void optimization_barrier(T& x) EIGEN_BARRIER_SME_ATTRIBUTES { + using Operand = typename optimization_barrier_unqualified<T>::type; + constexpr optimization_barrier_operand_kind kKind = + optimization_barrier_in_memory<Operand>::value ? optimization_barrier_operand_kind::Memory + : optimization_barrier_in_register<Operand>::value ? optimization_barrier_operand_kind::Register + : optimization_barrier_operand_kind::General; + optimization_barrier_impl<kKind>::run(x); +} +} // namespace internal +} // namespace Eigen +#undef EIGEN_BARRIER_GENERAL_CONSTRAINT +#undef EIGEN_BARRIER_REGISTER_CONSTRAINT +#undef EIGEN_BARRIER_VECTOR_BYTES +#undef EIGEN_BARRIER_SME_ATTRIBUTES +#define EIGEN_OPTIMIZATION_BARRIER(X) Eigen::internal::optimization_barrier(X); #elif EIGEN_ARCH_ARM_OR_ARM64 #ifdef __ARM_FP -// General, VFP or NEON. -// Clang doesn't like "r", -// error: non-trivial scalar-to-vector conversion, possible invalid -// constraint for vector typ #define EIGEN_OPTIMIZATION_BARRIER(X) __asm__("" : "+g,w"(X)); #else // Arm without VFP or NEON. @@ -1224,6 +1371,29 @@ #else #define EIGEN_OPTIMIZATION_BARRIER(X) #endif +#elif EIGEN_COMP_MSVC && defined(_M_FP_FAST) +// MSVC has no inline asm on x64 or Arm64, and /fp:fast folds across an empty barrier, e.g. x - q*c1 - q*c2 into +// x - q*(c1 + c2). Reloading x from an address OR-ed with a volatile zero hides its value from the optimizer at the +// cost of a store and a store-forwarded load. +namespace Eigen { +namespace internal { +template <typename Dummy = void> +struct msvc_optimization_barrier_zero { + static volatile decltype(sizeof(0)) value; +}; +template <typename Dummy> +volatile decltype(sizeof(0)) msvc_optimization_barrier_zero<Dummy>::value = 0; +template <typename T> +EIGEN_ALWAYS_INLINE void msvc_optimization_barrier(T& x) { + x = *reinterpret_cast<T*>(reinterpret_cast<decltype(sizeof(0))>(&x) | msvc_optimization_barrier_zero<>::value); +} +template <typename T, decltype(sizeof(0)) N> +EIGEN_ALWAYS_INLINE void msvc_optimization_barrier(T (&x)[N]) { + for (auto& element : x) msvc_optimization_barrier(element); +} +} // namespace internal +} // namespace Eigen +#define EIGEN_OPTIMIZATION_BARRIER(X) Eigen::internal::msvc_optimization_barrier(X); #else // Not implemented for other compilers. #define EIGEN_OPTIMIZATION_BARRIER(X) @@ -1234,8 +1404,7 @@ namespace internal { // A volatile asm that the optimizer must execute exactly where it is written, so a value passed through it is defined // only on its branch, and an operation that consumes it stays there instead of being evaluated speculatively. -// EIGEN_OPTIMIZATION_BARRIER is not volatile, and its "+g,w" and "+g,x" alternatives put a scalar float or double in -// memory on clang; the float and double overloads keep the value in a register. +// EIGEN_OPTIMIZATION_BARRIER is not volatile and can move across branches. template <typename T> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE void volatile_barrier(T& x) { EIGEN_UNUSED_VARIABLE(x);
diff --git a/Eigen/src/SVD/BDCSVDImpl.h b/Eigen/src/SVD/BDCSVDImpl.h index 5f0a7b7..c302eb9 100644 --- a/Eigen/src/SVD/BDCSVDImpl.h +++ b/Eigen/src/SVD/BDCSVDImpl.h
@@ -412,8 +412,8 @@ RealScalar secondQuotient = secondNumerator / secondDenominator; #if defined(__FAST_MATH__) || EIGEN_COMP_NVHPC // NVHPC does not expose a preprocessor macro for -fast, so retain the barriers in all NVHPC builds. - EIGEN_OPTIMIZATION_BARRIER(firstQuotient) - EIGEN_OPTIMIZATION_BARRIER(secondQuotient) + pbarrier(firstQuotient); + pbarrier(secondQuotient); #endif return firstQuotient * secondQuotient; } @@ -423,7 +423,7 @@ RealScalar numerator, RealScalar firstDenominator, RealScalar secondDenominator) { RealScalar firstQuotient = numerator / firstDenominator; #if defined(__FAST_MATH__) || EIGEN_COMP_NVHPC - EIGEN_OPTIMIZATION_BARRIER(firstQuotient) + pbarrier(firstQuotient); #endif return firstQuotient / secondDenominator; } @@ -656,7 +656,7 @@ // Unsafe FP reassociation may otherwise turn `mus + (shift - dk)` into // `(mus + shift) - dk`, losing `mus` when `shift` and `dk` cancel. RealScalar diff = shifts(lastIdx) - dk; - EIGEN_OPTIMIZATION_BARRIER(diff) + pbarrier(diff); RealScalar prod = (singVals(lastIdx) + dk) * (mus(lastIdx) + diff); for (Index l = 0; l < m; ++l) { @@ -671,7 +671,7 @@ } Index j = i < k ? i : perm(l - 1); diff = shifts(j) - dk; - EIGEN_OPTIMIZATION_BARRIER(diff) + pbarrier(diff); prod *= productOfQuotients(singVals(j) + dk, diag(i) + dk, mus(j) + diff, diag(i) - dk); } } @@ -706,9 +706,9 @@ const RealScalar diagonal = diag(i); const RealScalar z = zhat(i); RealScalar diff = diagonal - shift; - EIGEN_OPTIMIZATION_BARRIER(diff) + pbarrier(diff); diff -= mu; - EIGEN_OPTIMIZATION_BARRIER(diff) + pbarrier(diff); const RealScalar sum = diagonal + singularValue; U(i, k) = sequentialQuotient(z, diff, sum); if (m_compV && l > 0) V(i, k) = sequentialQuotient(diagonal * z, diff, sum);
diff --git a/benchmarks/Core/CMakeLists.txt b/benchmarks/Core/CMakeLists.txt index b02faf5..ae9b2df 100644 --- a/benchmarks/Core/CMakeLists.txt +++ b/benchmarks/Core/CMakeLists.txt
@@ -81,6 +81,8 @@ eigen_add_benchmark(bench_selfadjoint_l1norm bench_selfadjoint_l1norm.cpp) eigen_add_benchmark(bench_selfadjoint_permutation bench_selfadjoint_permutation.cpp) eigen_add_benchmark(bench_cwise_math bench_cwise_math.cpp) +eigen_add_benchmark(bench_trig_reduction bench_trig_reduction.cpp) +eigen_add_benchmark(bench_optimization_barrier bench_optimization_barrier.cpp) eigen_add_benchmark(bench_complex_sqrt bench_complex_sqrt.cpp) eigen_add_benchmark(bench_cast bench_cast.cpp) eigen_add_benchmark(bench_complex_log bench_complex_log.cpp)
diff --git a/benchmarks/Core/bench_optimization_barrier.cpp b/benchmarks/Core/bench_optimization_barrier.cpp new file mode 100644 index 0000000..eebc49a --- /dev/null +++ b/benchmarks/Core/bench_optimization_barrier.cpp
@@ -0,0 +1,60 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Core> + +using namespace Eigen; + +template <typename Scalar, typename Packet, bool TwoSum> +EIGEN_DONT_INLINE void barrier_kernel(const Scalar* x, const Scalar* y, Scalar* hi, Scalar* lo, Index size) { + constexpr int width = internal::unpacket_traits<Packet>::size; + for (Index i = 0; i < size; i += width) { + const Packet value = internal::ploadu<Packet>(x + i); + Packet high = value; + Packet low = internal::pset1<Packet>(Scalar(0)); + if (TwoSum) { + internal::twosum(value, internal::ploadu<Packet>(y + i), high, low); + } else { + EIGEN_OPTIMIZATION_BARRIER(high); + } + internal::pstoreu(hi + i, high); + internal::pstoreu(lo + i, low); + } +} + +template <typename Scalar, typename Packet, bool TwoSum> +static void BM_OptimizationBarrier(benchmark::State& state) { + const Index size = state.range(0); + const Scalar tiny = numext::ldexp(Scalar(1), -NumTraits<Scalar>::digits() - 1); + Array<Scalar, Dynamic, 1> x(size), y(size), hi(size), lo(size); + for (Index i = 0; i < size; ++i) { + x[i] = Scalar(1) + Scalar(i % 16) / Scalar(16); + y[i] = tiny; + } + barrier_kernel<Scalar, Packet, TwoSum>(x.data(), y.data(), hi.data(), lo.data(), size); + for (Index i = 0; i < size; ++i) { + if (hi[i] != x[i] || lo[i] != (TwoSum ? tiny : Scalar(0))) { + state.SkipWithError("optimization barrier changed the result"); + return; + } + } + for (auto _ : state) { + barrier_kernel<Scalar, Packet, TwoSum>(x.data(), y.data(), hi.data(), lo.data(), size); + benchmark::DoNotOptimize(hi.data()); + benchmark::DoNotOptimize(lo.data()); + benchmark::ClobberMemory(); + } + state.SetItemsProcessed(state.iterations() * size); +} + +using FloatPacket = internal::packet_traits<float>::type; +using DoublePacket = internal::packet_traits<double>::type; +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, float, float, false)->Name("MacroFloatScalar")->Arg(4096); +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, float, FloatPacket, false)->Name("MacroFloatPacket")->Arg(4096); +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, double, double, false)->Name("MacroDoubleScalar")->Arg(4096); +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, double, DoublePacket, false)->Name("MacroDoublePacket")->Arg(4096); +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, float, float, true)->Name("TwoSumFloatScalar")->Arg(4096); +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, float, FloatPacket, true)->Name("TwoSumFloatPacket")->Arg(4096); +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, double, double, true)->Name("TwoSumDoubleScalar")->Arg(4096); +BENCHMARK_TEMPLATE(BM_OptimizationBarrier, double, DoublePacket, true)->Name("TwoSumDoublePacket")->Arg(4096);
diff --git a/benchmarks/Core/bench_trig_reduction.cpp b/benchmarks/Core/bench_trig_reduction.cpp new file mode 100644 index 0000000..4a92abf --- /dev/null +++ b/benchmarks/Core/bench_trig_reduction.cpp
@@ -0,0 +1,66 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include <benchmark/benchmark.h> +#include <Eigen/Core> + +#include <cmath> +#include <random> + +using namespace Eigen; + +// Domains exercise the small reduction, double-word reduction, and libm fallback. +template <int Function> +static void BM_TrigReduction(benchmark::State& state) { + const Index size = state.range(0); + const int domain = static_cast<int>(state.range(1)); + const double lower[] = {-10.0, 15.0, 1e4, 1e15}; + const double upper[] = {10.0, 1e4, 1e13, 2e15}; + // MPFR 4.2.1 at 256 bits, rounded to binary64; checked outside the timed loop. + const double reference[][4] = { + {0.3, 0.29552020666133955, 0.95533648912560598, 0.30933624960962325}, + {1000.0, 0.82687954053200252, 0.56237907629070294, 1.4703241557027185}, + {1e12, -0.61123870237688949, 0.79144630185289022, -0.77230596813187613}, + {1e15, 0.85827279317023586, -0.51319373778697031, -1.672414782127583}, + }; + std::mt19937_64 random(0x3389); + ArrayXd x(size), y(size); + for (Index i = 0; i < size; ++i) { + const double fraction = double(random() >> 11) / 9007199254740992.0; + x(i) = lower[domain] + (upper[domain] - lower[domain]) * fraction; + if (i % 2) x(i) = -x(i); + } + x(0) = reference[domain][0]; + auto evaluate = [&]() { + EIGEN_IF_CONSTEXPR (Function == 0) y = x.sin(); + EIGEN_IF_CONSTEXPR (Function == 1) y = x.cos(); + EIGEN_IF_CONSTEXPR (Function == 2) y = x.tan(); + }; + evaluate(); + const double expected = reference[domain][Function + 1]; + const double bound = 4 * NumTraits<double>::epsilon() * (std::max)(1.0, std::abs(expected)); + if (!(std::abs(y(0) - expected) <= bound)) { + state.SkipWithError("trig result disagrees with MPFR reference"); + return; + } + for (auto _ : state) { + benchmark::DoNotOptimize(x.data()); + evaluate(); + benchmark::DoNotOptimize(y.data()); + benchmark::ClobberMemory(); + } + state.SetItemsProcessed(state.iterations() * size); +} + +BENCHMARK_TEMPLATE(BM_TrigReduction, 0) + ->Name("SinDoubleReduction") + ->ArgsProduct({{256, 4096, 65536}, {0, 1, 2, 3}}) + ->ArgNames({"size", "domain"}); +BENCHMARK_TEMPLATE(BM_TrigReduction, 1) + ->Name("CosDoubleReduction") + ->ArgsProduct({{256, 4096, 65536}, {0, 1, 2, 3}}) + ->ArgNames({"size", "domain"}); +BENCHMARK_TEMPLATE(BM_TrigReduction, 2) + ->Name("TanDoubleReduction") + ->ArgsProduct({{256, 4096, 65536}, {0, 1, 2, 3}}) + ->ArgNames({"size", "domain"});
diff --git a/cmake/EigenSmokeTestList.cmake b/cmake/EigenSmokeTestList.cmake index d21c252..0d65626 100644 --- a/cmake/EigenSmokeTestList.cmake +++ b/cmake/EigenSmokeTestList.cmake
@@ -152,6 +152,7 @@ triangular_1 tridiagonal_eigensolver_1 tridiagonal_eigensolver_3 + trig_reduction_fastmath type_alias lu_9 umeyama_3
diff --git a/contrib/Eigen/src/NumericalDiff/NumericalDiff.h b/contrib/Eigen/src/NumericalDiff/NumericalDiff.h index 2d0ffbe..f2aba41 100644 --- a/contrib/Eigen/src/NumericalDiff/NumericalDiff.h +++ b/contrib/Eigen/src/NumericalDiff/NumericalDiff.h
@@ -19,21 +19,6 @@ namespace Eigen { -namespace internal { - -// Keeps -ffast-math from folding the rounded evaluation points back into x and h, e.g. (x + h) - x -// into h. EIGEN_OPTIMIZATION_BARRIER is an asm operand constraint restricted to plain types (see -// Macros.h), so class-type scalars go unguarded. -template <typename Scalar, std::enable_if_t<std::is_floating_point<Scalar>::value, int> = 0> -EIGEN_STRONG_INLINE void numerical_diff_barrier(Scalar& x) { - EIGEN_UNUSED_VARIABLE(x); - EIGEN_OPTIMIZATION_BARRIER(x) -} -template <typename Scalar, std::enable_if_t<!std::is_floating_point<Scalar>::value, int> = 0> -EIGEN_STRONG_INLINE void numerical_diff_barrier(Scalar&) {} - -} // namespace internal - enum NumericalDiffMode { Forward, Central }; /** @@ -149,7 +134,8 @@ // of x[j] + h perturbs h by up to ulp(x[j]) <= epsilon |x[j]| <= (epsilon/rel_tol) * h <= h, comparable // to the error of the difference quotient itself. Scalar x_plus = _x[j] + h; - internal::numerical_diff_barrier(x_plus); + // Keeps -ffast-math from folding (x[j] + h) - x[j] back into h. + internal::pbarrier(x_plus); h = x_plus - _x[j]; switch (mode) { case Forward: @@ -166,7 +152,7 @@ // x[j] - h can round (a tie when x[j] < 0 has |x[j]| within h above a power of two), so // divide by the separation of the two evaluation points rather than by 2*h. Scalar x_minus = _x[j] - h; - internal::numerical_diff_barrier(x_minus); + internal::pbarrier(x_minus); x[j] = x_minus; Functor::operator()(x, val1); nfev++;
diff --git a/contrib/Eigen/src/NumericalDiff/Ridders.h b/contrib/Eigen/src/NumericalDiff/Ridders.h index 0c9042e..4b2109a 100644 --- a/contrib/Eigen/src/NumericalDiff/Ridders.h +++ b/contrib/Eigen/src/NumericalDiff/Ridders.h
@@ -182,8 +182,8 @@ for (int step = 0; step < m_parameters.maxSteps; ++step) { Scalar xPlus = input[j] + h; Scalar xMinus = input[j] - h; - internal::numerical_diff_barrier(xPlus); - internal::numerical_diff_barrier(xMinus); + internal::pbarrier(xPlus); + internal::pbarrier(xMinus); if (!(numext::isfinite)(xPlus) || !(numext::isfinite)(xMinus)) break; if (step > 0 && xPlus == previousPlus && xMinus == previousMinus) break;
diff --git a/contrib/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h b/contrib/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h index 331211f..a3e91a1 100644 --- a/contrib/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h +++ b/contrib/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h
@@ -192,7 +192,7 @@ #if defined(__FAST_MATH__) || defined(__ASSOCIATIVE_MATH__) || EIGEN_COMP_NVHPC // Unfortunately fast-math reassociates (exp(|x|/2) * exp(|x|/2)) * i0e(x) which reintroduces the overflow, // so we need to tell it not to. - EIGEN_OPTIMIZATION_BARRIER(scaled) + pbarrier(scaled); #endif // i0e(inf) = 0 would turn exp(inf) * i0e(inf) into NaN; i0(+-inf) = +inf. const T inf = pinf<T>(); @@ -353,7 +353,7 @@ T scaled = pmul(half_exp, i1e); #if defined(__FAST_MATH__) || defined(__ASSOCIATIVE_MATH__) || EIGEN_COMP_NVHPC // Reassociating back to (exp(|x|/2) * exp(|x|/2)) * i1e(x) would reintroduce the overflow. - EIGEN_OPTIMIZATION_BARRIER(scaled) + pbarrier(scaled); #endif // i1(+-inf) = +-inf; the scaled product would be inf * 0 = NaN. return pselect(pcmp_eq(ax, pinf<T>()), x, pmul(half_exp, scaled));
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 862a9f8..1f8c4b7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt
@@ -616,6 +616,8 @@ ei_add_test(householder_fastmath "${EIGEN_FASTMATH_FLAGS}") ei_add_test(trsm_packet_cancellation_fastmath "${EIGEN_FASTMATH_FLAGS}") ei_add_test(trig_sign_fastmath "${EIGEN_FASTMATH_FLAGS}") +ei_add_test(trig_reduction_fastmath "${EIGEN_FASTMATH_FLAGS}") +ei_add_test(optimization_barrier "${EIGEN_FASTMATH_FLAGS}") if(QT4_FOUND) ei_add_test(qtvector "" "${QT_QTCORE_LIBRARY}")
diff --git a/test/optimization_barrier.cpp b/test/optimization_barrier.cpp new file mode 100644 index 0000000..0bbc04e --- /dev/null +++ b/test/optimization_barrier.cpp
@@ -0,0 +1,139 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include "main.h" +#include <Eigen/Core> +#ifdef __ARM_FEATURE_SVE +#include <arm_sve.h> +#endif + +enum class BarrierEnum : unsigned long long { Value = 0xfedcba9876543210ull }; +struct BarrierRecord { + unsigned long long value; + bool operator==(const BarrierRecord& other) const { return value == other.value; } +}; +union BarrierUnion { + unsigned long long value; + double real; + bool operator==(const BarrierUnion& other) const { return value == other.value; } +}; +// GCC rejects register and general asm constraints for this struct on every target. +struct BarrierTriple { + double x, y, z; + bool operator==(const BarrierTriple& other) const { return x == other.x && y == other.y && z == other.z; } +}; + +template <typename T> +void check_barrier_value(T value) { + const T expected = value; + EIGEN_OPTIMIZATION_BARRIER(value); + VERIFY(value == expected); + internal::pbarrier(value); + VERIFY(value == expected); +} + +template <typename Scalar> +EIGEN_DONT_INLINE Scalar rounded_by_barrier(Scalar x) { + // A constant bias: MSVC /fp:fast folds (x + bias) - bias only when it can see the constant. + constexpr Scalar bias = Scalar(1ull << (NumTraits<Scalar>::digits() - 1)); + Scalar rounded = x + bias; + EIGEN_OPTIMIZATION_BARRIER(rounded); + return rounded - bias; +} + +EIGEN_DECLARE_TEST(optimization_barrier) { + check_barrier_value(false); + check_barrier_value(static_cast<signed char>(-17)); + check_barrier_value(static_cast<unsigned short>(65521)); + check_barrier_value(-1234567); + check_barrier_value(0xfedcba9876543210ull); + check_barrier_value(BarrierEnum::Value); + check_barrier_value(nullptr); + check_barrier_value(1.25f); + check_barrier_value(1.25); + check_barrier_value(1.25L); + int pointee = 7; + check_barrier_value(&pointee); + volatile decltype(nullptr) null_value = nullptr; + EIGEN_OPTIMIZATION_BARRIER(null_value); + VERIFY(null_value == nullptr); + volatile double real_value = 1.25; + EIGEN_OPTIMIZATION_BARRIER(real_value); + VERIFY(real_value == 1.25); + + check_barrier_value(BarrierRecord{0xfedcba9876543210ull}); + check_barrier_value(BarrierUnion{0xfedcba9876543210ull}); + check_barrier_value(BarrierTriple{1.25, -2.5, 3.75}); + check_barrier_value(std::complex<double>(1.25, -2.5)); + int array[] = {1, -2, 3, -4}; + EIGEN_OPTIMIZATION_BARRIER(array); + for (int i = 0; i < 4; ++i) VERIFY_IS_EQUAL(array[i], (i % 2 ? -1 : 1) * (i + 1)); + double wide_array[8] = {1, -2, 3, -4, 5, -6, 7, -8}; + EIGEN_OPTIMIZATION_BARRIER(wide_array); + for (int i = 0; i < 8; ++i) VERIFY_IS_EQUAL(wide_array[i], (i % 2 ? -1 : 1) * (i + 1)); + internal::eigen_packet_wrapper<int, 7> wrapped(-1234567); + internal::pbarrier(wrapped); + VERIFY_IS_EQUAL(static_cast<int&>(wrapped), -1234567); +#ifdef __ARM_FEATURE_SVE + // Sizeless: sizeof(svfloat64_t) is ill-formed, and there is no operator== for check_barrier_value. + svfloat64_t sizeless = svdup_n_f64(1.25); + EIGEN_OPTIMIZATION_BARRIER(sizeless); + VERIFY(!svptest_any(svptrue_b64(), svcmpne_n_f64(svptrue_b64(), sizeless, 1.25))); + internal::pbarrier(sizeless); + VERIFY(!svptest_any(svptrue_b64(), svcmpne_n_f64(svptrue_b64(), sizeless, 1.25))); +#endif + +#if !defined(EIGEN_GPU_COMPILE_PHASE) && !defined(SYCL_DEVICE_ONLY) && \ + (EIGEN_COMP_GNUC || (defined(__clang__) && defined(_MSC_VER)) || (EIGEN_COMP_MSVC && defined(_M_FP_FAST))) + VERIFY_IS_EQUAL(rounded_by_barrier(1.25f), 1.0f); + VERIFY_IS_EQUAL(rounded_by_barrier(1.25), 1.0); +#endif + +#if !defined(EIGEN_GPU_COMPILE_PHASE) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_CLANG || EIGEN_COMP_MINGW) && \ + (EIGEN_ARCH_PPC || EIGEN_ARCH_ARM_OR_ARM64 || EIGEN_ARCH_i386_OR_x86_64) + using internal::optimization_barrier_in_memory; + STATIC_CHECK(optimization_barrier_in_memory<BarrierRecord>::value); + STATIC_CHECK(optimization_barrier_in_memory<BarrierUnion>::value); + STATIC_CHECK(optimization_barrier_in_memory<BarrierTriple>::value); + STATIC_CHECK(optimization_barrier_in_memory<std::complex<double>>::value); + STATIC_CHECK(optimization_barrier_in_memory<Eigen::half>::value); + STATIC_CHECK(optimization_barrier_in_memory<int[4]>::value); + STATIC_CHECK(optimization_barrier_in_memory<double[8]>::value); + STATIC_CHECK(!optimization_barrier_in_memory<double>::value); + STATIC_CHECK(!optimization_barrier_in_memory<long double>::value); + STATIC_CHECK(!optimization_barrier_in_memory<int*>::value); + STATIC_CHECK(!optimization_barrier_in_memory<decltype(nullptr)>::value); + STATIC_CHECK(!optimization_barrier_in_memory<BarrierEnum>::value); +#endif + +#if !defined(EIGEN_GPU_COMPILE_PHASE) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_CLANG || EIGEN_COMP_MINGW) && \ + (EIGEN_ARCH_ARM_OR_ARM64 || EIGEN_ARCH_i386_OR_x86_64) + using internal::optimization_barrier_in_register; +#if (EIGEN_ARCH_ARM_OR_ARM64 && defined(__ARM_FP) && (__ARM_FP & 4)) || (EIGEN_ARCH_i386_OR_x86_64 && defined(__SSE2__)) + STATIC_CHECK(optimization_barrier_in_register<float>::value); +#endif +#if (EIGEN_ARCH_ARM_OR_ARM64 && defined(__ARM_FP) && (__ARM_FP & 8)) || (EIGEN_ARCH_i386_OR_x86_64 && defined(__SSE2__)) + STATIC_CHECK(optimization_barrier_in_register<double>::value); +#endif + STATIC_CHECK(!optimization_barrier_in_register<long double>::value); + STATIC_CHECK(!optimization_barrier_in_register<int>::value); + STATIC_CHECK(!optimization_barrier_in_register<int*>::value); + STATIC_CHECK(!optimization_barrier_in_register<decltype(nullptr)>::value); + STATIC_CHECK(!optimization_barrier_in_register<BarrierEnum>::value); + STATIC_CHECK(!optimization_barrier_in_register<BarrierRecord>::value); + STATIC_CHECK(!optimization_barrier_in_register<BarrierUnion>::value); + STATIC_CHECK(!optimization_barrier_in_register<Eigen::half>::value); + STATIC_CHECK(!optimization_barrier_in_register<Eigen::bfloat16>::value); + STATIC_CHECK(!optimization_barrier_in_register<std::complex<double>>::value); + STATIC_CHECK(!optimization_barrier_in_register<int[4]>::value); +#ifdef __ARM_FEATURE_SVE + STATIC_CHECK(!optimization_barrier_in_register<svfloat64_t>::value); +#endif +#if (EIGEN_ARCH_i386_OR_x86_64 && defined(__SSE2__)) || defined(__ARM_NEON) + using BarrierVector = float __attribute__((vector_size(16))); + using WideBarrierVector = float __attribute__((vector_size(128))); + STATIC_CHECK(optimization_barrier_in_register<BarrierVector>::value); + STATIC_CHECK(!optimization_barrier_in_register<WideBarrierVector>::value); +#endif +#endif +}
diff --git a/test/packetmath.cpp b/test/packetmath.cpp index d9bfb04..bebd03a 100644 --- a/test/packetmath.cpp +++ b/test/packetmath.cpp
@@ -736,7 +736,7 @@ } }; -// Ensure optimization barrier compiles and doesn't modify contents. +// Ensure the macro and packet optimization barriers preserve every lane. // Only applies to raw types, so will not work for std::complex, Eigen::half // or Eigen::bfloat16. For those you would need to refer to an underlying // storage element. @@ -750,11 +750,18 @@ Packet, std::enable_if_t<!NumTraits<Packet>::IsComplex && !std::is_same<Packet, Eigen::half>::value && !std::is_same<Packet, Eigen::bfloat16>::value>> { static void run() { - typedef typename internal::unpacket_traits<Packet>::type Scalar; - Scalar s = internal::random<Scalar>(); - Packet barrier = internal::pset1<Packet>(s); + using Scalar = typename internal::unpacket_traits<Packet>::type; + constexpr int PacketSize = internal::unpacket_traits<Packet>::size; + Scalar input[PacketSize], output[PacketSize]; + for (int i = 0; i < PacketSize; ++i) input[i] = internal::random<Scalar>(); + Packet barrier = internal::ploadu<Packet>(input); EIGEN_OPTIMIZATION_BARRIER(barrier); - eigen_assert(s == internal::pfirst(barrier) && "EIGEN_OPTIMIZATION_BARRIER"); + internal::pstoreu(output, barrier); + for (int i = 0; i < PacketSize; ++i) VERIFY_IS_EQUAL(input[i], output[i]); + // Deduced, as in the kernels, so wrapped packets take the unwrapping overload. + internal::pbarrier(barrier); + internal::pstoreu(output, barrier); + for (int i = 0; i < PacketSize; ++i) VERIFY_IS_EQUAL(input[i], output[i]); } };
diff --git a/test/trig_reduction_fastmath.cpp b/test/trig_reduction_fastmath.cpp new file mode 100644 index 0000000..036b416 --- /dev/null +++ b/test/trig_reduction_fastmath.cpp
@@ -0,0 +1,89 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#include "packetmath_test_shared.h" + +void check_trig_reduction() { + // MPFR 4.2.1, 256-bit precision, rounded to binary64. Include both sides of the medium/large reduction cutoffs. + const double cases[][4] = { + {1.5707963267948963, 1, 2.8327694488239898e-16, 3530114321217157.5}, + {1.5707963267948966, 1, 6.123233995736766e-17, 16331239353195370}, + {1.5707963267948968, 1, -1.6081226496766366e-16, -6218431163823738}, + {3.1415926535897927, 5.6655388976479796e-16, -1, -5.6655388976479796e-16}, + {3.1415926535897931, 1.2246467991473532e-16, -1, -1.2246467991473532e-16}, + {3.1415926535897936, -3.2162452993532732e-16, -1, 3.2162452993532732e-16}, + {14.137166941154067, 1, 2.3274478990165594e-15, 429655160239049.94}, + {14.137166941154069, 1, 5.5109105961630896e-16, 1814582150355041}, + {14.137166941154071, 1, -1.2252657797839415e-15, -816149456305174.88}, + {0.29999999999999999, 0.29552020666133955, 0.95533648912560598, 0.30933624960962325}, + {4.7123889803846897, -1, -1.8369701987210297e-16, 5443746451065123}, + {14.999999999999998, 0.65028784015711827, -0.75968791285882009, -0.85599340090852183}, + {15, 0.65028784015711683, -0.75968791285882131, -0.85599340090851872}, + {15.000000000000002, 0.65028784015711549, -0.75968791285882242, -0.85599340090851572}, + {20, 0.91294525072762767, 0.40808206181339196, 2.2371609442247422}, + {100, -0.50636564110975879, 0.86231887228768389, -0.58721391515692911}, + {123.456, -0.80393736857282394, -0.5947139710921574, 1.3518050821917635}, + {1000, 0.82687954053200252, 0.56237907629070294, 1.4703241557027185}, + {12345.678900000001, -0.70344192126325633, 0.71075274421521484, -0.98971397154430851}, + {123456.789, -0.99866408234322457, 0.051672532718701383, -19.326787943215852}, + {1000000, -0.34999350217129294, 0.93675212753314474, -0.37362445398759903}, + {10000000, 0.4205477931907825, -0.90727038618173961, -0.46353082785018906}, + {100000000, 0.93163902710972601, -0.36338508935569053, -2.5637789067283774}, + {1000000000, 0.54584344944869956, 0.83788718136390239, 0.65145220214514132}, + {10000000000, -0.48750602508751067, 0.87311962267685606, -0.55834963781124181}, + {100000000000, 0.92869366049659197, 0.37084779216471114, 2.5042448144982212}, + {1000000000000, -0.61123870237688949, 0.79144630185289022, -0.77230596813187613}, + {6907362400000, 0.5754525274479354, 0.81783518428457402, 0.70362896889955862}, + {10000000000000, -0.28888529481752512, 0.95736371690083999, -0.30175082856983471}, + {99999999999999.984, -0.19410480015836853, -0.9809807982603328, 0.1978680933435121}, + {100000000000000, -0.20940830749645231, -0.97782828796853249, 0.21415652428250395}, + {100000000000000.02, -0.22466069079959752, -0.97443705492425092, 0.23055433869667633}, + {1000000000000000, 0.85827279317023586, -0.51319373778697031, -1.672414782127583}, + // Medium-path inputs within 6e-19 to 9e-7 of k*pi/2, for k from 29 to 6e13. + {45.553093477052002, 1, -6.1898063658835771e-19, -1.6155594228467482e+18}, + {364.42474781641602, 4.9518450927068616e-18, 1, 4.9518450927068616e-18}, + {857879.41830474522, -1, -1.2354480663331819e-14, 80942293508783.953}, + {956274858.90980339, -1, 2.2360885370100136e-11, -44720948363.571968}, + {93293944418765.062, 9.1063643066365511e-07, 0.99999999999958533, 9.1063643066403267e-07}, + }; + const Index count = sizeof(cases) / sizeof(cases[0]); + ArrayXd x(2 * count), sin_ref(2 * count), cos_ref(2 * count), tan_ref(2 * count); + for (Index i = 0; i < count; ++i) { + x[2 * i] = cases[i][0]; + x[2 * i + 1] = -cases[i][0]; + sin_ref[2 * i] = cases[i][1]; + sin_ref[2 * i + 1] = -cases[i][1]; + cos_ref[2 * i] = cos_ref[2 * i + 1] = cases[i][2]; + tan_ref[2 * i] = cases[i][3]; + tan_ref[2 * i + 1] = -cases[i][3]; + } + for (Index offset = 0; offset < 2; ++offset) { + for (Index size = 1; size <= x.size() - offset; ++size) { + const ArrayXd sine = x.segment(offset, size).sin(); + const ArrayXd cosine = x.segment(offset, size).cos(); + const ArrayXd tangent = x.segment(offset, size).tan(); + VERIFY(test::areWithinUlps(sin_ref.data() + offset, sine.data(), static_cast<int>(size), 4)); + VERIFY(test::areWithinUlps(cos_ref.data() + offset, cosine.data(), static_cast<int>(size), 4)); + VERIFY(test::areWithinUlps(tan_ref.data() + offset, tangent.data(), static_cast<int>(size), 4)); + } + } +} + +template <typename Packet> +void check_error_free_transformations() { + const double low = std::ldexp(1.0, -54); + Packet hi, lo; + internal::twosum(internal::pset1<Packet>(low), internal::pset1<Packet>(1.0), hi, lo); + VERIFY_IS_EQUAL(internal::pfirst(hi), 1.0); + VERIFY_IS_EQUAL(internal::pfirst(lo), low); + const Packet x = internal::pset1<Packet>(1.0 + std::ldexp(1.0, -30)); + internal::twoprod<Packet, true>(x, x, hi, lo); + VERIFY_IS_EQUAL(internal::pfirst(hi), 1.0 + std::ldexp(1.0, -29)); + VERIFY_IS_EQUAL(internal::pfirst(lo), std::ldexp(1.0, -60)); +} + +EIGEN_DECLARE_TEST(trig_reduction_fastmath) { + CALL_SUBTEST(check_trig_reduction()); + CALL_SUBTEST(check_error_free_transformations<double>()); + CALL_SUBTEST(check_error_free_transformations<internal::packet_traits<double>::type>()); +}