Core: Search all(), any(), allFinite() and isZero() with a blocked predicate search libeigen/eigen!3324
diff --git a/Eigen/Core b/Eigen/Core index 7851324..4677123 100644 --- a/Eigen/Core +++ b/Eigen/Core
@@ -408,6 +408,7 @@ #include "src/Core/SkewSymmetricMatrix3.h" #include "src/Core/Redux.h" #include "src/Core/Visitor.h" +#include "src/Core/Predicator.h" #include "src/Core/FindCoeff.h" #include "src/Core/ApproxComparisons.h" #include "src/Core/Swap.h"
diff --git a/Eigen/src/Core/Predicator.h b/Eigen/src/Core/Predicator.h new file mode 100644 index 0000000..eaba4d6 --- /dev/null +++ b/Eigen/src/Core/Predicator.h
@@ -0,0 +1,555 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +#ifndef EIGEN_PREDICATOR_H +#define EIGEN_PREDICATOR_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +namespace internal { + +// Predicate reductions all_of and any_of. A predicate tests a scalar with operator() and a packet with +// packetOp, which returns a truth packet: each lane is zero where the predicate is false and, where it is true, the +// lane mask pcmp_* produces (all ones, or one for bool packets). Truth packets combine with pand and por, so the +// traversals below reduce and branch once per block of packets rather than once per packet. Results do not depend on +// the grouping, but a short-circuit may read past the deciding coefficient, up to the end of its block. +template <typename Scalar> +struct nonzero_predicate { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator()(const Scalar& x) const { return x != Scalar(0); } + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const { + return pandnot(ptrue(x), pcmp_eq(x, pzero(x))); + } +}; +template <typename Scalar> +struct functor_traits<nonzero_predicate<Scalar>> { + static constexpr bool PacketAccess = packet_traits<Scalar>::HasCmp; + // Accepts that ARMv7 NEON reads a subnormal float as zero: packets compare there already, and grouped scalar tests + // may be vectorized onto it. Testing float one at a time instead cost 0.7-0.85x at sizes with a scalar tail on + // Cortex-A72. + static constexpr bool GroupScalars = true; +}; + +template <typename Scalar> +struct isfinite_predicate { + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator()(const Scalar& x) const { return (numext::isfinite)(x); } + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const { + return pisfinite(x); + } +}; +template <typename Scalar> +struct functor_traits<isfinite_predicate<Scalar>> { + static constexpr bool PacketAccess = packet_traits<Scalar>::HasCmp; + static constexpr bool GroupScalars = true; +}; + +// Separate specializations keep custom scalars that support only one of isApprox and isMuchSmallerThan compiling. +template <typename Scalar, bool Approximate> +struct fuzzy_constant_predicate { + using RealScalar = typename NumTraits<Scalar>::Real; + Scalar value; + RealScalar precision; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator()(const Scalar& x) const { + return internal::isApprox(x, value, precision); + } + // Requires 0 <= precision <= 1 and a finite value. Then min(|x|, |v|) p = min(|x| p, |v| p) exactly, and |x| p is + // x copysign(p, v) where x has the sign of v; where it does not, |x - v| >= max(|x|, |v|) > p min(|x|, |v|) unless + // x = v = 0, and x copysign(p, v) <= 0 rejects it as well. + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const { + const Packet bound = pset1<Packet>(numext::abs(value) * precision); + const Packet scale = pset1<Packet>(numext::copysign(precision, value)); + return pcmp_le(pabs(psub(x, pset1<Packet>(value))), pmin(pmul(x, scale), bound)); + } + EIGEN_DEVICE_FUNC bool packetOpApplies() const { + return precision >= RealScalar(0) && precision <= RealScalar(1) && (numext::isfinite)(value); + } +}; +template <typename Scalar> +struct fuzzy_constant_predicate<Scalar, false> { + Scalar value; + typename NumTraits<Scalar>::Real precision; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator()(const Scalar& x) const { + return internal::isMuchSmallerThan(x, Scalar(1), precision); + } + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const { + return pcmp_le(pabs(x), pset1<Packet>(precision)); + } +}; +template <typename Scalar, bool Approximate> +struct functor_traits<fuzzy_constant_predicate<Scalar, Approximate>> { + // ARMv7 NEON flushes subnormal float operands and results; scalar VFP does not. + // Explicitly flushing would change exact-zero checks; normal operands can also have a subnormal difference. + static constexpr bool PacketAccess = + (std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value) && + !(EIGEN_ARCH_ARM && std::is_same<Scalar, float>::value) && packet_traits<Scalar>::HasAbs && + packet_traits<Scalar>::HasCmp && + (!Approximate || + (packet_traits<Scalar>::HasSub && packet_traits<Scalar>::HasMin && packet_traits<Scalar>::HasMul)); + // Grouped scalar compares may be vectorized onto the unit PacketAccess excludes. + static constexpr bool GroupScalars = PacketAccess; + // The approximate test saturates the vector units, so larger blocks barely speed up full scans but read further + // past a deciding coefficient: 2 packets per exit is never slower than master's exit per packet. + static constexpr int BlockPackets = Approximate ? 2 : 8; +}; + +// A predicate's scalar test alone, tested one coefficient at a time. +template <typename Predicate> +struct scalar_predicate { + Predicate pred; + template <typename Scalar> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator()(const Scalar& x) const { + return pred(x); + } +}; +template <typename Predicate> +struct functor_traits<scalar_predicate<Predicate>> { + static constexpr bool PacketAccess = false; +}; + +// Searches for a coefficient on which the predicate is Target: any_of searches for true, all_of for false. Packets +// combine truth masks and test the combined mask; scalars search for a true test_scalar by OR, which vectorizes +// better than an AND of test results. +struct packet_predicate { + template <typename Predicate, typename Packet> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet apply(const Predicate& pred, const Packet& x) { + return pred.packetOp(x); + } +}; +template <bool Target> +struct predicate_search; +template <> +struct predicate_search<true> : packet_predicate { + template <typename Packet> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet combine(const Packet& a, const Packet& b) { + return por(a, b); + } + template <typename Packet> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool found(const Packet& m) { + return predux_any(m); + } + template <typename Predicate, typename Scalar> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool test_scalar(const Predicate& pred, const Scalar& x) { + return pred(x); + } +}; +template <> +struct predicate_search<false> : packet_predicate { + template <typename Packet> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Packet combine(const Packet& a, const Packet& b) { + return pand(a, b); + } + // pandnot rather than pnot, which some targets implement only bytewise. + template <typename Packet> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool found(const Packet& m) { + return predux_any(pandnot(ptrue(m), m)); + } + template <typename Predicate, typename Scalar> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool test_scalar(const Predicate& pred, const Scalar& x) { + return !pred(x); + } +}; + +// all_of(nonzero) searches for zero lanes, x == 0 being the exact complement of x != 0: a single compare per packet +// and no complement of the combined mask, which some compilers otherwise lower through narrowing and widening. +// Scalars still test the predicate: a coefficient is found where it is false. +struct zero_lane_search : predicate_search<true> { + template <typename Predicate, typename Packet> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet apply(const Predicate&, const Packet& x) { + return pcmp_eq(x, pzero(x)); + } + template <typename Predicate, typename Scalar> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool test_scalar(const Predicate& pred, const Scalar& x) { + return !pred(x); + } +}; +// any_of on bool lanes reduces them directly, treating every nonzero byte as true as the scalar reduction does. +struct bool_lane_search : predicate_search<true> { + template <typename Predicate, typename Packet> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet apply(const Predicate&, const Packet& x) { + return x; + } +}; +template <bool Target, typename Predicate> +struct predicate_search_for { + using type = predicate_search<Target>; +}; +template <typename Scalar> +struct predicate_search_for<false, nonzero_predicate<Scalar>> { + using type = zero_lane_search; +}; +template <> +struct predicate_search_for<true, nonzero_predicate<bool>> { + using type = bool_lane_search; +}; + +// Coefficient access along one inner vector, or along the whole expression with linear access. +template <typename Evaluator, bool Linear = Evaluator::LinearAccess> +struct predicate_segment { + using Scalar = typename Evaluator::Scalar; + const Evaluator& eval; + Index length; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE predicate_segment(const Evaluator& evaluator, Index, Index size) + : eval(evaluator), length(size) {} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index size() const { return length; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename Evaluator::CoeffReturnType coeff(Index i) const { + return eval.coeff(i); + } + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packet(Index i) const { + return eval.template packet<Packet>(i); + } +}; +template <typename Evaluator> +struct predicate_segment<Evaluator, false> { + using Scalar = typename Evaluator::Scalar; + const Evaluator& eval; + Index outer; + Index length; + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE predicate_segment(const Evaluator& evaluator, Index j, Index size) + : eval(evaluator), outer(j), length(size) {} + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index size() const { return length; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename Evaluator::CoeffReturnType coeff(Index i) const { + return Evaluator::IsRowMajor ? eval.coeff(outer, i) : eval.coeff(i, outer); + } + template <typename Packet> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packet(Index i) const { + return Evaluator::IsRowMajor ? eval.template packet<Packet>(outer, i) : eval.template packet<Packet>(i, outer); + } +}; + +// Whether a predicate's scalar tests may be grouped; predicates that do not say are tested one at a time. +template <typename Predicate, typename = void> +struct predicate_groups_scalars : std::false_type {}; +template <typename Predicate> +struct predicate_groups_scalars<Predicate, void_t<decltype(functor_traits<Predicate>::GroupScalars)>> + : bool_constant<static_cast<bool>(functor_traits<Predicate>::GroupScalars)> {}; + +// Packets tested between exits on the long path; predicates that do not say test blocks of 8. +template <typename Predicate, typename = void> +struct predicate_block_packets : std::integral_constant<int, 8> {}; +template <typename Predicate> +struct predicate_block_packets<Predicate, void_t<decltype(functor_traits<Predicate>::BlockPackets)>> + : std::integral_constant<int, functor_traits<Predicate>::BlockPackets> {}; + +// Scalar searches OR their tests without a branch per test; spelling a | b on the calls directly draws clang's +// -Wbitwise-instead-of-logical. +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool predicate_or(bool a, bool b) { return a | b; } + +// At most three coefficients from i on, as straight-line tests. +template <typename Search, typename Segment, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool predicate_search_upto3(const Segment& seg, Index i, Predicate pred) { + const Index size = seg.size(); + if (i == size) return false; + bool b = Search::test_scalar(pred, seg.coeff(i)); + if (i + 1 < size) b = predicate_or(b, Search::test_scalar(pred, seg.coeff(i + 1))); + if (i + 2 < size) b = predicate_or(b, Search::test_scalar(pred, seg.coeff(i + 2))); + return b; +} + +template <typename Search, typename Segment, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool predicate_search_scalar(const Segment& seg, Index i, Predicate pred) { + const Index size = seg.size(); + for (; i + 3 < size; i += 4) { + bool b = predicate_or( + predicate_or(Search::test_scalar(pred, seg.coeff(i)), Search::test_scalar(pred, seg.coeff(i + 1))), + predicate_or(Search::test_scalar(pred, seg.coeff(i + 2)), Search::test_scalar(pred, seg.coeff(i + 3)))); + if EIGEN_PREDICT_FALSE (b) return true; + } + return predicate_search_upto3<Search>(seg, i, pred); +} + +// One coefficient at a time with an exit after each, a shape compilers do not vectorize. +template <typename Search, typename Segment, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool predicate_search_sequential(const Segment& seg, Index i, Predicate pred) { + for (const Index size = seg.size(); i < size; ++i) + if (Search::test_scalar(pred, seg.coeff(i))) return true; + return false; +} + +// Scalar tests from i on: grouped only where the predicate allows it, since grouped scalar tests may be vectorized +// onto the unit a predicate excludes (ARMv7 NEON, which flushes subnormals, turns grouped float compares into vceq). +template <typename Search, typename Segment, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool predicate_search_tail(const Segment& seg, Index i, Predicate pred) { + if (predicate_groups_scalars<Predicate>::value) return predicate_search_scalar<Search>(seg, i, pred); + return predicate_search_sequential<Search>(seg, i, pred); +} + +// Scalar blocks OR their tests into an unsigned integer as wide as the scalar. gcc vectorizes that reduction but not an +// OR of bools, and a wider accumulator widens every test. +template <typename Scalar, typename Unsigned = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type> +struct predicate_block_accumulator { + using type = Unsigned; +}; +template <typename Scalar> +struct predicate_block_accumulator<Scalar, void> { + using type = bool; +}; + +template <typename Search, typename Packet, bool Vectorize> +struct predicate_search_segment { + static constexpr Index ShortSize = 0; + // Without packets, blocks of 64 tests without an exit, which compilers vectorize, and one exit check per block; + // starting each block from the identity keeps its trip count regular once inlined. Only where the predicate allows + // grouped tests, as for predicate_search_tail. + template <typename Segment, typename Predicate> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool run(const Segment& seg, Predicate pred) { + if (!predicate_groups_scalars<Predicate>::value) return predicate_search_sequential<Search>(seg, 0, pred); + const Index size = seg.size(); + // The shortest segments first, so they skip the block and group guards. + if (EIGEN_PREDICT_FALSE(size < 4)) return predicate_search_upto3<Search>(seg, 0, pred); + Index i = 0; + using Acc = typename predicate_block_accumulator<typename Segment::Scalar>::type; + for (; i + 64 <= size; i += 64) { + Acc b = 0; + for (Index k = 0; k < 64; ++k) b |= Acc(Search::test_scalar(pred, seg.coeff(i + k))); + if EIGEN_PREDICT_FALSE (b != Acc(0)) return true; + } + return predicate_search_scalar<Search>(seg, i, pred); + } + template <typename Segment, typename Predicate> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool run_short(const Segment& seg, Predicate pred) { + return run(seg, pred); + } +}; +template <typename Search, typename Packet> +struct predicate_search_segment<Search, Packet, true> { + // One helper per packet, taking its operands by reference and returning the mask by value: if the early inliner + // leaves it as a call, the caller needs no addressable packet temporaries. + template <typename Segment, typename Predicate> + static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet test(const Segment& seg, Index i, const Predicate& pred) { + return Search::apply(pred, seg.template packet<Packet>(i)); + } + static constexpr Index PacketSize = unpacket_traits<Packet>::size; + // Segments shorter than ShortSize take run_short. + static constexpr Index ShortSize = 4 * PacketSize; + + // Up to three packets (plus a ragged tail), the common short call: test them directly. + template <typename Segment, typename Predicate> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool run_short(const Segment& seg, Predicate pred) { + const Index size = seg.size(); + Index i = 0; + if (PacketSize <= size) { + if (size < 2 * PacketSize) { + if (Search::found(test(seg, 0, pred))) return true; + i = PacketSize; + } else if (size < 3 * PacketSize) { + if (Search::found(Search::combine(test(seg, 0, pred), test(seg, PacketSize, pred)))) return true; + i = 2 * PacketSize; + } else { + if (Search::found(Search::combine(Search::combine(test(seg, 0, pred), test(seg, PacketSize, pred)), + test(seg, 2 * PacketSize, pred)))) + return true; + i = 3 * PacketSize; + } + if (EIGEN_PREDICT_TRUE(i == size)) return false; + } else if (size < 4 && predicate_groups_scalars<Predicate>::value) { + // Shorter than one packet and than a group of four: only these sizes reach this check. + return predicate_search_upto3<Search>(seg, 0, pred); + } + return predicate_search_tail<Search>(seg, i, pred); + } + + template <int Block, typename Segment, typename Predicate> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool run_blocks(const Segment& seg, Predicate pred) { + const Index size = seg.size(); + Index i = 0; + for (; i + Block * PacketSize <= size; i += Block * PacketSize) { + Packet m = test(seg, i, pred); + for (int k = 1; k < Block; ++k) m = Search::combine(m, test(seg, i + k * PacketSize, pred)); + if EIGEN_PREDICT_FALSE (Search::found(m)) return true; + } + if (Block > 1 && i + PacketSize <= size) { + Packet m = test(seg, i, pred); + for (i += PacketSize; i + PacketSize <= size; i += PacketSize) m = Search::combine(m, test(seg, i, pred)); + if (Search::found(m)) return true; + } + if (EIGEN_PREDICT_TRUE(i == size)) return false; + return predicate_search_tail<Search>(seg, i, pred); + } + + template <typename Segment, typename Predicate> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool run(const Segment& seg, Predicate pred) { + const Index size = seg.size(); + if (size < ShortSize) return run_short(seg, pred); + EIGEN_IF_CONSTEXPR (predicate_block_packets<Predicate>::value < 8) { + return run_blocks<predicate_block_packets<Predicate>::value>(seg, pred); + } + Index i = 0; + // Longer segments: the block loop is laid out off the fall-through path. + if (EIGEN_PREDICT_FALSE(8 * PacketSize <= size)) { + do { + Packet m = Search::combine( + Search::combine(Search::combine(test(seg, i + 0 * PacketSize, pred), test(seg, i + 1 * PacketSize, pred)), + Search::combine(test(seg, i + 2 * PacketSize, pred), test(seg, i + 3 * PacketSize, pred))), + Search::combine(Search::combine(test(seg, i + 4 * PacketSize, pred), test(seg, i + 5 * PacketSize, pred)), + Search::combine(test(seg, i + 6 * PacketSize, pred), test(seg, i + 7 * PacketSize, pred)))); + if EIGEN_PREDICT_FALSE (Search::found(m)) return true; + i += 8 * PacketSize; + } while (i + 8 * PacketSize <= size); + } + if (i + 4 * PacketSize <= size) { + Packet m = + Search::combine(Search::combine(test(seg, i + 0 * PacketSize, pred), test(seg, i + 1 * PacketSize, pred)), + Search::combine(test(seg, i + 2 * PacketSize, pred), test(seg, i + 3 * PacketSize, pred))); + if (Search::found(m)) return true; + i += 4 * PacketSize; + } + if (i + PacketSize <= size) { + Packet m = test(seg, i, pred); + for (i += PacketSize; i + PacketSize <= size; i += PacketSize) m = Search::combine(m, test(seg, i, pred)); + if (Search::found(m)) return true; + } + if (EIGEN_PREDICT_TRUE(i == size)) return false; + return predicate_search_tail<Search>(seg, i, pred); + } +}; + +template <typename Derived, typename Predicate> +struct predicate_reduction { + using Evaluator = visitor_evaluator<Derived>; + using Scalar = typename Derived::Scalar; + static constexpr bool LinearAccess = Evaluator::LinearAccess; + // A segment's compile-time length, so a short fixed-size segment uses a packet that fits it, as assignment does. + static constexpr int SegmentSizeAtCompileTime = + LinearAccess ? int(Derived::SizeAtCompileTime) + : int(Evaluator::IsRowMajor ? Derived::ColsAtCompileTime : Derived::RowsAtCompileTime); + using Packet = typename find_largest_packet<Scalar, SegmentSizeAtCompileTime>::type; + static constexpr bool Vectorize = Evaluator::PacketAccess && packet_traits<Scalar>::Vectorizable && + (unpacket_traits<Packet>::size > 1) && + static_cast<bool>(functor_traits<Predicate>::PacketAccess); + using Segment = predicate_segment<Evaluator>; + + template <bool Target> + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool search(const Derived& xpr, Predicate pred) { + using impl = predicate_search_segment<typename predicate_search_for<Target, Predicate>::type, Packet, Vectorize>; + Evaluator eval(xpr); + EIGEN_IF_CONSTEXPR (LinearAccess) { + return impl::run(Segment(eval, 0, eval.size()), pred); + } else { + // The inner size is read once: through the evaluator's reference to the expression, compilers reload it for + // every inner vector. + const Index innerSize = Evaluator::IsRowMajor ? eval.cols() : eval.rows(); + const Index outerSize = Evaluator::IsRowMajor ? eval.rows() : eval.cols(); + // Choosing the segment path once keeps the length dispatch out of the loop over short inner vectors. + if (innerSize < impl::ShortSize) { + for (Index j = 0; j < outerSize; ++j) + if (impl::run_short(Segment(eval, j, innerSize), pred)) return true; + } else { + for (Index j = 0; j < outerSize; ++j) + if (impl::run(Segment(eval, j, innerSize), pred)) return true; + } + return false; + } + } +}; + +/** \internal \returns true if \a pred holds for every coefficient of \a xpr. */ +template <typename Derived, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool all_of(const Derived& xpr, const Predicate& pred) { + return !predicate_reduction<Derived, Predicate>::template search<false>(xpr, pred); +} + +/** \internal \returns true if \a pred holds for some coefficient of \a xpr. */ +template <typename Derived, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool any_of(const Derived& xpr, const Predicate& pred) { + return predicate_reduction<Derived, Predicate>::template search<true>(xpr, pred); +} + +// The packet test of fuzzy_constant_predicate holds only where packetOpApplies; elsewhere it tests scalars. +template <typename Derived, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool fuzzy_constant_all_of(const Derived& xpr, const Predicate& pred, + std::true_type) { + if (EIGEN_PREDICT_FALSE(!pred.packetOpApplies())) return all_of(xpr, scalar_predicate<Predicate>{pred}); + return all_of(xpr, pred); +} +template <typename Derived, typename Predicate> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool fuzzy_constant_all_of(const Derived& xpr, const Predicate& pred, + std::false_type) { + return all_of(xpr, pred); +} + +template <typename Derived, bool AlwaysTrue = NumTraits<typename traits<Derived>::Scalar>::IsInteger> +struct all_finite_impl { + static EIGEN_DEVICE_FUNC inline bool run(const Derived& /*derived*/) { return true; } +}; +#if !defined(__FINITE_MATH_ONLY__) || !(__FINITE_MATH_ONLY__) +template <typename Derived> +struct all_finite_impl<Derived, false> { + static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool run(const Derived& derived) { + return all_of(derived, isfinite_predicate<typename traits<Derived>::Scalar>()); + } +}; +#endif + +} // end namespace internal + +/** \returns true if all coefficients are true + * + * Example: \include MatrixBase_all.cpp + * Output: \verbinclude MatrixBase_all.out + * + * \sa any(), Cwise::operator<() + */ +template <typename Derived> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool DenseBase<Derived>::all() const { + return internal::all_of(derived(), internal::nonzero_predicate<Scalar>()); +} + +/** \returns true if at least one coefficient is true + * + * \sa all() + */ +template <typename Derived> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool DenseBase<Derived>::any() const { + return internal::any_of(derived(), internal::nonzero_predicate<Scalar>()); +} + +template <typename Derived> +EIGEN_DEVICE_FUNC inline bool DenseBase<Derived>::hasNaN() const { + return derived().cwiseTypedNotEqual(derived()).any(); +} + +/** \returns true if \c *this contains only finite numbers, i.e., no NaN and no +/-INF values. + * + * \sa hasNaN() + */ +template <typename Derived> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool DenseBase<Derived>::allFinite() const { + return internal::all_finite_impl<Derived>::run(derived()); +} + +/** \returns true if all coefficients in this matrix are approximately equal to \a val, to within precision \a prec */ +template <typename Derived> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool DenseBase<Derived>::isApproxToConstant(const Scalar& val, + const RealScalar& prec) const { + using Predicate = internal::fuzzy_constant_predicate<Scalar, true>; + return internal::fuzzy_constant_all_of( + derived(), Predicate{val, prec}, + std::integral_constant<bool, internal::functor_traits<Predicate>::PacketAccess>()); +} + +/** \returns true if *this is approximately equal to the zero matrix, + * within the precision given by \a prec. + * + * Example: \include MatrixBase_isZero.cpp + * Output: \verbinclude MatrixBase_isZero.out + * + * \sa class CwiseNullaryOp, Zero() + */ +template <typename Derived> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool DenseBase<Derived>::isZero(const RealScalar& prec) const { + return internal::all_of(derived(), internal::fuzzy_constant_predicate<Scalar, false>{Scalar(0), prec}); +} + +} // end namespace Eigen + +#endif // EIGEN_PREDICATOR_H
diff --git a/Eigen/src/Core/Visitor.h b/Eigen/src/Core/Visitor.h index 10893e5..0fb61c9 100644 --- a/Eigen/src/Core/Visitor.h +++ b/Eigen/src/Core/Visitor.h
@@ -397,167 +397,6 @@ } }; -} // end namespace internal - -/** Applies the visitor \a visitor to the whole coefficients of the matrix or vector. - * - * The template parameter \a Visitor is the type of the visitor and provides the following interface: - * \code - * struct MyVisitor { - * // called for the first coefficient - * void init(const Scalar& value, Index i, Index j); - * // called for all other coefficients - * void operator() (const Scalar& value, Index i, Index j); - * }; - * \endcode - * - * \note compared to one or two \em for \em loops, visitors offer automatic - * unrolling for small fixed size matrix. - * - * \note if the matrix is empty, then the visitor is left unchanged. - * - * \sa minCoeff(Index*,Index*), maxCoeff(Index*,Index*), DenseBase::redux() - */ -template <typename Derived> -template <typename Visitor> -EIGEN_DEVICE_FUNC void DenseBase<Derived>::visit(Visitor& visitor) const { - using impl = internal::visit_impl<Derived, Visitor, /*ShortCircuitEvaluation*/ false>; - impl::run(derived(), visitor); -} - -namespace internal { - -template <typename Scalar, bool Approximate> -struct fuzzy_constant_visitor { - Scalar value; - typename NumTraits<Scalar>::Real precision; - bool result = true; - - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void operator()(const Scalar& x, Index, Index = 0) { - result = result && (Approximate ? internal::isApprox(x, value, precision) - : internal::isMuchSmallerThan(x, Scalar(1), precision)); - } - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void init(const Scalar& x, Index r, Index c = 0) { (*this)(x, r, c); } - EIGEN_DEVICE_FUNC bool done() const { return !result; } - template <typename Packet> - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void packet(const Packet& x, Index, Index = 0) { - const Packet p = pset1<Packet>(precision); - Packet mask; - EIGEN_IF_CONSTEXPR (Approximate) { - const Packet v = pset1<Packet>(value); - mask = pcmp_le(pabs(psub(x, v)), pmul(pmin(pabs(x), pabs(v)), p)); - } else { - mask = pcmp_le(pabs(x), p); - } - // Reduce the comparison mask directly, without comparing its lanes to zero again. - result = result && !predux_any(pandnot(ptrue(mask), mask)); - } - template <typename Packet> - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void initpacket(const Packet& x, Index r, Index c = 0) { - packet(x, r, c); - } -}; - -template <typename Scalar, bool Approximate> -struct functor_traits<fuzzy_constant_visitor<Scalar, Approximate>> { - static constexpr bool AlreadyInitialized = true; - static constexpr bool LinearAccess = true; - static constexpr int Cost = 4 * NumTraits<Scalar>::AddCost + NumTraits<Scalar>::MulCost; - // ARMv7 NEON flushes subnormal float operands and results; scalar VFP does not. - // Explicitly flushing would change exact-zero checks; normal operands can also have a subnormal difference. - static constexpr bool PacketAccess = - (std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value) && - !(EIGEN_ARCH_ARM && std::is_same<Scalar, float>::value) && packet_traits<Scalar>::HasAbs && - packet_traits<Scalar>::HasCmp && - (!Approximate || - (packet_traits<Scalar>::HasSub && packet_traits<Scalar>::HasMin && packet_traits<Scalar>::HasMul)); -}; - -template <typename Scalar> -using use_fuzzy_constant_visitor = - bool_constant<std::is_same<Scalar, float>::value || std::is_same<Scalar, double>::value>; - -template <bool Approximate, typename Derived, - std::enable_if_t<use_fuzzy_constant_visitor<typename Derived::Scalar>::value, int> = 0> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool fuzzy_constant_all(const Derived& matrix, - const typename Derived::Scalar& value, - const typename Derived::RealScalar& precision) { - fuzzy_constant_visitor<typename Derived::Scalar, Approximate> visitor{value, precision}; - visit_impl<Derived, decltype(visitor), true>::run(matrix, visitor); - return visitor.result; -} - -// Separate overloads keep custom scalars from instantiating unused comparisons in C++14. -template <bool Approximate, typename Derived, - std::enable_if_t<!use_fuzzy_constant_visitor<typename Derived::Scalar>::value && Approximate, int> = 0> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool fuzzy_constant_all(const Derived& matrix, - const typename Derived::Scalar& value, - const typename Derived::RealScalar& precision) { - for (Index j = 0; j < matrix.cols(); ++j) - for (Index i = 0; i < matrix.rows(); ++i) - if (!internal::isApprox(matrix.coeff(i, j), value, precision)) return false; - return true; -} - -template <bool Approximate, typename Derived, - std::enable_if_t<!use_fuzzy_constant_visitor<typename Derived::Scalar>::value && !Approximate, int> = 0> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool fuzzy_constant_all(const Derived& matrix, const typename Derived::Scalar&, - const typename Derived::RealScalar& precision) { - using Scalar = typename Derived::Scalar; - for (Index j = 0; j < matrix.cols(); ++j) - for (Index i = 0; i < matrix.rows(); ++i) - if (!internal::isMuchSmallerThan(matrix.coeff(i, j), Scalar(1), precision)) return false; - return true; -} - -template <typename Scalar> -struct all_visitor { - using result_type = bool; - using Packet = typename packet_traits<Scalar>::type; - EIGEN_DEVICE_FUNC inline void init(const Scalar& value, Index, Index) { res = (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline void init(const Scalar& value, Index) { res = (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline bool all_predux(const Packet& p) const { return predux_all(p); } - EIGEN_DEVICE_FUNC inline void initpacket(const Packet& p, Index, Index) { res = all_predux(p); } - EIGEN_DEVICE_FUNC inline void initpacket(const Packet& p, Index) { res = all_predux(p); } - EIGEN_DEVICE_FUNC inline void operator()(const Scalar& value, Index, Index) { res = res && (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline void operator()(const Scalar& value, Index) { res = res && (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline void packet(const Packet& p, Index, Index) { res = res && all_predux(p); } - EIGEN_DEVICE_FUNC inline void packet(const Packet& p, Index) { res = res && all_predux(p); } - EIGEN_DEVICE_FUNC inline bool done() const { return !res; } - bool res = true; -}; -template <typename Scalar> -struct functor_traits<all_visitor<Scalar>> { - enum { Cost = NumTraits<Scalar>::ReadCost, LinearAccess = true, PacketAccess = packet_traits<Scalar>::HasCmp }; -}; - -template <typename Scalar> -struct any_visitor { - using result_type = bool; - using Packet = typename packet_traits<Scalar>::type; - EIGEN_DEVICE_FUNC inline void init(const Scalar& value, Index, Index) { res = (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline void init(const Scalar& value, Index) { res = (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline bool any_predux(const Packet& p) const { - return predux_any(pandnot(ptrue(p), pcmp_eq(p, pzero(p)))); - } - EIGEN_DEVICE_FUNC inline void initpacket(const Packet& p, Index, Index) { res = any_predux(p); } - EIGEN_DEVICE_FUNC inline void initpacket(const Packet& p, Index) { res = any_predux(p); } - EIGEN_DEVICE_FUNC inline void operator()(const Scalar& value, Index, Index) { res = res || (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline void operator()(const Scalar& value, Index) { res = res || (value != Scalar(0)); } - EIGEN_DEVICE_FUNC inline void packet(const Packet& p, Index, Index) { res = res || any_predux(p); } - EIGEN_DEVICE_FUNC inline void packet(const Packet& p, Index) { res = res || any_predux(p); } - EIGEN_DEVICE_FUNC inline bool done() const { return res; } - bool res = false; -}; -template <> -EIGEN_DEVICE_FUNC inline bool any_visitor<bool>::any_predux(const Packet& p) const { - return predux(p); -} -template <typename Scalar> -struct functor_traits<any_visitor<Scalar>> { - enum { Cost = NumTraits<Scalar>::ReadCost, LinearAccess = true, PacketAccess = packet_traits<Scalar>::HasCmp }; -}; - template <typename Scalar> struct count_visitor { using result_type = Index; @@ -587,81 +426,32 @@ }; }; -// Reduces pisfinite masks directly: their lanes are all-ones or zero, so the coefficients are all finite exactly when -// no lane of the complement is set. predux_all, which all() on isFiniteTyped() uses, compares each lane against zero -// instead, which a compiler can only fold away when it can prove the operand is a mask. -template <typename Scalar> -struct all_finite_visitor { - using result_type = bool; - using Packet = typename packet_traits<Scalar>::type; - EIGEN_DEVICE_FUNC inline bool finite_predux(const Packet& p) const { return !predux_any(pnot(pisfinite(p))); } - EIGEN_DEVICE_FUNC inline void init(const Scalar& value, Index, Index) { res = (numext::isfinite)(value); } - EIGEN_DEVICE_FUNC inline void init(const Scalar& value, Index) { res = (numext::isfinite)(value); } - EIGEN_DEVICE_FUNC inline void initpacket(const Packet& p, Index, Index) { res = finite_predux(p); } - EIGEN_DEVICE_FUNC inline void initpacket(const Packet& p, Index) { res = finite_predux(p); } - EIGEN_DEVICE_FUNC inline void operator()(const Scalar& value, Index, Index) { - res = res && (numext::isfinite)(value); - } - EIGEN_DEVICE_FUNC inline void operator()(const Scalar& value, Index) { res = res && (numext::isfinite)(value); } - EIGEN_DEVICE_FUNC inline void packet(const Packet& p, Index, Index) { res = res && finite_predux(p); } - EIGEN_DEVICE_FUNC inline void packet(const Packet& p, Index) { res = res && finite_predux(p); } - EIGEN_DEVICE_FUNC inline bool done() const { return !res; } - bool res = true; -}; -template <typename Scalar> -struct functor_traits<all_finite_visitor<Scalar>> { - enum { - Cost = NumTraits<Scalar>::ReadCost + NumTraits<Scalar>::MulCost, - LinearAccess = true, - PacketAccess = packet_traits<Scalar>::HasCmp - }; -}; - -template <typename Derived, bool AlwaysTrue = NumTraits<typename traits<Derived>::Scalar>::IsInteger> -struct all_finite_impl { - static EIGEN_DEVICE_FUNC inline bool run(const Derived& /*derived*/) { return true; } -}; -#if !defined(__FINITE_MATH_ONLY__) || !(__FINITE_MATH_ONLY__) -template <typename Derived> -struct all_finite_impl<Derived, false> { - static EIGEN_DEVICE_FUNC inline bool run(const Derived& derived) { - using Visitor = all_finite_visitor<typename traits<Derived>::Scalar>; - Visitor visitor; - visit_impl<Derived, Visitor, /*ShortCircuitEvaluation*/ true>::run(derived, visitor); - return visitor.res; - } -}; -#endif - } // end namespace internal -/** \returns true if all coefficients are true +/** Applies the visitor \a visitor to the whole coefficients of the matrix or vector. * - * Example: \include MatrixBase_all.cpp - * Output: \verbinclude MatrixBase_all.out + * The template parameter \a Visitor is the type of the visitor and provides the following interface: + * \code + * struct MyVisitor { + * // called for the first coefficient + * void init(const Scalar& value, Index i, Index j); + * // called for all other coefficients + * void operator() (const Scalar& value, Index i, Index j); + * }; + * \endcode * - * \sa any(), Cwise::operator<() + * \note compared to one or two \em for \em loops, visitors offer automatic + * unrolling for small fixed size matrix. + * + * \note if the matrix is empty, then the visitor is left unchanged. + * + * \sa minCoeff(Index*,Index*), maxCoeff(Index*,Index*), DenseBase::redux() */ template <typename Derived> -EIGEN_DEVICE_FUNC inline bool DenseBase<Derived>::all() const { - using Visitor = internal::all_visitor<Scalar>; - using impl = internal::visit_impl<Derived, Visitor, /*ShortCircuitEvaluation*/ true>; - Visitor visitor; +template <typename Visitor> +EIGEN_DEVICE_FUNC void DenseBase<Derived>::visit(Visitor& visitor) const { + using impl = internal::visit_impl<Derived, Visitor, /*ShortCircuitEvaluation*/ false>; impl::run(derived(), visitor); - return visitor.res; -} - -/** \returns true if at least one coefficient is true - * - * \sa all() - */ -template <typename Derived> -EIGEN_DEVICE_FUNC inline bool DenseBase<Derived>::any() const { - using Visitor = internal::any_visitor<Scalar>; - using impl = internal::visit_impl<Derived, Visitor, /*ShortCircuitEvaluation*/ true>; - Visitor visitor; - impl::run(derived(), visitor); - return visitor.res; } /** \returns the number of coefficients which evaluate to true @@ -677,42 +467,6 @@ return visitor.res; } -template <typename Derived> -EIGEN_DEVICE_FUNC inline bool DenseBase<Derived>::hasNaN() const { - return derived().cwiseTypedNotEqual(derived()).any(); -} - -/** \returns true if \c *this contains only finite numbers, i.e., no NaN and no +/-INF values. - * - * \sa hasNaN() - */ -template <typename Derived> -EIGEN_DEVICE_FUNC inline bool DenseBase<Derived>::allFinite() const { - return internal::all_finite_impl<Derived>::run(derived()); -} - -/** \returns true if all coefficients in this matrix are approximately equal to \a val, to within precision \a prec */ -template <typename Derived> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool DenseBase<Derived>::isApproxToConstant(const Scalar& val, - const RealScalar& prec) const { - typename internal::nested_eval<Derived, 1>::type self(derived()); - return internal::fuzzy_constant_all<true>(self, val, prec); -} - -/** \returns true if *this is approximately equal to the zero matrix, - * within the precision given by \a prec. - * - * Example: \include MatrixBase_isZero.cpp - * Output: \verbinclude MatrixBase_isZero.out - * - * \sa class CwiseNullaryOp, Zero() - */ -template <typename Derived> -EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool DenseBase<Derived>::isZero(const RealScalar& prec) const { - typename internal::nested_eval<Derived, 1>::type self(derived()); - return internal::fuzzy_constant_all<false>(self, Scalar(0), prec); -} - } // end namespace Eigen #endif // EIGEN_VISITOR_H
diff --git a/benchmarks/Core/CMakeLists.txt b/benchmarks/Core/CMakeLists.txt index 94929cd..e87a4c8 100644 --- a/benchmarks/Core/CMakeLists.txt +++ b/benchmarks/Core/CMakeLists.txt
@@ -35,6 +35,26 @@ eigen_add_benchmark(bench_unary_pow bench_unary_pow.cpp) eigen_add_benchmark(bench_pgather bench_pgather.cpp) eigen_add_benchmark(bench_boolean_reductions bench_boolean_reductions.cpp) +# One binary per family: in a shared TU, inlining and alignment shift with unrelated benchmark bodies. +foreach(family float:All float:Any float:AllFinite float:IsZero float:ApproxConstant + double:All double:AllFinite double:IsZero int:All int:Any + int16_t:All int16_t:Any bool:All bool:Any) + string(REPLACE ":" ";" family ${family}) + list(GET family 0 scalar) + list(GET family 1 operation) + string(TOLOWER ${operation} operation_name) + eigen_add_benchmark(bench_predicate_${operation_name}_${scalar} bench_predicate_reductions.cpp + DEFINITIONS EIGEN_BENCH_PREDICATE_SCALAR=${scalar} EIGEN_BENCH_PREDICATE_OP=${operation}) +endforeach() +foreach(family Array4f:Any Array16f:AllFinite Array33f:IsZero Matrix4d:IsZero Vector3f:All) + string(REPLACE ":" ";" family ${family}) + list(GET family 0 xpr) + list(GET family 1 operation) + string(TOLOWER "${operation}_${xpr}" target_suffix) + eigen_add_benchmark(bench_predicate_${target_suffix} bench_predicate_reductions.cpp + DEFINITIONS EIGEN_BENCH_PREDICATE_XPR=${xpr} EIGEN_BENCH_PREDICATE_OP=${operation}) +endforeach() +eigen_add_benchmark(bench_predicate_block_iszero bench_predicate_reductions.cpp) # Separate operations and storage families to avoid moving unrelated hot loops. foreach(operation All Any Count) string(TOLOWER ${operation} operation_name)
diff --git a/benchmarks/Core/bench_predicate_reductions.cpp b/benchmarks/Core/bench_predicate_reductions.cpp new file mode 100644 index 0000000..6a373dc --- /dev/null +++ b/benchmarks/Core/bench_predicate_reductions.cpp
@@ -0,0 +1,134 @@ +// SPDX-FileCopyrightText: The Eigen Authors +// SPDX-License-Identifier: MPL-2.0 + +// all(), any(), allFinite(), isZero() and isApproxToConstant(). The plain cases are full scans: no coefficient decides +// the result early. The "/first" cases are decided by coefficient 0, which measures the reads past the deciding one. +// Dynamic sizes hide the size and data pointer from the compiler on every iteration, as a caller whose sizes vary +// would; fixed sizes and a block of a matrix cover the other traversals. +// +// Each binary measures one family, selected by EIGEN_BENCH_PREDICATE_SCALAR and EIGEN_BENCH_PREDICATE_OP (dynamic), +// EIGEN_BENCH_PREDICATE_XPR and EIGEN_BENCH_PREDICATE_OP (fixed), or neither (block). + +#include <benchmark/benchmark.h> + +#include <Eigen/Core> + +#include <cstdint> +#include <limits> +#include <string> + +namespace Eigen { +namespace { + +using std::int16_t; + +enum Op { All, Any, AllFinite, IsZero, ApproxConstant }; + +template <typename T, Op op> +T fill() { + return op == Any || op == IsZero ? T(0) : op == ApproxConstant ? T(2) : T(1); +} + +// A coefficient that decides the result on its own. +template <typename T, Op op> +T decider() { + return op == All ? T(0) : op == AllFinite ? std::numeric_limits<T>::infinity() : op == ApproxConstant ? T(3) : T(1); +} + +template <Op op, typename Xpr> +bool apply(const Xpr& x) { + using Scalar = typename Xpr::Scalar; + switch (op) { + case All: + return x.all(); + case Any: + return x.any(); + case AllFinite: + return x.allFinite(); + case IsZero: + return x.isZero(); + default: + return x.isApproxToConstant(Scalar(2)); + } +} + +template <typename T, Op op, bool First> +void BM_Dynamic(benchmark::State& state) { + const Index n = state.range(0); + ArrayX<T> a = ArrayX<T>::Constant(n, fill<T, op>()); + if (First) a(0) = decider<T, op>(); + for (auto _ : state) { + Index m = n; + benchmark::DoNotOptimize(m); + const T* p = a.data(); + benchmark::DoNotOptimize(p); + bool r = apply<op>(Map<const ArrayX<T>>(p, m)); + benchmark::DoNotOptimize(r); + } + state.SetItemsProcessed(state.iterations() * n); +} + +template <typename Xpr, Op op> +void BM_Fixed(benchmark::State& state) { + Xpr a = Xpr::Constant(fill<typename Xpr::Scalar, op>()); + for (auto _ : state) { + benchmark::DoNotOptimize(a); + bool r = apply<op>(a); + benchmark::DoNotOptimize(r); + } +} + +// isZero() on the top rows of a 64-column matrix: one inner vector per column, no linear access. +void BM_BlockIsZero(benchmark::State& state) { + const Index n = state.range(0); + MatrixXf m = MatrixXf::Zero(n + 3, 64); + for (auto _ : state) { + Index k = n; + benchmark::DoNotOptimize(k); + benchmark::DoNotOptimize(m.data()); + bool r = m.topRows(k).isZero(); + benchmark::DoNotOptimize(r); + } + state.SetItemsProcessed(state.iterations() * n * 64); +} + +using Array16f = Array<float, 16, 1>; +using Array33f = Array<float, 33, 1>; + +const char* const opNames[] = {"all", "any", "allFinite", "isZero", "isApproxToConstant"}; + +#if defined(EIGEN_BENCH_PREDICATE_SCALAR) +const bool registered = [] { + using T = EIGEN_BENCH_PREDICATE_SCALAR; + constexpr Op op = EIGEN_BENCH_PREDICATE_OP; + const std::string name = std::string(opNames[op]) + "<" EIGEN_MAKESTRING(EIGEN_BENCH_PREDICATE_SCALAR) ">"; + benchmark::RegisterBenchmark(name.c_str(), BM_Dynamic<T, op, false>) + ->DenseRange(1, 16) + ->Arg(32) + ->Arg(64) + ->Arg(1024) + ->Arg(16384); + benchmark::RegisterBenchmark((name + "/first").c_str(), BM_Dynamic<T, op, true>) + ->Arg(16) + ->Arg(64) + ->Arg(256) + ->Arg(1024) + ->Arg(16384); + return true; +}(); +#elif defined(EIGEN_BENCH_PREDICATE_XPR) +const bool registered = [] { + constexpr Op op = EIGEN_BENCH_PREDICATE_OP; + const std::string name = std::string(opNames[op]) + "<" EIGEN_MAKESTRING(EIGEN_BENCH_PREDICATE_XPR) ">"; + benchmark::RegisterBenchmark(name.c_str(), BM_Fixed<EIGEN_BENCH_PREDICATE_XPR, op>); + return true; +}(); +#else +const bool registered = [] { + benchmark::RegisterBenchmark("blockIsZero<float>", BM_BlockIsZero)->Arg(4)->Arg(16)->Arg(64)->Arg(256); + return true; +}(); +#endif + +} // namespace +} // namespace Eigen
diff --git a/test/fuzzy_matrix_checks.cpp b/test/fuzzy_matrix_checks.cpp index 05ce641..f450c7f 100644 --- a/test/fuzzy_matrix_checks.cpp +++ b/test/fuzzy_matrix_checks.cpp
@@ -100,12 +100,14 @@ template <typename Scalar, int Order> void check_fuzzy_boundaries() { using Mat = Matrix<Scalar, Dynamic, Dynamic, Order>; - using Visitor = internal::fuzzy_constant_visitor<Scalar, false>; - STATIC_CHECK((internal::visit_impl<Mat, Visitor, true>::LinearAccess)); - STATIC_CHECK((!internal::visit_impl<Block<Mat>, Visitor, true>::LinearAccess)); + using Predicate = internal::fuzzy_constant_predicate<Scalar, false>; + STATIC_CHECK((internal::predicate_reduction<Mat, Predicate>::LinearAccess)); + STATIC_CHECK((!internal::predicate_reduction<Block<Mat>, Predicate>::LinearAccess)); constexpr int packetSize = internal::packet_traits<Scalar>::size; - for (Index size : - {Index(1), Index(packetSize - 1), Index(packetSize), Index(packetSize + 1), Index(2 * packetSize + 1)}) { + // Sizes up to 9 * packetSize + 3 run the long path's 8-packet blocks, 4-packet block and remainder for every packet + // size. + for (Index size : {Index(1), Index(packetSize - 1), Index(packetSize), Index(packetSize + 1), + Index(2 * packetSize + 1), Index(4 * packetSize), Index(9 * packetSize + 3)}) { // Unaligned linear access, including a packet crossing a column boundary. std::vector<Scalar> data(3 * size + 1); Map<Mat, Unaligned> matrix(data.data() + 1, 3, size); @@ -120,17 +122,20 @@ VERIFY(!matrix.isApproxToConstant(Scalar(0.5))); } } - for (Index inner : {Index(1), Index(packetSize), Index(packetSize + 1), Index(2 * packetSize + 1)}) { - Mat storage = Mat::Zero(inner + 2, inner + 2); - auto block = storage.block(1, 1, inner, inner); - for (Index k = 0; k < inner * inner; ++k) { + for (Index inner : {Index(1), Index(packetSize), Index(packetSize + 1), Index(2 * packetSize + 1), + Index(4 * packetSize), Index(9 * packetSize + 3)}) { + // Three inner vectors of the given length. + const Index rows = Order == RowMajor ? 3 : inner, cols = Order == RowMajor ? inner : 3; + Mat storage = Mat::Zero(rows + 2, cols + 2); + auto block = storage.block(1, 1, rows, cols); + for (Index k = 0; k < rows * cols; ++k) { block.setZero(); VERIFY(block.isZero()); - block(k / inner, k % inner) = Scalar(1); + block(k / cols, k % cols) = Scalar(1); VERIFY(!block.isZero()); block.setConstant(Scalar(0.5)); VERIFY(block.isApproxToConstant(Scalar(0.5))); - block(k / inner, k % inner) = Scalar(1); + block(k / cols, k % cols) = Scalar(1); VERIFY(!block.isApproxToConstant(Scalar(0.5))); } } @@ -162,6 +167,30 @@ VERIFY(!strided.isZero()); } +// The packet test of isApproxToConstant() holds for 0 <= prec <= 1 and a finite constant; other precisions and +// constants fall back to scalar tests. Coefficients of either sign, signed zeros and non-finite values at every +// position of a vector long enough for packet blocks and a scalar tail. +template <typename Scalar> +void check_fuzzy_constant_signs() { + const Scalar inf = NumTraits<Scalar>::infinity(), nan = NumTraits<Scalar>::quiet_NaN(); + constexpr int packetSize = internal::packet_traits<Scalar>::size; + Matrix<Scalar, Dynamic, 1> v(9 * packetSize + 3); + for (Scalar value : {Scalar(0.5), Scalar(-0.5), Scalar(0), -Scalar(0), inf, -inf, NumTraits<Scalar>::highest()}) { + for (Scalar precision : {Scalar(0), Scalar(1e-3), Scalar(0.5), Scalar(1), Scalar(1.5), Scalar(3), inf}) { + for (Scalar x : + {-value, value * Scalar(1.25), value * Scalar(-0.25), Scalar(0), -Scalar(0), inf, -inf, nan, Scalar(1)}) { + for (Index k = 0; k < v.size(); ++k) { + v.setConstant(value); + v(k) = x; + bool expected = true; + for (Index i = 0; i < v.size(); ++i) expected &= internal::isApprox(v(i), value, precision); + VERIFY_IS_EQUAL(v.isApproxToConstant(value, precision), expected); + } + } + } + } +} + template <typename Scalar, int Order> void check_fuzzy_subnormals() { const ScopedFlushToZero gradual(FlushToZeroMode::None); @@ -169,9 +198,9 @@ if (subnormalInputProbe<Scalar>() == Scalar(0) || underflowProbe<Scalar>() == Scalar(0)) return; const Scalar subnormal = (std::numeric_limits<Scalar>::min)() / Scalar(4); using Mat = Matrix<Scalar, Dynamic, Dynamic, Order>; - using Visitor = internal::fuzzy_constant_visitor<Scalar, false>; + using Predicate = internal::fuzzy_constant_predicate<Scalar, false>; STATIC_CHECK( - (!EIGEN_ARCH_ARM || !std::is_same<Scalar, float>::value || !internal::functor_traits<Visitor>::PacketAccess)); + (!EIGEN_ARCH_ARM || !std::is_same<Scalar, float>::value || !internal::functor_traits<Predicate>::PacketAccess)); for (Index size : {1, 2, 4, 8, 9, 17}) { Matrix<Scalar, Dynamic, 1> v = Matrix<Scalar, Dynamic, 1>::Zero(size); for (Index k = 0; k < size; ++k) { @@ -218,5 +247,7 @@ CALL_SUBTEST_4((check_fuzzy_loops<double, RowMajor>())); CALL_SUBTEST_4((check_fuzzy_subnormals<double, RowMajor>())); CALL_SUBTEST_5(check_fuzzy_custom_scalars()); + CALL_SUBTEST_6(check_fuzzy_constant_signs<float>()); + CALL_SUBTEST_6(check_fuzzy_constant_signs<double>()); } }
diff --git a/test/visitor.cpp b/test/visitor.cpp index 66d3ed9..3c0153f 100644 --- a/test/visitor.cpp +++ b/test/visitor.cpp
@@ -392,12 +392,16 @@ template <typename Scalar, int Options> void checkBooleanVisitorTraversal() { - STATIC_CHECK((internal::visitor_has_linear_access<internal::all_visitor<Scalar>>::value)); - STATIC_CHECK((internal::visitor_has_linear_access<internal::any_visitor<Scalar>>::value)); - STATIC_CHECK((internal::visitor_has_linear_access<internal::count_visitor<Scalar>>::value)); constexpr int PacketSize = internal::packet_traits<Scalar>::size; using MatrixType = Matrix<Scalar, Dynamic, Dynamic, Options>; - const Index sizes[] = {0, 1, PacketSize - 1, PacketSize, PacketSize + 1, 2 * PacketSize + 1}; + using Reduction = internal::predicate_reduction<MatrixType, internal::nonzero_predicate<Scalar>>; + STATIC_CHECK((Reduction::LinearAccess)); + STATIC_CHECK((!internal::predicate_reduction<Block<MatrixType>, internal::nonzero_predicate<Scalar>>::LinearAccess)); + STATIC_CHECK((internal::visitor_has_linear_access<internal::count_visitor<Scalar>>::value)); + // 4 and 9 packets plus a tail reach the blocked loop, the merged remainder packets, and the scalar tail; 131 reaches + // two 64-coefficient blocks and a tail without packets. + const Index sizes[] = { + 0, 1, PacketSize - 1, PacketSize, PacketSize + 1, 2 * PacketSize + 1, 4 * PacketSize, 9 * PacketSize + 3, 131}; for (Index inner : sizes) { MatrixType matrix(Options == RowMajor ? 3 : inner, Options == RowMajor ? inner : 3); for (int value = 0; value <= 1; ++value) { @@ -449,7 +453,11 @@ const Index inner = Vectorize ? 1 : 7; const Index outer = 35; MatrixType matrix(Options == RowMajor ? outer : inner, Options == RowMajor ? inner : outer); - for (Index index : {Index(0), Index(1), matrix.size() - 1}) { + // A short-circuit stops at the end of the block holding the deciding coefficient: up to eight packets, or 64 + // coefficients without packets. + const Index blockSize = Vectorize ? 8 * internal::packet_traits<float>::size : 64; + auto maxReads = [&](Index index) { return (std::min)(matrix.size(), index + blockSize); }; + for (Index index : {Index(0), Index(1), matrix.size() / 2, matrix.size() - 1}) { Index reads = 0; auto counted = matrix.unaryExpr(CountVisitorReads<Vectorize>{&reads}); for (int value = 0; value <= 1; ++value) { @@ -457,11 +465,11 @@ matrix(index) = float(!value); reads = 0; VERIFY_IS_EQUAL(value ? counted.all() : counted.any(), !value); - VERIFY_IS_EQUAL(reads, index + 1); + VERIFY(reads > index && reads <= maxReads(index)); reads = 0; auto block = counted.block(0, 0, matrix.rows(), matrix.cols()); VERIFY_IS_EQUAL(value ? block.all() : block.any(), !value); - VERIFY_IS_EQUAL(reads, index + 1); + VERIFY(reads > index && reads <= maxReads(index)); } } } @@ -665,7 +673,11 @@ CALL_SUBTEST_14((checkBooleanVisitorTraversal<float, ColMajor>())); CALL_SUBTEST_14((checkBooleanVisitorTraversal<double, RowMajor>())); CALL_SUBTEST_14((checkBooleanVisitorTraversal<int, ColMajor>())); + CALL_SUBTEST_14((checkBooleanVisitorTraversal<int64_t, RowMajor>())); + CALL_SUBTEST_14((checkBooleanVisitorTraversal<Eigen::half, ColMajor>())); CALL_SUBTEST_14((checkBooleanVisitorTraversal<std::complex<float>, RowMajor>())); + CALL_SUBTEST_14((checkBooleanVisitorTraversal<int16_t, ColMajor>())); + CALL_SUBTEST_14((checkBooleanVisitorTraversal<long double, RowMajor>())); CALL_SUBTEST_14(checkVisitorShortCircuit<ColMajor>()); CALL_SUBTEST_14(checkVisitorShortCircuit<RowMajor>()); CALL_SUBTEST_14((checkVisitorShortCircuit<ColMajor, true>()));