diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h index 79beb7bbc..230197f1f 100644 --- a/Eigen/src/Core/AssignEvaluator.h +++ b/Eigen/src/Core/AssignEvaluator.h @@ -631,6 +631,41 @@ struct dense_assignment_loop_impl using head_loop = unaligned_dense_assignment_loop; using tail_loop = unaligned_dense_assignment_loop; + // All inner slices share one alignment offset: the loop bounds are outer-invariant and stay + // hoisted out of the outer loop. Keeping that loop free of per-slice bookkeeping is what makes + // slice vectorization competitive with compiler-vectorized scalar code when slices are short. + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE constexpr void runInvariant(Kernel& kernel, Index alignedStart, + Index innerSize, Index outerSize) { + const Index alignedEnd = alignedStart + numext::round_down(innerSize - alignedStart, PacketSize); + for (Index outer = 0; outer < outerSize; ++outer) { + head_loop::run(kernel, outer, 0, alignedStart); + + // do the vectorizable part of the assignment + for (Index inner = alignedStart; inner < alignedEnd; inner += PacketSize) + kernel.template assignPacketByOuterInner(outer, inner); + + tail_loop::run(kernel, outer, alignedEnd, innerSize); + } + } + +#if EIGEN_UNALIGNED_VECTORIZE + // The slice alignment offset varies from one outer index to the next. Unaligned stores with + // outer-invariant bounds beat chasing the aligned position of every slice. + using unaligned_tail_loop = + unaligned_dense_assignment_loop; + + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE constexpr void runUnaligned(Kernel& kernel, Index innerSize, + Index outerSize) { + const Index packetEnd = numext::round_down(innerSize, PacketSize); + for (Index outer = 0; outer < outerSize; ++outer) { + for (Index inner = 0; inner < packetEnd; inner += PacketSize) + kernel.template assignPacketByOuterInner(outer, inner); + + unaligned_tail_loop::run(kernel, outer, packetEnd, innerSize); + } + } +#endif + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE constexpr void run(Kernel& kernel) { const Scalar* dst_ptr = kernel.dstDataPtr(); const Index innerSize = kernel.innerSize(); @@ -638,6 +673,13 @@ struct dense_assignment_loop_impl const Index alignedStep = Alignable ? (PacketSize - kernel.outerStride() % PacketSize) % PacketSize : 0; Index alignedStart = ((!Alignable) || DstIsAligned) ? 0 : internal::first_aligned(dst_ptr, innerSize); + if (alignedStep == 0) { + runInvariant(kernel, alignedStart, innerSize, outerSize); + return; + } +#if EIGEN_UNALIGNED_VECTORIZE + runUnaligned(kernel, innerSize, outerSize); +#else for (Index outer = 0; outer < outerSize; ++outer) { const Index alignedEnd = alignedStart + numext::round_down(innerSize - alignedStart, PacketSize); @@ -651,6 +693,7 @@ struct dense_assignment_loop_impl alignedStart = numext::mini((alignedStart + alignedStep) % PacketSize, innerSize); } +#endif } }; diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h index c8acaa3a3..889f46adb 100644 --- a/Eigen/src/Core/CoreEvaluators.h +++ b/Eigen/src/Core/CoreEvaluators.h @@ -1547,7 +1547,15 @@ struct unary_evaluator> enum { CoeffReadCost = evaluator::CoeffReadCost, LinearAccessMask = XprType::IsVectorAtCompileTime ? LinearAccessBit : 0, - Flags = (evaluator::Flags & (HereditaryBits | LinearAccessMask) & ~RowMajorBit) | + // The packet paths below load from a single copy of the nested expression, so they are valid + // exactly when a packet cannot cross a replication boundary: the inner (storage-order) + // direction must not be replicated. The outer coordinate's modulo then maps any packet into + // the nested expression unchanged. When the inner direction is replicated, serving a packet + // would need a broadcast (or a wrap-around load) the methods below do not perform. + InnerFactor = traits::IsRowMajor ? ColFactor : RowFactor, + MaskPacketAccessBit = InnerFactor == 1 ? PacketAccessBit : 0, + Flags = (evaluator::Flags & (HereditaryBits | LinearAccessMask | MaskPacketAccessBit) & + ~RowMajorBit) | (traits::Flags & RowMajorBit), Alignment = evaluator::Alignment diff --git a/benchmarks/Core/bench_broadcasting.cpp b/benchmarks/Core/bench_broadcasting.cpp index 71af7e774..142b82729 100644 --- a/benchmarks/Core/bench_broadcasting.cpp +++ b/benchmarks/Core/bench_broadcasting.cpp @@ -2,6 +2,11 @@ // // Tests vectorwise reductions (sum, mean, norm, minCoeff, maxCoeff) and // broadcasting arithmetic (rowwise += vec, colwise -= vec, rowwise *= vec). +// The BroadcastSubExp pair isolates Eigen's packet path through the Replicate +// evaluator: the scalar fallback calls libm exp per element, so unlike the +// plain add/mul broadcasts it cannot be rescued by compiler auto-vectorization. +// The PerCol variant performs the same operation with an explicit column loop +// and serves as an upper-bound reference. // SPDX-FileCopyrightText: The Eigen Authors // SPDX-License-Identifier: MPL-2.0 @@ -158,6 +163,42 @@ static void BM_RowwiseBroadcastMul(benchmark::State& state) { state.SetBytesProcessed(state.iterations() * rows * cols * sizeof(Scalar) * 2); } +// Broadcast feeding a transcendental (softmax-style normalization): +// out = exp(m - v replicated across columns). +template +static void BM_ColwiseBroadcastSubExp(benchmark::State& state) { + const Index rows = state.range(0); + const Index cols = state.range(1); + using Mat = Matrix; + using Vec = Matrix; + Mat m = Mat::Random(rows, cols); + Vec v = Vec::Random(rows); + Mat out(rows, cols); + for (auto _ : state) { + out = (m.colwise() - v).array().exp(); + benchmark::DoNotOptimize(out.data()); + } + state.SetBytesProcessed(state.iterations() * rows * cols * sizeof(Scalar) * 2); +} + +// Same operation with an explicit per-column loop: upper-bound reference for +// BM_ColwiseBroadcastSubExp. +template +static void BM_PerColBroadcastSubExp(benchmark::State& state) { + const Index rows = state.range(0); + const Index cols = state.range(1); + using Mat = Matrix; + using Vec = Matrix; + Mat m = Mat::Random(rows, cols); + Vec v = Vec::Random(rows); + Mat out(rows, cols); + for (auto _ : state) { + for (Index j = 0; j < cols; ++j) out.col(j) = (m.col(j) - v).array().exp(); + benchmark::DoNotOptimize(out.data()); + } + state.SetBytesProcessed(state.iterations() * rows * cols * sizeof(Scalar) * 2); +} + // --- Size configurations --- // clang-format off // Square matrices; tall-thin (many rows, few cols); short-wide (few rows, many cols). @@ -176,6 +217,8 @@ BENCHMARK(BM_RowwiseNorm) BROADCAST_SIZES ->Name("RowwiseNorm_float"); BENCHMARK(BM_RowwiseBroadcastAdd) BROADCAST_SIZES ->Name("RowwiseBroadcastAdd_float"); BENCHMARK(BM_ColwiseBroadcastAdd) BROADCAST_SIZES ->Name("ColwiseBroadcastAdd_float"); BENCHMARK(BM_RowwiseBroadcastMul) BROADCAST_SIZES ->Name("RowwiseBroadcastMul_float"); +BENCHMARK(BM_ColwiseBroadcastSubExp) BROADCAST_SIZES ->Name("ColwiseBroadcastSubExp_float"); +BENCHMARK(BM_PerColBroadcastSubExp) BROADCAST_SIZES ->Name("PerColBroadcastSubExp_float"); // --- Register: double --- BENCHMARK(BM_ColwiseSum) BROADCAST_SIZES ->Name("ColwiseSum_double"); @@ -188,6 +231,8 @@ BENCHMARK(BM_RowwiseNorm) BROADCAST_SIZES ->Name("RowwiseNorm_double"); BENCHMARK(BM_RowwiseBroadcastAdd) BROADCAST_SIZES ->Name("RowwiseBroadcastAdd_double"); BENCHMARK(BM_ColwiseBroadcastAdd) BROADCAST_SIZES ->Name("ColwiseBroadcastAdd_double"); BENCHMARK(BM_RowwiseBroadcastMul) BROADCAST_SIZES ->Name("RowwiseBroadcastMul_double"); +BENCHMARK(BM_ColwiseBroadcastSubExp) BROADCAST_SIZES ->Name("ColwiseBroadcastSubExp_double"); +BENCHMARK(BM_PerColBroadcastSubExp) BROADCAST_SIZES ->Name("PerColBroadcastSubExp_double"); #undef BROADCAST_SIZES // clang-format on diff --git a/test/array_replicate.cpp b/test/array_replicate.cpp index d29223a41..0e7587c41 100644 --- a/test/array_replicate.cpp +++ b/test/array_replicate.cpp @@ -62,6 +62,69 @@ void replicate(const MatrixType& m) { VERIFY_IS_APPROX(vx1, v1.colwise().replicate(f2)); } +// A Replicate packet cannot cross a replication boundary, so the evaluator may serve packets +// exactly when the inner (storage-order) direction is not replicated; without the flag, colwise +// and rowwise broadcast operations silently fall back to scalar traversal. +template +void check_replicate_evaluator_flags() { + // Storage orders are pinned so the checks keep their meaning under EIGEN_DEFAULT_TO_ROW_MAJOR. + typedef Matrix Mat; + typedef Matrix RowMat; + typedef Matrix Vec; + typedef Matrix RowVec; + enum { BasePacket = int(internal::evaluator::Flags) & PacketAccessBit }; + + // Inner direction not replicated: packets serve from a single copy of the argument. The first + // two are the shapes colwise (col-major) and rowwise (row-major) operations expand to. + STATIC_CHECK((int(internal::evaluator >::Flags) & PacketAccessBit) == int(BasePacket)); + STATIC_CHECK((int(internal::evaluator >::Flags) & PacketAccessBit) == int(BasePacket)); + STATIC_CHECK((int(internal::evaluator >::Flags) & PacketAccessBit) == int(BasePacket)); + STATIC_CHECK((int(internal::evaluator >::Flags) & PacketAccessBit) == int(BasePacket)); + STATIC_CHECK(int(internal::evaluator >::Alignment) == + int(internal::evaluator::Alignment)); + + // Replicated inner direction (or factors unknown at compile time): a packet could cross a copy + // boundary, so there is no packet access. + STATIC_CHECK((int(internal::evaluator >::Flags) & PacketAccessBit) == 0); + STATIC_CHECK((int(internal::evaluator >::Flags) & PacketAccessBit) == 0); + STATIC_CHECK((int(internal::evaluator >::Flags) & PacketAccessBit) == 0); +} + +// Exercise the (possibly vectorized) broadcast kernels with sizes that have partial-packet tails. +template +void replicate_broadcasts(Index rows, Index cols) { + typedef Matrix Mat; + typedef Matrix Vec; + typedef Matrix RowVec; + + Mat m = Mat::Random(rows, cols); + Vec v = Vec::Random(rows); + RowVec rv = RowVec::Random(cols); + + Mat c = m; + c.colwise() += v; + for (Index j = 0; j < cols; ++j) + for (Index i = 0; i < rows; ++i) VERIFY_IS_EQUAL(c(i, j), Scalar(m(i, j) + v(i))); + + Mat r = m; + r.rowwise() += rv; + for (Index j = 0; j < cols; ++j) + for (Index i = 0; i < rows; ++i) VERIFY_IS_EQUAL(r(i, j), Scalar(m(i, j) + rv(j))); + + Mat h = m.template replicate<1, 3>(); + for (Index j = 0; j < 3 * cols; ++j) + for (Index i = 0; i < rows; ++i) VERIFY_IS_EQUAL(h(i, j), m(i, j % cols)); + + Mat ver = m.template replicate<3, 1>(); + for (Index j = 0; j < cols; ++j) + for (Index i = 0; i < 3 * rows; ++i) VERIFY_IS_EQUAL(ver(i, j), m(i % rows, j)); + + // A replicate nested inside a larger coefficient-wise expression. + Mat sum = m + v.rowwise().replicate(cols); + for (Index j = 0; j < cols; ++j) + for (Index i = 0; i < rows; ++i) VERIFY_IS_EQUAL(sum(i, j), Scalar(m(i, j) + v(i))); +} + EIGEN_DECLARE_TEST(array_replicate) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(replicate(Matrix())); @@ -70,5 +133,11 @@ EIGEN_DECLARE_TEST(array_replicate) { CALL_SUBTEST_4(replicate(Vector4f())); CALL_SUBTEST_5(replicate(VectorXf(16))); CALL_SUBTEST_6(replicate(VectorXcd(10))); + CALL_SUBTEST_7(check_replicate_evaluator_flags()); + CALL_SUBTEST_7(replicate_broadcasts(internal::random(1, 64), internal::random(1, 64))); + CALL_SUBTEST_7(replicate_broadcasts(17, 19)); + CALL_SUBTEST_8(check_replicate_evaluator_flags()); + CALL_SUBTEST_8(replicate_broadcasts(internal::random(1, 64), internal::random(1, 64))); + CALL_SUBTEST_8(replicate_broadcasts >(9, 5)); } }