Core: Share empty-input dispatch for reductions libeigen/eigen!3358
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index ea2c729..933be95 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h
@@ -87,7 +87,8 @@ using HasWritableRealView = bool_constant<!NumTraits<Scalar>::IsComplex || complex_array_access<Scalar>::value>; EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Accumulator max_abs(const VectorType& vec) { - return vec.realView().template cast<Accumulator>().cwiseAbs().template maxCoeff<PropagateNaN>(); + return redux_with_identity(vec.realView().template cast<Accumulator>().cwiseAbs(), + scalar_max_op<Accumulator, Accumulator, PropagateNaN>(), 0); } EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Accumulator scaled_squared_norm(const VectorType& vec, @@ -502,9 +503,7 @@ struct lpNorm_selector<Derived, Infinity> { using RealScalar = typename NumTraits<typename traits<Derived>::Scalar>::Real; EIGEN_DEVICE_FUNC static inline RealScalar run(const MatrixBase<Derived>& m) { - if (Derived::SizeAtCompileTime == 0 || (Derived::SizeAtCompileTime == Dynamic && m.size() == 0)) - return RealScalar(0); - return m.cwiseAbs().maxCoeff(); + return redux_with_identity(m.cwiseAbs(), scalar_max_op<RealScalar, RealScalar>(), 0); } };
diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h index 972d01a..b07e049 100644 --- a/Eigen/src/Core/Redux.h +++ b/Eigen/src/Core/Redux.h
@@ -740,6 +740,36 @@ } }; +// Convert identity to Scalar only for empty inputs; otherwise reduce without inserting it. +// Tag dispatch avoids instantiating redux for statically empty inputs in C++14. +// TODO(C++17): replace tag dispatch with if constexpr. +template <typename Derived, typename Func, typename Identity> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE typename Derived::Scalar redux_with_identity_impl(const DenseBase<Derived>&, + const Func&, + Identity&& identity, + true_type) { + return typename Derived::Scalar(std::forward<Identity>(identity)); +} + +template <typename Derived, typename Func, typename Identity> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE typename Derived::Scalar redux_with_identity_impl(const DenseBase<Derived>& xpr, + const Func& func, + Identity&& identity, + false_type) { + EIGEN_IF_CONSTEXPR (Derived::SizeAtCompileTime == Dynamic) { + if (xpr.size() == 0) return typename Derived::Scalar(std::forward<Identity>(identity)); + } + return xpr.derived().redux(func); +} + +template <typename Derived, typename Func, typename Identity> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE typename Derived::Scalar redux_with_identity(const DenseBase<Derived>& xpr, + const Func& func, + Identity&& identity) { + return redux_with_identity_impl(xpr, func, std::forward<Identity>(identity), + bool_constant<Derived::MaxSizeAtCompileTime == 0>()); +} + } // end namespace internal /*************************************************************************** @@ -808,14 +838,7 @@ */ template <typename Derived> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar DenseBase<Derived>::sum() const { - EIGEN_IF_CONSTEXPR (MaxSizeAtCompileTime == 0 || SizeAtCompileTime == 0) { - return Scalar(0); - } else { - EIGEN_IF_CONSTEXPR (SizeAtCompileTime == Dynamic) { - if (size() == 0) return Scalar(0); - } - return derived().redux(Eigen::internal::scalar_sum_op<Scalar, Scalar>()); - } + return internal::redux_with_identity(*this, internal::scalar_sum_op<Scalar, Scalar>(), 0); } /** \returns the mean of all coefficients of *this @@ -843,14 +866,7 @@ */ template <typename Derived> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::traits<Derived>::Scalar DenseBase<Derived>::prod() const { - EIGEN_IF_CONSTEXPR (MaxSizeAtCompileTime == 0 || SizeAtCompileTime == 0) { - return Scalar(1); - } else { - EIGEN_IF_CONSTEXPR (SizeAtCompileTime == Dynamic) { - if (size() == 0) return Scalar(1); - } - return derived().redux(Eigen::internal::scalar_product_op<Scalar>()); - } + return internal::redux_with_identity(*this, internal::scalar_product_op<Scalar>(), 1); } /** \returns the trace of \c *this, i.e. the sum of the coefficients on the main diagonal.
diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 2a43961..654d78c 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h
@@ -148,7 +148,14 @@ } }; -template <typename VectorType, std::enable_if_t<VectorType::IsVectorAtCompileTime, int> = 0> +// Avoid instantiating block reductions for statically empty inputs, including in C++14. +template <typename Derived, std::enable_if_t<Derived::MaxSizeAtCompileTime == 0, int> = 0> +typename Derived::RealScalar stable_norm_impl(const Derived&) { + return typename Derived::RealScalar(0); +} + +template <typename VectorType, + std::enable_if_t<VectorType::IsVectorAtCompileTime && VectorType::MaxSizeAtCompileTime != 0, int> = 0> typename VectorType::RealScalar stable_norm_impl(const VectorType& vec) { using std::sqrt; @@ -166,7 +173,8 @@ return RealScalar(scale * sqrt(ssq)); } -template <typename MatrixType, std::enable_if_t<!MatrixType::IsVectorAtCompileTime, int> = 0> +template <typename MatrixType, + std::enable_if_t<!MatrixType::IsVectorAtCompileTime && MatrixType::MaxSizeAtCompileTime != 0, int> = 0> typename MatrixType::RealScalar stable_norm_impl(const MatrixType& mat) { using std::sqrt; @@ -329,10 +337,9 @@ template <typename Derived> inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real MatrixBase<Derived>::hypotNorm() const { using Accumulator = typename internal::stable_norm_accumulator<RealScalar>::type; - if (size() == 0) return RealScalar(0); // Component reduction avoids rounded complex magnitudes and permits promoted accumulation. - return RealScalar( - derived().realView().template cast<Accumulator>().cwiseAbs().redux(internal::scalar_hypot_op<Accumulator>())); + return RealScalar(internal::redux_with_identity(derived().realView().template cast<Accumulator>().cwiseAbs(), + internal::scalar_hypot_op<Accumulator>(), 0)); } } // end namespace Eigen
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 95efe61..0582a1e 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -208,6 +208,10 @@ namespace internal { template <typename Lhs, typename Rhs> struct product_type; +template <typename Derived, typename Func, typename Identity> +EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE typename Derived::Scalar redux_with_identity(const DenseBase<Derived>& xpr, + const Func& func, + Identity&& identity); /** \internal * \class product_evaluator
diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index 5373ca0..365adc6 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h
@@ -100,6 +100,23 @@ template <typename Func> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::result_of<Func(Scalar, Scalar)>::type redux( const Func& func) const { + return redux_impl(func, internal::bool_constant<MatrixType::MaxSizeAtCompileTime == 0>()); + } + + protected: + template <typename Func> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::result_of<Func(Scalar, Scalar)>::type redux_impl( + const Func& func, std::true_type) const { + // Reduce the appended coefficients, not the empty underlying expression. + return Base::redux(func); + } + + template <typename Func> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename internal::result_of<Func(Scalar, Scalar)>::type redux_impl( + const Func& func, std::false_type) const { + EIGEN_IF_CONSTEXPR (MatrixType::SizeAtCompileTime == Dynamic) { + if (m_matrix.size() == 0) return Base::redux(func); + } // One appended coefficient per column (Vertical) or per row (Horizontal) of m_matrix. using Ones = std::conditional_t<int(Direction) == Vertical, typename internal::plain_row_type<MatrixType>::type, typename internal::plain_col_type<MatrixType>::type>; @@ -107,7 +124,6 @@ Ones::Ones(int(Direction) == Vertical ? m_matrix.cols() : m_matrix.rows()).redux(func)); } - protected: typename MatrixType::Nested m_matrix; };
diff --git a/test/geo_homogeneous.cpp b/test/geo_homogeneous.cpp index fc88086..155ae86 100644 --- a/test/geo_homogeneous.cpp +++ b/test/geo_homogeneous.cpp
@@ -40,6 +40,7 @@ VERIFY_IS_APPROX(v0, hv0.hnormalized()); VERIFY_IS_APPROX(v0.homogeneous().sum(), hv0.sum()); + VERIFY_IS_APPROX(v0.homogeneous().prod(), hv0.prod()); VERIFY_IS_APPROX(v0.homogeneous().minCoeff(), hv0.minCoeff()); VERIFY_IS_APPROX(v0.homogeneous().maxCoeff(), hv0.maxCoeff()); @@ -138,7 +139,59 @@ } } +template <typename Expression> +void check_empty_homogeneous_reductions(const Expression& expression, Index appended_ones) { + using Scalar = typename Expression::Scalar; + VERIFY_IS_EQUAL(expression.size(), appended_ones); + VERIFY_IS_EQUAL(expression.sum(), Scalar(appended_ones)); + VERIFY_IS_EQUAL(expression.prod(), Scalar(1)); + if (appended_ones > 0) { + VERIFY_IS_EQUAL(expression.minCoeff(), Scalar(1)); + VERIFY_IS_EQUAL(expression.maxCoeff(), Scalar(1)); + // This associative functor has identity -1, not a guessed sum/product identity. + VERIFY_IS_EQUAL(expression.redux([](Scalar a, Scalar b) { return a + b + Scalar(1); }), + Scalar(2 * appended_ones - 1)); + } else { + VERIFY_RAISES_ASSERT(expression.minCoeff()); + VERIFY_RAISES_ASSERT(expression.maxCoeff()); + } +} + +template <typename Scalar> +void homogeneous_empty_reductions() { + using DynamicMatrix = Matrix<Scalar, Dynamic, Dynamic>; + check_empty_homogeneous_reductions(Matrix<Scalar, 0, 1>().homogeneous(), 1); + check_empty_homogeneous_reductions(Matrix<Scalar, 1, 0>().homogeneous(), 1); + check_empty_homogeneous_reductions(Matrix<Scalar, Dynamic, 1>(0).homogeneous(), 1); + check_empty_homogeneous_reductions(Matrix<Scalar, 1, Dynamic>(0).homogeneous(), 1); + check_empty_homogeneous_reductions(Matrix<Scalar, 0, 3>().colwise().homogeneous(), 3); + check_empty_homogeneous_reductions(Matrix<Scalar, 0, 3>().rowwise().homogeneous(), 0); + check_empty_homogeneous_reductions(Matrix<Scalar, 3, 0>().colwise().homogeneous(), 0); + check_empty_homogeneous_reductions(Matrix<Scalar, 3, 0>().rowwise().homogeneous(), 3); + check_empty_homogeneous_reductions(Matrix<Scalar, 0, 0>().colwise().homogeneous(), 0); + check_empty_homogeneous_reductions(Matrix<Scalar, 0, 0>().rowwise().homogeneous(), 0); + check_empty_homogeneous_reductions(DynamicMatrix(0, 3).colwise().homogeneous(), 3); + check_empty_homogeneous_reductions(DynamicMatrix(0, 3).rowwise().homogeneous(), 0); + check_empty_homogeneous_reductions(DynamicMatrix(3, 0).colwise().homogeneous(), 0); + check_empty_homogeneous_reductions(DynamicMatrix(3, 0).rowwise().homogeneous(), 3); + check_empty_homogeneous_reductions(Matrix<Scalar, Dynamic, Dynamic, 0, 0, 0>().colwise().homogeneous(), 0); + check_empty_homogeneous_reductions(Matrix<Scalar, Dynamic, Dynamic, 0, 0, 0>().rowwise().homogeneous(), 0); +} + +template <typename Scalar> +void homogeneous_empty_complex_reductions() { + VERIFY_IS_EQUAL((Matrix<Scalar, 0, 1>().homogeneous().sum()), Scalar(1)); + VERIFY_IS_EQUAL((Matrix<Scalar, 1, 0>().homogeneous().prod()), Scalar(1)); + VERIFY_IS_EQUAL((Matrix<Scalar, Dynamic, 1>(0).homogeneous().prod()), Scalar(1)); + VERIFY_IS_EQUAL((Matrix<Scalar, 1, Dynamic>(0).homogeneous().sum()), Scalar(1)); + VERIFY_IS_EQUAL((Matrix<Scalar, 0, 3>().colwise().homogeneous().sum()), Scalar(3)); + VERIFY_IS_EQUAL((Matrix<Scalar, 3, 0>().rowwise().homogeneous().sum()), Scalar(3)); +} + EIGEN_DECLARE_TEST(geo_homogeneous) { + CALL_SUBTEST_1(homogeneous_empty_reductions<float>()); + CALL_SUBTEST_2(homogeneous_empty_reductions<double>()); + CALL_SUBTEST_2(homogeneous_empty_complex_reductions<std::complex<double>>()); for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1((homogeneous<float, 1>())); CALL_SUBTEST_2((homogeneous<double, 3>()));
diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp index a428c1a..5b6d679 100644 --- a/test/stable_norm.cpp +++ b/test/stable_norm.cpp
@@ -230,6 +230,35 @@ VERIFY_IS_EQUAL(empty.stableNorm(), 0.0f); } +template <typename VectorType> +void check_empty_stable_normalization(const VectorType& input) { + VectorType empty(input); + using RealScalar = typename VectorType::RealScalar; + VERIFY_IS_EQUAL((internal::stable_normalization_dispatch<VectorType, RealScalar>::max_abs(empty)), RealScalar(0)); + VERIFY_IS_EQUAL(empty.stableNormalized().size(), 0); + empty.stableNormalize(); + VERIFY_IS_EQUAL(empty.size(), 0); +} + +template <typename Scalar> +void stable_normalization_empty_reductions() { + check_empty_stable_normalization(Matrix<Scalar, 0, 1>()); + check_empty_stable_normalization(Matrix<Scalar, 1, 0>()); + check_empty_stable_normalization(Matrix<Scalar, Dynamic, 1>(0)); + check_empty_stable_normalization(Matrix<Scalar, 1, Dynamic>(0)); +} + +template <typename RealScalar> +void stable_normalization_max_abs_special_values() { + using VectorType = Matrix<RealScalar, 2, 1>; + using Dispatch = internal::stable_normalization_dispatch<VectorType, RealScalar>; + VectorType special; + special << RealScalar(1), NumTraits<RealScalar>::quiet_NaN(); + VERIFY((numext::isnan)(Dispatch::max_abs(special))); + special << -NumTraits<RealScalar>::infinity(), RealScalar(1); + VERIFY_IS_EQUAL(Dispatch::max_abs(special), NumTraits<RealScalar>::infinity()); +} + template <typename RealScalar> void stable_normalize_extremes() { typedef Matrix<RealScalar, 2, 1> Vector2; @@ -832,6 +861,11 @@ EIGEN_DECLARE_TEST(stable_norm) { CALL_SUBTEST_1(test_empty<>()); + CALL_SUBTEST_1(stable_normalization_empty_reductions<float>()); + CALL_SUBTEST_2(stable_normalization_empty_reductions<double>()); + CALL_SUBTEST_2(stable_normalization_empty_reductions<std::complex<double>>()); + CALL_SUBTEST_1(stable_normalization_max_abs_special_values<float>()); + CALL_SUBTEST_2(stable_normalization_max_abs_special_values<double>()); for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_3(test_hypot<double>());
diff --git a/test/zerosized.cpp b/test/zerosized.cpp index ca3ebd4..7393a36 100644 --- a/test/zerosized.cpp +++ b/test/zerosized.cpp
@@ -9,6 +9,65 @@ // SPDX-License-Identifier: MPL-2.0 #include "main.h" +#include "MovableScalar.h" + +template <typename Scalar> +struct ReductionIdentity { + ReductionIdentity(double initialValue, int& lvalueCount, int& rvalueCount) + : value(initialValue), lvalueConversions(lvalueCount), rvalueConversions(rvalueCount) {} + ReductionIdentity(const ReductionIdentity&) = delete; + ReductionIdentity(ReductionIdentity&&) = delete; + + explicit operator Scalar() const& { + ++lvalueConversions; + return Scalar(value); + } + explicit operator Scalar() && { + ++rvalueConversions; + return Scalar(value); + } + + double value; + int& lvalueConversions; + int& rvalueConversions; +}; + +template <typename Scalar> +void reductionIdentityArguments() { + const internal::scalar_sum_op<Scalar, Scalar> sum; + const Matrix<Scalar, 0, 1> fixedEmpty; + const Matrix<Scalar, Dynamic, 1, 0, 0, 1> zeroCapacity(0); + Matrix<Scalar, Dynamic, 1> dynamic(0); + const double runtimeIdentity = 0.25; + VERIFY_IS_EQUAL(double(internal::redux_with_identity(fixedEmpty, sum, runtimeIdentity)), runtimeIdentity); + VERIFY_IS_EQUAL(double(internal::redux_with_identity(dynamic, sum, runtimeIdentity)), runtimeIdentity); + + int lvalueConversions = 0, rvalueConversions = 0; + ReductionIdentity<Scalar> identity(runtimeIdentity, lvalueConversions, rvalueConversions); + const auto& constIdentity = identity; + VERIFY_IS_EQUAL(double(internal::redux_with_identity(fixedEmpty, sum, identity)), runtimeIdentity); + VERIFY_IS_EQUAL(double(internal::redux_with_identity(zeroCapacity, sum, identity)), runtimeIdentity); + VERIFY_IS_EQUAL(double(internal::redux_with_identity(dynamic, sum, constIdentity)), runtimeIdentity); + VERIFY_IS_EQUAL(lvalueConversions, 3); + VERIFY_IS_EQUAL(rvalueConversions, 0); + + Matrix<Scalar, 2, 1> fixed; + fixed << Scalar(2), Scalar(3); + dynamic = fixed; + VERIFY_IS_EQUAL(double(internal::redux_with_identity(fixed, sum, identity)), 5.0); + VERIFY_IS_EQUAL(double(internal::redux_with_identity(dynamic, sum, std::move(identity))), 5.0); + VERIFY_IS_EQUAL(lvalueConversions, 3); + VERIFY_IS_EQUAL(rvalueConversions, 0); + + dynamic.resize(0); + VERIFY_IS_EQUAL(double(internal::redux_with_identity(dynamic, sum, std::move(identity))), runtimeIdentity); + VERIFY_IS_EQUAL( + double(internal::redux_with_identity( + fixedEmpty, sum, ReductionIdentity<Scalar>(runtimeIdentity, lvalueConversions, rvalueConversions))), + runtimeIdentity); + VERIFY_IS_EQUAL(lvalueConversions, 3); + VERIFY_IS_EQUAL(rvalueConversions, 2); +} template <typename MatrixType> void zeroReduction(const MatrixType& m) { @@ -19,6 +78,11 @@ VERIFY(m.sum() == 0); VERIFY(m.norm() == 0); VERIFY(m.squaredNorm() == 0); + EIGEN_IF_CONSTEXPR (!NumTraits<typename MatrixType::Scalar>::IsInteger) { + VERIFY(m.template lpNorm<Infinity>() == 0); + VERIFY(m.hypotNorm() == 0); + VERIFY(m.stableNorm() == 0); + } VERIFY(m.count() == 0); VERIFY(m.allFinite()); VERIFY(!m.hasNaN()); @@ -82,6 +146,17 @@ } } +// Complex extrema are unordered; check only reductions with defined empty identities. +template <typename MatrixType> +void zeroComplexReduction(const MatrixType& m) { + using Scalar = typename MatrixType::Scalar; + VERIFY_IS_EQUAL(m.sum(), Scalar(0)); + VERIFY_IS_EQUAL(m.prod(), Scalar(1)); + VERIFY_IS_EQUAL(m.template lpNorm<Infinity>(), 0.0); + VERIFY_IS_EQUAL(m.hypotNorm(), 0.0); + VERIFY_IS_EQUAL(m.stableNorm(), 0.0); +} + template <typename VectorType> void zeroSizedVector() { VectorType t1; @@ -98,6 +173,8 @@ } EIGEN_DECLARE_TEST(zerosized) { + reductionIdentityArguments<double>(); + reductionIdentityArguments<MovableScalar<double>>(); zeroSizedMatrix<Matrix2d>(); zeroSizedMatrix<Matrix3i>(); zeroSizedMatrix<Matrix<float, 2, Dynamic> >(); @@ -114,4 +191,9 @@ zeroSizedVector<VectorXf>(); zeroSizedVector<Matrix<float, 0, 1> >(); zeroSizedVector<Matrix<float, 1, 0> >(); + zeroComplexReduction(MatrixXcd(0, 3)); + zeroComplexReduction(MatrixXcd(3, 0)); + zeroComplexReduction(Matrix<std::complex<double>, 0, 0>()); + zeroComplexReduction(Matrix<std::complex<double>, 0, 3>()); + zeroComplexReduction(Matrix<std::complex<double>, 3, 0>()); }