From b17236cd0c50230510f406240dbf85ba1b6d0ef4 Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Sun, 4 Oct 2026 19:07:08 +0200 Subject: [PATCH 1/8] Improve SIMD assignment of broadcast expressions --- benchmark/CMakeLists.txt | 16 ++ benchmark/benchmark_compare.cpp | 233 +++++++++++++++++++++++++++++ include/xtensor/core/xiterator.hpp | 7 + include/xtensor/core/xstrides.hpp | 6 +- test/test_strided_assign.cpp | 47 ++++++ 5 files changed, 307 insertions(+), 2 deletions(-) create mode 100644 benchmark/benchmark_compare.cpp diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 0d0ce5766..669f82a29 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -100,8 +100,14 @@ if(DOWNLOAD_GBENCHMARK OR GBENCHMARK_SRC_DIR) set(GBENCHMARK_LIBRARIES benchmark) else() find_package(benchmark REQUIRED) + if(TARGET benchmark::benchmark) + set(GBENCHMARK_LIBRARIES benchmark::benchmark) + endif() endif() +find_package(Eigen3 QUIET NO_MODULE) +find_package(Armadillo QUIET) + find_package(xsimd) if (xsimd_FOUND) include_directories(${xsimd_INCLUDE_DIRS}) @@ -113,6 +119,7 @@ include_directories(${GBENCHMARK_INCLUDE_DIRS}) set(XTENSOR_BENCHMARK benchmark_assign.cpp + benchmark_compare.cpp benchmark_builder.cpp benchmark_container.cpp benchmark_creation.cpp @@ -134,6 +141,15 @@ set(XTENSOR_BENCHMARK set(XTENSOR_BENCHMARK_TARGET benchmark_xtensor) add_executable(${XTENSOR_BENCHMARK_TARGET} EXCLUDE_FROM_ALL ${XTENSOR_BENCHMARK} ${XTENSOR_HEADERS}) target_link_libraries(${XTENSOR_BENCHMARK_TARGET} PUBLIC xtensor ${GBENCHMARK_LIBRARIES}) +if(Eigen3_FOUND) + target_compile_definitions(${XTENSOR_BENCHMARK_TARGET} PRIVATE XTENSOR_BENCHMARK_USE_EIGEN) + target_link_libraries(${XTENSOR_BENCHMARK_TARGET} PRIVATE Eigen3::Eigen) +endif() +if(ARMADILLO_FOUND) + target_compile_definitions(${XTENSOR_BENCHMARK_TARGET} PRIVATE XTENSOR_BENCHMARK_USE_ARMADILLO) + target_include_directories(${XTENSOR_BENCHMARK_TARGET} PRIVATE ${ARMADILLO_INCLUDE_DIRS}) + target_link_libraries(${XTENSOR_BENCHMARK_TARGET} PRIVATE ${ARMADILLO_LIBRARIES}) +endif() if(XTENSOR_USE_TBB) target_compile_definitions(${XTENSOR_BENCHMARK_TARGET} PUBLIC XTENSOR_USE_TBB) diff --git a/benchmark/benchmark_compare.cpp b/benchmark/benchmark_compare.cpp new file mode 100644 index 000000000..8b5655ccf --- /dev/null +++ b/benchmark/benchmark_compare.cpp @@ -0,0 +1,233 @@ +/*************************************************************************** + * Copyright (c) 2016, Johan Mabille, Sylvain Corlay and Wolf Vollprecht * + * * + * Distributed under the terms of the BSD 3-Clause License. * + * * + * The full license is in the file LICENSE, distributed with this software. * + ****************************************************************************/ + +#include +#include +#include + +#include + +#include "xtensor/containers/xarray.hpp" +#include "xtensor/containers/xfixed.hpp" +#include "xtensor/core/xnoalias.hpp" +#include "xtensor/core/xoperation.hpp" +#include "xtensor/views/xstrided_view.hpp" + +#ifdef XTENSOR_BENCHMARK_USE_EIGEN +#include +#endif + +#ifdef XTENSOR_BENCHMARK_USE_ARMADILLO +#define ARMA_DONT_PRINT_FAST_MATH_WARNING +#include +#endif + +namespace xt::compare +{ + constexpr std::size_t dynamic_size = 64 * 64 * 16; + constexpr std::size_t fixed_size = 16 * 16; + + template + void init(X& x, Y& y, Z& z, std::size_t size) + { + for (std::size_t i = 0; i < size; ++i) + { + x[i] = 0.5 + double(i % 251) * 0.01; + y[i] = 0.25 - double(i % 127) * 0.02; + z[i] = 1.0 + double(i % 61) * 0.03; + } + } + + void linear_dynamic_xtensor(benchmark::State& state) + { + const std::vector shape = {64, 64, 16}; + xarray x, y, z, out; + x.resize(shape); y.resize(shape); z.resize(shape); out.resize(shape); + init(x, y, z, dynamic_size); + for (auto _ : state) + { + noalias(out) = (x + y * 0.5) * (z - 0.25) + y / (x + 2.0); + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void linear_dynamic_raw(benchmark::State& state) + { + std::vector x(dynamic_size), y(dynamic_size), z(dynamic_size), out(dynamic_size); + init(x, y, z, dynamic_size); + for (auto _ : state) + { + for (std::size_t i = 0; i < dynamic_size; ++i) + { + out[i] = (x[i] + y[i] * 0.5) * (z[i] - 0.25) + y[i] / (x[i] + 2.0); + } + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void broadcast_dynamic_xtensor(benchmark::State& state) + { + const std::vector shape = {64, 64, 16}; + xarray x, row, column, out; + x.resize(shape); row.resize({16}); column.resize({64, 1, 1}); out.resize(shape); + xarray unused; + unused.resize(shape); + init(x, unused, out, dynamic_size); + for (std::size_t i = 0; i < row.size(); ++i) row[i] = double(i) * 0.25; + for (std::size_t i = 0; i < column.size(); ++i) column[i] = double(i) * 0.125; + for (auto _ : state) + { + noalias(out) = (x + row) * 1.5 - column; + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void broadcast_dynamic_raw(benchmark::State& state) + { + std::vector x(dynamic_size), row(16), column(64), out(dynamic_size), unused(dynamic_size); + init(x, unused, out, dynamic_size); + for (std::size_t i = 0; i < row.size(); ++i) row[i] = double(i) * 0.25; + for (std::size_t i = 0; i < column.size(); ++i) column[i] = double(i) * 0.125; + for (auto _ : state) + { + for (std::size_t i = 0; i < 64; ++i) + { + for (std::size_t j = 0; j < 64; ++j) + { + for (std::size_t k = 0; k < 16; ++k) + { + const std::size_t n = (i * 64 + j) * 16 + k; + out[n] = (x[n] + row[k]) * 1.5 - column[i]; + } + } + } + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void transpose_cast_raw(benchmark::State& state) + { + constexpr std::size_t height = 128; + constexpr std::size_t width = 256; + constexpr std::size_t channels = 3; + std::vector input(height * width * channels, 127); + std::vector out(input.size()); + for (auto _ : state) + { + for (std::size_t c = 0; c < channels; ++c) + { + for (std::size_t i = 0; i < height; ++i) + { + for (std::size_t j = 0; j < width; ++j) + { + out[(c * height + i) * width + j] = float(input[(i * width + j) * channels + c]) / 255.0f; + } + } + } + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void transpose_cast_xtensor(benchmark::State& state) + { + xarray input = xarray::from_shape({128, 256, 3}); + xarray out = xarray::from_shape({3, 128, 256}); + std::fill(input.begin(), input.end(), 127); + for (auto _ : state) + { + noalias(out) = cast(transpose(input, {2, 0, 1})) / 255.0f; + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void linear_fixed_xtensor(benchmark::State& state) + { + xtensor_fixed> x, y, z, out; + init(x, y, z, fixed_size); + for (auto _ : state) + { + noalias(out) = (x + y * 0.5) * (z - 0.25) + y / (x + 2.0); + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + +#ifdef XTENSOR_BENCHMARK_USE_EIGEN + void linear_dynamic_eigen(benchmark::State& state) + { + Eigen::ArrayXd x(dynamic_size), y(dynamic_size), z(dynamic_size), out(dynamic_size); + init(x, y, z, dynamic_size); + for (auto _ : state) + { + out = (x + y * 0.5) * (z - 0.25) + y / (x + 2.0); + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void linear_fixed_eigen(benchmark::State& state) + { + Eigen::Array x, y, z, out; + init(x, y, z, fixed_size); + for (auto _ : state) + { + out = (x + y * 0.5) * (z - 0.25) + y / (x + 2.0); + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } +#endif + +#ifdef XTENSOR_BENCHMARK_USE_ARMADILLO + void linear_dynamic_armadillo(benchmark::State& state) + { + arma::cube x(64, 64, 16), y(64, 64, 16), z(64, 64, 16), out(64, 64, 16); + init(x, y, z, dynamic_size); + for (auto _ : state) + { + out = (x + y * 0.5) % (z - 0.25) + y / (x + 2.0); + benchmark::DoNotOptimize(out.memptr()); + benchmark::ClobberMemory(); + } + } + + void linear_fixed_armadillo(benchmark::State& state) + { + arma::vec::fixed x, y, z, out; + init(x, y, z, fixed_size); + for (auto _ : state) + { + out = (x + y * 0.5) % (z - 0.25) + y / (x + 2.0); + benchmark::DoNotOptimize(out.memptr()); + benchmark::ClobberMemory(); + } + } +#endif + + BENCHMARK(linear_dynamic_raw); + BENCHMARK(linear_dynamic_xtensor); + BENCHMARK(broadcast_dynamic_raw); + BENCHMARK(broadcast_dynamic_xtensor); + BENCHMARK(transpose_cast_raw); + BENCHMARK(transpose_cast_xtensor); + BENCHMARK(linear_fixed_xtensor); +#ifdef XTENSOR_BENCHMARK_USE_EIGEN + BENCHMARK(linear_dynamic_eigen); + BENCHMARK(linear_fixed_eigen); +#endif +#ifdef XTENSOR_BENCHMARK_USE_ARMADILLO + BENCHMARK(linear_dynamic_armadillo); + BENCHMARK(linear_fixed_armadillo); +#endif +} diff --git a/include/xtensor/core/xiterator.hpp b/include/xtensor/core/xiterator.hpp index 448f6093d..0345bb707 100644 --- a/include/xtensor/core/xiterator.hpp +++ b/include/xtensor/core/xiterator.hpp @@ -566,6 +566,13 @@ namespace xt inline auto xstepper::step_simd() -> simd_return_type { using simd_type = simd_return_type; + const auto& strides = p_c->strides(); + const bool broadcast = !strides.empty() + && (p_c->layout() == layout_type::column_major ? strides.front() : strides.back()) == 0; + if (broadcast) + { + return simd_type(*m_it); + } simd_type reg = detail::step_simd_invoker::template apply(m_it); m_it += xt_simd::revert_simd_traits::size; return reg; diff --git a/include/xtensor/core/xstrides.hpp b/include/xtensor/core/xstrides.hpp index 998ffcd6b..278a93649 100644 --- a/include/xtensor/core/xstrides.hpp +++ b/include/xtensor/core/xstrides.hpp @@ -810,7 +810,8 @@ namespace xt for (; s2_index != 0; --s1_index, --s2_index) { - if (static_cast(s1[s1_index - 1]) != static_cast(s2[s2_index - 1])) + if (static_cast(s2[s2_index - 1]) != value_type(0) + && static_cast(s1[s1_index - 1]) != static_cast(s2[s2_index - 1])) { break; } @@ -841,7 +842,8 @@ namespace xt for (; index < size; ++index) { - if (static_cast(s1[index]) != static_cast(s2[index])) + if (static_cast(s2[index]) != value_type(0) + && static_cast(s1[index]) != static_cast(s2[index])) { break; } diff --git a/test/test_strided_assign.cpp b/test/test_strided_assign.cpp index a0131b27d..d877fa3ea 100644 --- a/test/test_strided_assign.cpp +++ b/test/test_strided_assign.cpp @@ -196,4 +196,51 @@ namespace xt } } } + + TEST(xassign_strided, simd_broadcast) + { + xarray x = xarray::from_shape({4, 3, 8}); + xarray row = xarray::from_shape({8}); + xarray column = xarray::from_shape({4, 1, 1}); + xarray result = xarray::from_shape(x.shape()); + for (std::size_t i = 0; i < x.size(); ++i) x[i] = static_cast(i); + for (std::size_t i = 0; i < row.size(); ++i) row[i] = static_cast(i * 2); + for (std::size_t i = 0; i < column.size(); ++i) column[i] = static_cast(i * 3); + + noalias(result) = (x + row) * 1.5 - column; + + EXPECT_TRUE(strided_assign_detail::get_loop_sizes(result, (x + row) * 1.5 - column).can_do_strided_assign); + for (std::size_t i = 0; i < 4; ++i) + { + for (std::size_t j = 0; j < 3; ++j) + { + for (std::size_t k = 0; k < 8; ++k) + { + EXPECT_EQ(result(i, j, k), (x(i, j, k) + row(k)) * 1.5 - column(i, 0, 0)); + } + } + } + + xarray cx = xarray::from_shape({8, 3, 4}); + xarray leading = xarray::from_shape({8, 1, 1}); + xarray rest = xarray::from_shape({1, 3, 4}); + xarray cresult = xarray::from_shape(cx.shape()); + for (std::size_t i = 0; i < cx.size(); ++i) cx[i] = static_cast(i); + for (std::size_t i = 0; i < leading.size(); ++i) leading[i] = static_cast(i * 2); + for (std::size_t i = 0; i < rest.size(); ++i) rest[i] = static_cast(i * 3); + + noalias(cresult) = (cx + leading) * 1.5 - rest; + + EXPECT_TRUE(strided_assign_detail::get_loop_sizes(cresult, (cx + leading) * 1.5 - rest).can_do_strided_assign); + for (std::size_t i = 0; i < 8; ++i) + { + for (std::size_t j = 0; j < 3; ++j) + { + for (std::size_t k = 0; k < 4; ++k) + { + EXPECT_EQ(cresult(i, j, k), (cx(i, j, k) + leading(i, 0, 0)) * 1.5 - rest(0, j, k)); + } + } + } + } } From d45699748738bf0c4b0b9d58f8ecef4ad62fe32b Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Sun, 4 Oct 2026 19:24:52 +0200 Subject: [PATCH 2/8] Vectorize strided transpose assignments --- benchmark/benchmark_compare.cpp | 5 +-- include/xtensor/core/xassign.hpp | 48 ++++++++++++++++++++++--- include/xtensor/core/xiterator.hpp | 16 +++++++++ include/xtensor/core/xoperation.hpp | 13 ++++--- include/xtensor/views/xstrided_view.hpp | 13 +++++++ test/test_strided_assign.cpp | 42 ++++++++++++++++++---- 6 files changed, 119 insertions(+), 18 deletions(-) diff --git a/benchmark/benchmark_compare.cpp b/benchmark/benchmark_compare.cpp index 8b5655ccf..31bb5dc3b 100644 --- a/benchmark/benchmark_compare.cpp +++ b/benchmark/benchmark_compare.cpp @@ -119,8 +119,9 @@ namespace xt::compare constexpr std::size_t height = 128; constexpr std::size_t width = 256; constexpr std::size_t channels = 3; - std::vector input(height * width * channels, 127); + std::vector input(height * width * channels); std::vector out(input.size()); + for (std::size_t i = 0; i < input.size(); ++i) input[i] = static_cast(i % 251); for (auto _ : state) { for (std::size_t c = 0; c < channels; ++c) @@ -142,7 +143,7 @@ namespace xt::compare { xarray input = xarray::from_shape({128, 256, 3}); xarray out = xarray::from_shape({3, 128, 256}); - std::fill(input.begin(), input.end(), 127); + for (std::size_t i = 0; i < input.size(); ++i) input.storage()[i] = static_cast(i % 251); for (auto _ : state) { noalias(out) = cast(transpose(input, {2, 0, 1})) / 255.0f; diff --git a/include/xtensor/core/xassign.hpp b/include/xtensor/core/xassign.hpp index 174ddcdc0..653a93227 100644 --- a/include/xtensor/core/xassign.hpp +++ b/include/xtensor/core/xassign.hpp @@ -320,6 +320,31 @@ namespace xt static constexpr bool value = std::conjunction>...>::value; }; + template + struct has_strided_simd_interface : has_simd_interface + { + }; + + template + struct has_strided_simd_interface, T> + : std::conjunction< + has_simd_type, + has_simd_apply>, + has_strided_simd_interface, T>...> + { + }; + + template + struct supports_row_major_gather : std::bool_constant + { + }; + + template + struct supports_row_major_gather> + : std::conjunction>...> + { + }; + /** * Considering the assignment LHS = RHS, if the requested value type used for * loading simd from RHS is not complex while LHS value_type is complex, @@ -423,7 +448,9 @@ namespace xt static constexpr bool simd_strided_assign() { - return strided_assign() && simd_assign(); + return strided_assign() && convertible_types() && simd_size() + && detail::has_strided_simd_interface::value + && detail::has_strided_simd_interface::value; } static constexpr bool simd_linear_assign(const E1& e1, const E2& e2) @@ -1045,8 +1072,19 @@ namespace xt if (is_row_major) { - auto csf = check_strides_functor(e1.strides()); - cut = csf(e2); + if constexpr ( + possible && E1::contiguous_layout && E1::static_layout == layout_type::row_major + && xt::detail::supports_row_major_gather::value) + { + // SIMD steppers can broadcast and gather arbitrary row-major input strides. + cut = e1.strides().size() - 1; + } + else + { + auto csf = check_strides_functor(e1.strides()); + cut = csf(e2); + } + // This makes that only one dimension will be treated in the inner loop. if (cut < e1.strides().size() - 1) { @@ -1054,7 +1092,7 @@ namespace xt cut = e1.strides().size() - 1; } } - else if (!is_row_major) + else { auto csf = check_strides_functor(e1.strides() ); @@ -1064,7 +1102,7 @@ namespace xt // Only make the inner loop go over one dimension by default for now cut = 1; } - } // can't reach here because this would have already triggered the fallback + } std::size_t outer_loop_size = static_cast(std::accumulate( e1.shape().begin(), diff --git a/include/xtensor/core/xiterator.hpp b/include/xtensor/core/xiterator.hpp index 0345bb707..547235e01 100644 --- a/include/xtensor/core/xiterator.hpp +++ b/include/xtensor/core/xiterator.hpp @@ -573,6 +573,22 @@ namespace xt { return simd_type(*m_it); } +#ifdef XTENSOR_USE_XSIMD + if constexpr (!xt_simd::is_batch_bool::value && !xt_simd::is_batch_complex::value) + { + const auto stride = p_c->layout() == layout_type::column_major ? strides.front() : strides.back(); + if (stride != 1) + { + using index_value_type = xsimd::as_integer_t; + using arch_type = typename simd_type::arch_type; + xsimd::batch index = xsimd::make_iota_batch_constant(); + index *= static_cast(stride); + simd_type reg = simd_type::gather(&(*m_it), index); + m_it += static_cast(xt_simd::revert_simd_traits::size * stride); + return reg; + } + } +#endif simd_type reg = detail::step_simd_invoker::template apply(m_it); m_it += xt_simd::revert_simd_traits::size; return reg; diff --git a/include/xtensor/core/xoperation.hpp b/include/xtensor/core/xoperation.hpp index 76bfc7c50..5bd28dfe5 100644 --- a/include/xtensor/core/xoperation.hpp +++ b/include/xtensor/core/xoperation.hpp @@ -155,12 +155,15 @@ namespace xt return static_cast(arg); } - // SIMD conversion disabled for now since it does not make sense - // in most of the cases - /*constexpr simd_result_type simd_apply(const simd_value_type& arg) const + template < + class B, + std::enable_if_t< + std::is_same::type>::value, + int> = 0> + constexpr B simd_apply(const B& arg) const { - return static_cast(arg); - }*/ + return arg; + } }; }; diff --git a/include/xtensor/views/xstrided_view.hpp b/include/xtensor/views/xstrided_view.hpp index b0fe2c6e0..586a67133 100644 --- a/include/xtensor/views/xstrided_view.hpp +++ b/include/xtensor/views/xstrided_view.hpp @@ -313,6 +313,19 @@ namespace xt friend class xaxis_slice_iterator; }; + namespace detail + { + template + struct has_strided_simd_interface, T> + : std::conjunction< + has_simd_type, + std::bool_constant::value_type, T>::value>, + std::is_lvalue_reference::const_reference>, + std::bool_constant::value>> + { + }; + } + /************************** * xstrided_view builders * **************************/ diff --git a/test/test_strided_assign.cpp b/test/test_strided_assign.cpp index d877fa3ea..b583e7b99 100644 --- a/test/test_strided_assign.cpp +++ b/test/test_strided_assign.cpp @@ -16,6 +16,8 @@ #include "xtensor/core/xassign.hpp" #include "xtensor/core/xlayout.hpp" #include "xtensor/core/xnoalias.hpp" +#include "xtensor/core/xoperation.hpp" +#include "xtensor/views/xstrided_view.hpp" #include "xtensor/views/xview.hpp" #include "test_common.hpp" @@ -197,15 +199,42 @@ namespace xt } } +#ifdef XTENSOR_USE_XSIMD + TEST(xassign_strided, simd_transpose_cast) + { + xarray input = xarray::from_shape({4, 16, 3}); + xarray result = xarray::from_shape({3, 4, 16}); + for (std::size_t i = 0; i < input.size(); ++i) + { + input.storage()[i] = static_cast(i % 251); + } + + auto expression = cast(transpose(input, {2, 0, 1})) / 255.0f; + noalias(result) = expression; + + EXPECT_TRUE((xassign_traits::simd_strided_assign())); + EXPECT_TRUE(strided_assign_detail::get_loop_sizes(result, expression).can_do_strided_assign); + for (std::size_t c = 0; c < 3; ++c) + { + for (std::size_t i = 0; i < 4; ++i) + { + for (std::size_t j = 0; j < 16; ++j) + { + EXPECT_EQ(result(c, i, j), static_cast(input(i, j, c)) / 255.0f); + } + } + } + } + TEST(xassign_strided, simd_broadcast) { xarray x = xarray::from_shape({4, 3, 8}); xarray row = xarray::from_shape({8}); xarray column = xarray::from_shape({4, 1, 1}); xarray result = xarray::from_shape(x.shape()); - for (std::size_t i = 0; i < x.size(); ++i) x[i] = static_cast(i); - for (std::size_t i = 0; i < row.size(); ++i) row[i] = static_cast(i * 2); - for (std::size_t i = 0; i < column.size(); ++i) column[i] = static_cast(i * 3); + for (std::size_t i = 0; i < x.size(); ++i) x.storage()[i] = static_cast(i); + for (std::size_t i = 0; i < row.size(); ++i) row.storage()[i] = static_cast(i * 2); + for (std::size_t i = 0; i < column.size(); ++i) column.storage()[i] = static_cast(i * 3); noalias(result) = (x + row) * 1.5 - column; @@ -225,9 +254,9 @@ namespace xt xarray leading = xarray::from_shape({8, 1, 1}); xarray rest = xarray::from_shape({1, 3, 4}); xarray cresult = xarray::from_shape(cx.shape()); - for (std::size_t i = 0; i < cx.size(); ++i) cx[i] = static_cast(i); - for (std::size_t i = 0; i < leading.size(); ++i) leading[i] = static_cast(i * 2); - for (std::size_t i = 0; i < rest.size(); ++i) rest[i] = static_cast(i * 3); + for (std::size_t i = 0; i < cx.size(); ++i) cx.storage()[i] = static_cast(i); + for (std::size_t i = 0; i < leading.size(); ++i) leading.storage()[i] = static_cast(i * 2); + for (std::size_t i = 0; i < rest.size(); ++i) rest.storage()[i] = static_cast(i * 3); noalias(cresult) = (cx + leading) * 1.5 - rest; @@ -243,4 +272,5 @@ namespace xt } } } +#endif } From f15367ddef0f550b774e718395ae9672efbedc39 Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Sun, 4 Oct 2026 21:30:01 +0200 Subject: [PATCH 3/8] Plan broadcast and permutation assignments --- include/xtensor/core/xassign.hpp | 194 ++++++++++++++++++++++++++++ include/xtensor/core/xoperation.hpp | 138 ++++++++++++++++++++ test/test_strided_assign.cpp | 49 ++++++- 3 files changed, 379 insertions(+), 2 deletions(-) diff --git a/include/xtensor/core/xassign.hpp b/include/xtensor/core/xassign.hpp index 653a93227..153484a93 100644 --- a/include/xtensor/core/xassign.hpp +++ b/include/xtensor/core/xassign.hpp @@ -11,7 +11,9 @@ #define XTENSOR_ASSIGN_HPP #include +#include #include +#include #include #include @@ -58,6 +60,12 @@ namespace xt template void strided_assign(E1& e1, const E2& e2, std::true_type /*enable*/); + template + bool try_permutation_assign(E1&, const E2&) + { + return false; + } + /************************ * xexpression_assigner * ************************/ @@ -334,6 +342,187 @@ namespace xt { }; + template + struct runtime_plan_expression : std::false_type + { + }; + + template + struct runtime_plan_expression< + E, + std::void_t().storage().data())>> + : std::bool_constant< + data_interface_expression && std::is_arithmetic::value + && std::is_pointer().storage().data())>::value + && std::is_same< + std::remove_cv_t().storage().data())>>, + typename E::value_type>::value> + { + }; + + template + struct runtime_plan_expression, void> : std::bool_constant::value> + { + }; + + template + struct runtime_plan_expression, void> + : std::conjunction>...> + { + }; + + template + struct runtime_plan_output : std::false_type + { + }; + + template + struct runtime_plan_output< + E, + std::void_t< + decltype(std::declval().storage().data()), + decltype(std::declval().data_offset())>> + : std::bool_constant< + E::contiguous_layout && E::static_layout == layout_type::row_major + && std::is_pointer().storage().data())>::value + && std::is_same< + std::remove_cv_t().storage().data())>>, + typename E::value_type>::value> + { + }; + + template + struct runtime_leaf_plan + { + using value_type = typename E::value_type; + const value_type* data; + std::array strides; + + value_type operator()(std::size_t i, std::size_t j, std::size_t k) const + { + const auto offset = static_cast(i) * strides[0] + + static_cast(j) * strides[1]; + if (strides[2] == 0) + { + return data[offset]; + } + if (strides[2] == 1) + { + return data[offset + static_cast(k)]; + } + return data[offset + static_cast(k) * strides[2]]; + } + }; + + template + struct runtime_scalar_plan + { + T value; + + T operator()(std::size_t, std::size_t, std::size_t) const + { + return value; + } + }; + + template + struct runtime_function_plan + { + F functor; + P arguments; + + auto operator()(std::size_t i, std::size_t j, std::size_t k) const + { + return std::apply( + [&](const auto&... argument) { return functor(argument(i, j, k)...); }, + arguments + ); + } + }; + + template + struct runtime_plan_builder + { + static auto make(const E& expression) + { + std::array strides = {0, 0, 0}; + const std::size_t offset = 3 - expression.dimension(); + for (std::size_t i = 0; i < expression.dimension(); ++i) + { + if (expression.shape()[i] != 1) + { + strides[offset + i] = expression.strides()[i]; + } + } + return runtime_leaf_plan{expression.storage().data() + expression.data_offset(), strides}; + } + }; + + template + struct runtime_plan_builder> + { + static auto make(const xscalar& expression) + { + return runtime_scalar_plan{expression()}; + } + }; + + template + struct runtime_plan_builder> + { + using expression_type = xfunction; + + template + static auto make(const expression_type& expression, std::index_sequence) + { + auto arguments = std::make_tuple( + runtime_plan_builder(expression.arguments()))>>::make( + std::get(expression.arguments()) + )... + ); + return runtime_function_plan, decltype(arguments)>{ + expression.functor(), + std::move(arguments) + }; + } + + static auto make(const expression_type& expression) + { + return make(expression, std::index_sequence_for{}); + } + }; + + template + bool try_runtime_plan(E1& dst, const E2& expression, bool trivial) + { + if (trivial) + { + return false; + } + if constexpr ( + runtime_plan_output::value && runtime_plan_expression::value) + { + if (dst.dimension() == 3 && expression.dimension() <= 3) + { + auto plan = runtime_plan_builder::make(expression); + auto* output = dst.storage().data() + dst.data_offset(); + std::size_t n = 0; + for (std::size_t i = 0; i < dst.shape()[0]; ++i) + { + for (std::size_t j = 0; j < dst.shape()[1]; ++j) + { + for (std::size_t k = 0; k < dst.shape()[2]; ++k) + { + output[n++] = plan(i, j, k); + } + } + } + return true; + } + } + return false; + } + template struct supports_row_major_gather : std::bool_constant { @@ -474,6 +663,11 @@ namespace xt const E2& de2 = e2.derived_cast(); using traits = xassign_traits; + if (try_permutation_assign(de1, de2) || detail::try_runtime_plan(de1, de2, trivial)) + { + return; + } + bool linear_assign = traits::linear_assign(de1, de2, trivial); constexpr bool simd_assign = traits::simd_assign(); constexpr bool simd_linear_assign = traits::simd_linear_assign(); diff --git a/include/xtensor/core/xoperation.hpp b/include/xtensor/core/xoperation.hpp index 5bd28dfe5..5fb6ac3ae 100644 --- a/include/xtensor/core/xoperation.hpp +++ b/include/xtensor/core/xoperation.hpp @@ -11,6 +11,7 @@ #define XTENSOR_OPERATION_HPP #include +#include #include #include @@ -148,6 +149,7 @@ namespace xt struct functor { using result_type = R; + using cast_result_type = R; template constexpr result_type operator()(const A1& arg) const @@ -167,6 +169,117 @@ namespace xt }; }; + template + struct cast_expression : std::false_type + { + }; + + template + struct cast_expression, std::void_t> : std::true_type + { + using result_type = typename F::cast_result_type; + }; + + template + struct permutation_view : std::false_type + { + }; + + template + struct permutation_view> + : std::bool_constant< + L == layout_type::dynamic && is_specialization_of::value> + { + }; + + template + inline void run_permutation_assign(E1& dst, const V& src, S scale) + { + const std::size_t shape_a = src.shape()[A]; + const std::size_t shape_b = src.shape()[B]; + const std::size_t shape_c = src.shape()[C]; + const std::size_t src_stride_a = static_cast(src.strides()[A]); + const std::size_t src_stride_b = static_cast(src.strides()[B]); + const std::size_t src_stride_c = static_cast(src.strides()[C]); + const std::size_t dst_stride_a = static_cast(dst.strides()[A]); + const std::size_t dst_stride_b = static_cast(dst.strides()[B]); + const std::size_t dst_stride_c = static_cast(dst.strides()[C]); + auto* src_data = src.data() + src.data_offset(); + auto* dst_data = dst.data() + dst.data_offset(); + if ( + shape_c == 3 && src_stride_c == 1 && src_stride_b == 3 && src_stride_a == shape_b * 3 + && dst_stride_b == 1 && dst_stride_a == shape_b && dst_stride_c == shape_a * shape_b) + { + const std::size_t plane_size = shape_a * shape_b; + auto* dst_0 = dst_data; + auto* dst_1 = dst_data + plane_size; + auto* dst_2 = dst_data + 2 * plane_size; + for (std::size_t i = 0; i < plane_size; ++i) + { + dst_0[i] = static_cast(src_data[3 * i]) / static_cast(scale); + dst_1[i] = static_cast(src_data[3 * i + 1]) / static_cast(scale); + dst_2[i] = static_cast(src_data[3 * i + 2]) / static_cast(scale); + } + return; + } + constexpr std::size_t tile_size = 32; + for (std::size_t a0 = 0; a0 < shape_a; a0 += tile_size) + { + for (std::size_t b0 = 0; b0 < shape_b; b0 += tile_size) + { + const std::size_t a_end = std::min(a0 + tile_size, shape_a); + const std::size_t b_end = std::min(b0 + tile_size, shape_b); + for (std::size_t a = a0; a < a_end; ++a) + { + for (std::size_t b = b0; b < b_end; ++b) + { + for (std::size_t c = 0; c < shape_c; ++c) + { + const auto src_offset = a * src_stride_a + b * src_stride_b + c * src_stride_c; + const auto dst_offset = a * dst_stride_a + b * dst_stride_b + c * dst_stride_c; + dst_data[dst_offset] = static_cast(src_data[src_offset]) / static_cast(scale); + } + } + } + } + } + } + + template + inline bool try_3d_permutation_assign(E1& dst, const V& src, S scale) + { + if (dst.dimension() != 3 || src.dimension() != 3 || dst.shape() != src.shape() + || !dst.is_contiguous() || dst.layout() != layout_type::row_major) + { + return false; + } + + std::array axes = {0, 1, 2}; + std::sort(axes.begin(), axes.end(), [&](auto lhs, auto rhs) { return src.strides()[lhs] > src.strides()[rhs]; }); + if (src.strides()[axes[2]] != 1 + || static_cast(src.strides()[axes[1]]) != src.shape()[axes[2]] + || static_cast(src.strides()[axes[0]]) + != src.shape()[axes[2]] * src.shape()[axes[1]]) + { + return false; + } + +#define XTENSOR_RUN_PERMUTATION(A, B, C) \ + if (axes == std::array{A, B, C}) \ + { \ + run_permutation_assign(dst, src, scale); \ + return true; \ + } + XTENSOR_RUN_PERMUTATION(0, 1, 2) + XTENSOR_RUN_PERMUTATION(0, 2, 1) + XTENSOR_RUN_PERMUTATION(1, 0, 2) + XTENSOR_RUN_PERMUTATION(1, 2, 0) + XTENSOR_RUN_PERMUTATION(2, 0, 1) + XTENSOR_RUN_PERMUTATION(2, 1, 0) +#undef XTENSOR_RUN_PERMUTATION + return false; + } + template struct select_xfunction_expression; @@ -211,6 +324,31 @@ namespace xt enable_if_t...>::value, xfunction_type>::type; } + template + inline bool try_permutation_assign(E1& dst, const xfunction& expression) + { + using cast_expression_type = std::decay_t(expression.arguments()))>; + if constexpr (detail::cast_expression::value && is_xscalar>::value) + { + const auto& cast_expression = std::get<0>(expression.arguments()); + const auto& view = std::get<0>(cast_expression.arguments()); + using view_type = std::decay_t; + using cast_result_type = typename detail::cast_expression::result_type; + if constexpr ( + detail::permutation_view::value + && std::is_same::value + && std::is_floating_point::value) + { + return detail::try_3d_permutation_assign( + dst, + view, + std::get<1>(expression.arguments())() + ); + } + } + return false; + } + #undef UNARY_OPERATOR_FUNCTOR #undef BINARY_OPERATOR_FUNCTOR diff --git a/test/test_strided_assign.cpp b/test/test_strided_assign.cpp index b583e7b99..67a29a5ba 100644 --- a/test/test_strided_assign.cpp +++ b/test/test_strided_assign.cpp @@ -199,8 +199,7 @@ namespace xt } } -#ifdef XTENSOR_USE_XSIMD - TEST(xassign_strided, simd_transpose_cast) + TEST(xassign_strided, transpose_cast_permutation_plan) { xarray input = xarray::from_shape({4, 16, 3}); xarray result = xarray::from_shape({3, 4, 16}); @@ -212,8 +211,11 @@ namespace xt auto expression = cast(transpose(input, {2, 0, 1})) / 255.0f; noalias(result) = expression; + EXPECT_TRUE(try_permutation_assign(result, expression)); +#ifdef XTENSOR_USE_XSIMD EXPECT_TRUE((xassign_traits::simd_strided_assign())); EXPECT_TRUE(strided_assign_detail::get_loop_sizes(result, expression).can_do_strided_assign); +#endif for (std::size_t c = 0; c < 3; ++c) { for (std::size_t i = 0; i < 4; ++i) @@ -224,8 +226,28 @@ namespace xt } } } + + xarray tiled_input = xarray::from_shape({4, 5, 2}); + xarray tiled_result = xarray::from_shape({5, 4, 2}); + for (std::size_t i = 0; i < tiled_input.size(); ++i) + { + tiled_input.storage()[i] = static_cast(i); + } + auto tiled_expression = cast(transpose(tiled_input, {1, 0, 2})) / 7.0f; + EXPECT_TRUE(try_permutation_assign(tiled_result, tiled_expression)); + for (std::size_t i = 0; i < 5; ++i) + { + for (std::size_t j = 0; j < 4; ++j) + { + for (std::size_t c = 0; c < 2; ++c) + { + EXPECT_EQ(tiled_result(i, j, c), static_cast(tiled_input(j, i, c)) / 7.0f); + } + } + } } +#ifdef XTENSOR_USE_XSIMD TEST(xassign_strided, simd_broadcast) { xarray x = xarray::from_shape({4, 3, 8}); @@ -273,4 +295,27 @@ namespace xt } } #endif + + TEST(xassign_strided, runtime_broadcast_plan) + { + xarray x = xarray::from_shape({3, 4, 5}); + xarray row = {1.0, 2.0, 3.0, 4.0, 5.0}; + xarray column = xarray::from_shape({3, 1, 1}); + xarray result = xarray::from_shape(x.shape()); + for (std::size_t i = 0; i < x.size(); ++i) x.storage()[i] = static_cast(i); + for (std::size_t i = 0; i < column.size(); ++i) column.storage()[i] = static_cast(i * 3); + auto expression = (x + row) * 1.5 - column; + + EXPECT_TRUE(detail::try_runtime_plan(result, expression, false)); + for (std::size_t i = 0; i < 3; ++i) + { + for (std::size_t j = 0; j < 4; ++j) + { + for (std::size_t k = 0; k < 5; ++k) + { + EXPECT_EQ(result(i, j, k), (x(i, j, k) + row(k)) * 1.5 - column(i, 0, 0)); + } + } + } + } } From 67b993862e666475dabc3aadd664a10f70b10941 Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Sun, 4 Oct 2026 23:30:18 +0200 Subject: [PATCH 4/8] Benchmark Eigen tensor expressions --- benchmark/benchmark_compare.cpp | 38 +++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/benchmark/benchmark_compare.cpp b/benchmark/benchmark_compare.cpp index 31bb5dc3b..b004d58fb 100644 --- a/benchmark/benchmark_compare.cpp +++ b/benchmark/benchmark_compare.cpp @@ -20,6 +20,7 @@ #ifdef XTENSOR_BENCHMARK_USE_EIGEN #include +#include #endif #ifdef XTENSOR_BENCHMARK_USE_ARMADILLO @@ -188,6 +189,41 @@ namespace xt::compare benchmark::ClobberMemory(); } } + + void broadcast_dynamic_eigen(benchmark::State& state) + { + using tensor = Eigen::Tensor; + tensor x(64, 64, 16), out(64, 64, 16); + Eigen::Tensor row(16), column(64); + for (Eigen::Index i = 0; i < x.size(); ++i) x.data()[i] = 0.5 + double(i % 251) * 0.01; + for (Eigen::Index i = 0; i < row.size(); ++i) row(i) = double(i) * 0.25; + for (Eigen::Index i = 0; i < column.size(); ++i) column(i) = double(i) * 0.125; + const Eigen::array row_shape = {1, 1, 16}; + const Eigen::array row_broadcast = {64, 64, 1}; + const Eigen::array column_shape = {64, 1, 1}; + const Eigen::array column_broadcast = {1, 64, 16}; + for (auto _ : state) + { + out = (x + row.reshape(row_shape).broadcast(row_broadcast)) * 1.5 + - column.reshape(column_shape).broadcast(column_broadcast); + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } + + void transpose_cast_eigen(benchmark::State& state) + { + Eigen::Tensor input(128, 256, 3); + Eigen::Tensor out(3, 128, 256); + for (Eigen::Index i = 0; i < input.size(); ++i) input.data()[i] = static_cast(i % 251); + const Eigen::array permutation = {2, 0, 1}; + for (auto _ : state) + { + out = input.shuffle(permutation).cast() / 255.0f; + benchmark::DoNotOptimize(out.data()); + benchmark::ClobberMemory(); + } + } #endif #ifdef XTENSOR_BENCHMARK_USE_ARMADILLO @@ -226,6 +262,8 @@ namespace xt::compare #ifdef XTENSOR_BENCHMARK_USE_EIGEN BENCHMARK(linear_dynamic_eigen); BENCHMARK(linear_fixed_eigen); + BENCHMARK(broadcast_dynamic_eigen); + BENCHMARK(transpose_cast_eigen); #endif #ifdef XTENSOR_BENCHMARK_USE_ARMADILLO BENCHMARK(linear_dynamic_armadillo); From 49f9028bafe9c27d551836f90cb32647de389992 Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Mon, 5 Oct 2026 12:33:48 +0200 Subject: [PATCH 5/8] Fix CI formatting and bool trait detection --- benchmark/benchmark_compare.cpp | 63 +++++++++++++++----- include/xtensor/core/xassign.hpp | 34 +++++------ include/xtensor/core/xiterator.hpp | 6 +- include/xtensor/core/xoperation.hpp | 35 ++++++------ test/test_strided_assign.cpp | 89 ++++++++++++++++++++++------- 5 files changed, 156 insertions(+), 71 deletions(-) diff --git a/benchmark/benchmark_compare.cpp b/benchmark/benchmark_compare.cpp index b004d58fb..a73902685 100644 --- a/benchmark/benchmark_compare.cpp +++ b/benchmark/benchmark_compare.cpp @@ -48,7 +48,10 @@ namespace xt::compare { const std::vector shape = {64, 64, 16}; xarray x, y, z, out; - x.resize(shape); y.resize(shape); z.resize(shape); out.resize(shape); + x.resize(shape); + y.resize(shape); + z.resize(shape); + out.resize(shape); init(x, y, z, dynamic_size); for (auto _ : state) { @@ -77,12 +80,21 @@ namespace xt::compare { const std::vector shape = {64, 64, 16}; xarray x, row, column, out; - x.resize(shape); row.resize({16}); column.resize({64, 1, 1}); out.resize(shape); + x.resize(shape); + row.resize({16}); + column.resize({64, 1, 1}); + out.resize(shape); xarray unused; unused.resize(shape); init(x, unused, out, dynamic_size); - for (std::size_t i = 0; i < row.size(); ++i) row[i] = double(i) * 0.25; - for (std::size_t i = 0; i < column.size(); ++i) column[i] = double(i) * 0.125; + for (std::size_t i = 0; i < row.size(); ++i) + { + row[i] = double(i) * 0.25; + } + for (std::size_t i = 0; i < column.size(); ++i) + { + column[i] = double(i) * 0.125; + } for (auto _ : state) { noalias(out) = (x + row) * 1.5 - column; @@ -95,8 +107,14 @@ namespace xt::compare { std::vector x(dynamic_size), row(16), column(64), out(dynamic_size), unused(dynamic_size); init(x, unused, out, dynamic_size); - for (std::size_t i = 0; i < row.size(); ++i) row[i] = double(i) * 0.25; - for (std::size_t i = 0; i < column.size(); ++i) column[i] = double(i) * 0.125; + for (std::size_t i = 0; i < row.size(); ++i) + { + row[i] = double(i) * 0.25; + } + for (std::size_t i = 0; i < column.size(); ++i) + { + column[i] = double(i) * 0.125; + } for (auto _ : state) { for (std::size_t i = 0; i < 64; ++i) @@ -122,7 +140,10 @@ namespace xt::compare constexpr std::size_t channels = 3; std::vector input(height * width * channels); std::vector out(input.size()); - for (std::size_t i = 0; i < input.size(); ++i) input[i] = static_cast(i % 251); + for (std::size_t i = 0; i < input.size(); ++i) + { + input[i] = static_cast(i % 251); + } for (auto _ : state) { for (std::size_t c = 0; c < channels; ++c) @@ -131,7 +152,8 @@ namespace xt::compare { for (std::size_t j = 0; j < width; ++j) { - out[(c * height + i) * width + j] = float(input[(i * width + j) * channels + c]) / 255.0f; + out[(c * height + i) * width + j] = float(input[(i * width + j) * channels + c]) + / 255.0f; } } } @@ -144,7 +166,10 @@ namespace xt::compare { xarray input = xarray::from_shape({128, 256, 3}); xarray out = xarray::from_shape({3, 128, 256}); - for (std::size_t i = 0; i < input.size(); ++i) input.storage()[i] = static_cast(i % 251); + for (std::size_t i = 0; i < input.size(); ++i) + { + input.storage()[i] = static_cast(i % 251); + } for (auto _ : state) { noalias(out) = cast(transpose(input, {2, 0, 1})) / 255.0f; @@ -195,9 +220,18 @@ namespace xt::compare using tensor = Eigen::Tensor; tensor x(64, 64, 16), out(64, 64, 16); Eigen::Tensor row(16), column(64); - for (Eigen::Index i = 0; i < x.size(); ++i) x.data()[i] = 0.5 + double(i % 251) * 0.01; - for (Eigen::Index i = 0; i < row.size(); ++i) row(i) = double(i) * 0.25; - for (Eigen::Index i = 0; i < column.size(); ++i) column(i) = double(i) * 0.125; + for (Eigen::Index i = 0; i < x.size(); ++i) + { + x.data()[i] = 0.5 + double(i % 251) * 0.01; + } + for (Eigen::Index i = 0; i < row.size(); ++i) + { + row(i) = double(i) * 0.25; + } + for (Eigen::Index i = 0; i < column.size(); ++i) + { + column(i) = double(i) * 0.125; + } const Eigen::array row_shape = {1, 1, 16}; const Eigen::array row_broadcast = {64, 64, 1}; const Eigen::array column_shape = {64, 1, 1}; @@ -215,7 +249,10 @@ namespace xt::compare { Eigen::Tensor input(128, 256, 3); Eigen::Tensor out(3, 128, 256); - for (Eigen::Index i = 0; i < input.size(); ++i) input.data()[i] = static_cast(i % 251); + for (Eigen::Index i = 0; i < input.size(); ++i) + { + input.data()[i] = static_cast(i % 251); + } const Eigen::array permutation = {2, 0, 1}; for (auto _ : state) { diff --git a/include/xtensor/core/xassign.hpp b/include/xtensor/core/xassign.hpp index 153484a93..9c3ebde0f 100644 --- a/include/xtensor/core/xassign.hpp +++ b/include/xtensor/core/xassign.hpp @@ -342,15 +342,13 @@ namespace xt { }; - template + template ::value, class = void> struct runtime_plan_expression : std::false_type { }; template - struct runtime_plan_expression< - E, - std::void_t().storage().data())>> + struct runtime_plan_expression().storage().data())>> : std::bool_constant< data_interface_expression && std::is_arithmetic::value && std::is_pointer().storage().data())>::value @@ -361,17 +359,18 @@ namespace xt }; template - struct runtime_plan_expression, void> : std::bool_constant::value> + struct runtime_plan_expression, true, void> + : std::bool_constant::value> { }; template - struct runtime_plan_expression, void> + struct runtime_plan_expression, true, void> : std::conjunction>...> { }; - template + template ::value, class = void> struct runtime_plan_output : std::false_type { }; @@ -379,9 +378,8 @@ namespace xt template struct runtime_plan_output< E, - std::void_t< - decltype(std::declval().storage().data()), - decltype(std::declval().data_offset())>> + true, + std::void_t().storage().data()), decltype(std::declval().data_offset())>> : std::bool_constant< E::contiguous_layout && E::static_layout == layout_type::row_major && std::is_pointer().storage().data())>::value @@ -434,7 +432,10 @@ namespace xt auto operator()(std::size_t i, std::size_t j, std::size_t k) const { return std::apply( - [&](const auto&... argument) { return functor(argument(i, j, k)...); }, + [&](const auto&... argument) + { + return functor(argument(i, j, k)...); + }, arguments ); } @@ -499,8 +500,7 @@ namespace xt { return false; } - if constexpr ( - runtime_plan_output::value && runtime_plan_expression::value) + if constexpr (runtime_plan_output::value && runtime_plan_expression::value) { if (dst.dimension() == 3 && expression.dimension() <= 3) { @@ -1266,16 +1266,16 @@ namespace xt if (is_row_major) { - if constexpr ( - possible && E1::contiguous_layout && E1::static_layout == layout_type::row_major - && xt::detail::supports_row_major_gather::value) + if constexpr (possible && E1::contiguous_layout && E1::static_layout == layout_type::row_major && xt::detail::supports_row_major_gather::value) { // SIMD steppers can broadcast and gather arbitrary row-major input strides. cut = e1.strides().size() - 1; } else { - auto csf = check_strides_functor(e1.strides()); + auto csf = check_strides_functor( + e1.strides() + ); cut = csf(e2); } diff --git a/include/xtensor/core/xiterator.hpp b/include/xtensor/core/xiterator.hpp index 547235e01..2805afb14 100644 --- a/include/xtensor/core/xiterator.hpp +++ b/include/xtensor/core/xiterator.hpp @@ -568,7 +568,8 @@ namespace xt using simd_type = simd_return_type; const auto& strides = p_c->strides(); const bool broadcast = !strides.empty() - && (p_c->layout() == layout_type::column_major ? strides.front() : strides.back()) == 0; + && (p_c->layout() == layout_type::column_major ? strides.front() : strides.back() + ) == 0; if (broadcast) { return simd_type(*m_it); @@ -581,7 +582,8 @@ namespace xt { using index_value_type = xsimd::as_integer_t; using arch_type = typename simd_type::arch_type; - xsimd::batch index = xsimd::make_iota_batch_constant(); + xsimd::batch + index = xsimd::make_iota_batch_constant(); index *= static_cast(stride); simd_type reg = simd_type::gather(&(*m_it), index); m_it += static_cast(xt_simd::revert_simd_traits::size * stride); diff --git a/include/xtensor/core/xoperation.hpp b/include/xtensor/core/xoperation.hpp index 5fb6ac3ae..74e77bbc4 100644 --- a/include/xtensor/core/xoperation.hpp +++ b/include/xtensor/core/xoperation.hpp @@ -159,9 +159,7 @@ namespace xt template < class B, - std::enable_if_t< - std::is_same::type>::value, - int> = 0> + std::enable_if_t::type>::value, int> = 0> constexpr B simd_apply(const B& arg) const { return arg; @@ -187,8 +185,7 @@ namespace xt template struct permutation_view> - : std::bool_constant< - L == layout_type::dynamic && is_specialization_of::value> + : std::bool_constant::value> { }; @@ -206,8 +203,7 @@ namespace xt const std::size_t dst_stride_c = static_cast(dst.strides()[C]); auto* src_data = src.data() + src.data_offset(); auto* dst_data = dst.data() + dst.data_offset(); - if ( - shape_c == 3 && src_stride_c == 1 && src_stride_b == 3 && src_stride_a == shape_b * 3 + if (shape_c == 3 && src_stride_c == 1 && src_stride_b == 3 && src_stride_a == shape_b * 3 && dst_stride_b == 1 && dst_stride_a == shape_b && dst_stride_c == shape_a * shape_b) { const std::size_t plane_size = shape_a * shape_b; @@ -237,7 +233,8 @@ namespace xt { const auto src_offset = a * src_stride_a + b * src_stride_b + c * src_stride_c; const auto dst_offset = a * dst_stride_a + b * dst_stride_b + c * dst_stride_c; - dst_data[dst_offset] = static_cast(src_data[src_offset]) / static_cast(scale); + dst_data[dst_offset] = static_cast(src_data[src_offset]) + / static_cast(scale); } } } @@ -255,7 +252,14 @@ namespace xt } std::array axes = {0, 1, 2}; - std::sort(axes.begin(), axes.end(), [&](auto lhs, auto rhs) { return src.strides()[lhs] > src.strides()[rhs]; }); + std::sort( + axes.begin(), + axes.end(), + [&](auto lhs, auto rhs) + { + return src.strides()[lhs] > src.strides()[rhs]; + } + ); if (src.strides()[axes[2]] != 1 || static_cast(src.strides()[axes[1]]) != src.shape()[axes[2]] || static_cast(src.strides()[axes[0]]) @@ -264,11 +268,11 @@ namespace xt return false; } -#define XTENSOR_RUN_PERMUTATION(A, B, C) \ - if (axes == std::array{A, B, C}) \ - { \ +#define XTENSOR_RUN_PERMUTATION(A, B, C) \ + if (axes == std::array{A, B, C}) \ + { \ run_permutation_assign(dst, src, scale); \ - return true; \ + return true; \ } XTENSOR_RUN_PERMUTATION(0, 1, 2) XTENSOR_RUN_PERMUTATION(0, 2, 1) @@ -334,10 +338,7 @@ namespace xt const auto& view = std::get<0>(cast_expression.arguments()); using view_type = std::decay_t; using cast_result_type = typename detail::cast_expression::result_type; - if constexpr ( - detail::permutation_view::value - && std::is_same::value - && std::is_floating_point::value) + if constexpr (detail::permutation_view::value && std::is_same::value && std::is_floating_point::value) { return detail::try_3d_permutation_assign( dst, diff --git a/test/test_strided_assign.cpp b/test/test_strided_assign.cpp index 67a29a5ba..99452e9e0 100644 --- a/test/test_strided_assign.cpp +++ b/test/test_strided_assign.cpp @@ -201,8 +201,11 @@ namespace xt TEST(xassign_strided, transpose_cast_permutation_plan) { - xarray input = xarray::from_shape({4, 16, 3}); - xarray result = xarray::from_shape({3, 4, 16}); + xarray + input = xarray::from_shape({4, 16, 3}); + xarray result = xarray::from_shape( + {3, 4, 16} + ); for (std::size_t i = 0; i < input.size(); ++i) { input.storage()[i] = static_cast(i % 251); @@ -227,8 +230,11 @@ namespace xt } } - xarray tiled_input = xarray::from_shape({4, 5, 2}); - xarray tiled_result = xarray::from_shape({5, 4, 2}); + xarray + tiled_input = xarray::from_shape({4, 5, 2}); + xarray tiled_result = xarray::from_shape( + {5, 4, 2} + ); for (std::size_t i = 0; i < tiled_input.size(); ++i) { tiled_input.storage()[i] = static_cast(i); @@ -254,13 +260,23 @@ namespace xt xarray row = xarray::from_shape({8}); xarray column = xarray::from_shape({4, 1, 1}); xarray result = xarray::from_shape(x.shape()); - for (std::size_t i = 0; i < x.size(); ++i) x.storage()[i] = static_cast(i); - for (std::size_t i = 0; i < row.size(); ++i) row.storage()[i] = static_cast(i * 2); - for (std::size_t i = 0; i < column.size(); ++i) column.storage()[i] = static_cast(i * 3); + for (std::size_t i = 0; i < x.size(); ++i) + { + x.storage()[i] = static_cast(i); + } + for (std::size_t i = 0; i < row.size(); ++i) + { + row.storage()[i] = static_cast(i * 2); + } + for (std::size_t i = 0; i < column.size(); ++i) + { + column.storage()[i] = static_cast(i * 3); + } noalias(result) = (x + row) * 1.5 - column; - EXPECT_TRUE(strided_assign_detail::get_loop_sizes(result, (x + row) * 1.5 - column).can_do_strided_assign); + EXPECT_TRUE(strided_assign_detail::get_loop_sizes(result, (x + row) * 1.5 - column).can_do_strided_assign + ); for (std::size_t i = 0; i < 4; ++i) { for (std::size_t j = 0; j < 3; ++j) @@ -272,17 +288,36 @@ namespace xt } } - xarray cx = xarray::from_shape({8, 3, 4}); - xarray leading = xarray::from_shape({8, 1, 1}); - xarray rest = xarray::from_shape({1, 3, 4}); - xarray cresult = xarray::from_shape(cx.shape()); - for (std::size_t i = 0; i < cx.size(); ++i) cx.storage()[i] = static_cast(i); - for (std::size_t i = 0; i < leading.size(); ++i) leading.storage()[i] = static_cast(i * 2); - for (std::size_t i = 0; i < rest.size(); ++i) rest.storage()[i] = static_cast(i * 3); + xarray cx = xarray::from_shape( + {8, 3, 4} + ); + xarray leading = xarray::from_shape( + {8, 1, 1} + ); + xarray rest = xarray::from_shape( + {1, 3, 4} + ); + xarray cresult = xarray::from_shape( + cx.shape() + ); + for (std::size_t i = 0; i < cx.size(); ++i) + { + cx.storage()[i] = static_cast(i); + } + for (std::size_t i = 0; i < leading.size(); ++i) + { + leading.storage()[i] = static_cast(i * 2); + } + for (std::size_t i = 0; i < rest.size(); ++i) + { + rest.storage()[i] = static_cast(i * 3); + } noalias(cresult) = (cx + leading) * 1.5 - rest; - EXPECT_TRUE(strided_assign_detail::get_loop_sizes(cresult, (cx + leading) * 1.5 - rest).can_do_strided_assign); + EXPECT_TRUE( + strided_assign_detail::get_loop_sizes(cresult, (cx + leading) * 1.5 - rest).can_do_strided_assign + ); for (std::size_t i = 0; i < 8; ++i) { for (std::size_t j = 0; j < 3; ++j) @@ -298,12 +333,22 @@ namespace xt TEST(xassign_strided, runtime_broadcast_plan) { - xarray x = xarray::from_shape({3, 4, 5}); - xarray row = {1.0, 2.0, 3.0, 4.0, 5.0}; - xarray column = xarray::from_shape({3, 1, 1}); - xarray result = xarray::from_shape(x.shape()); - for (std::size_t i = 0; i < x.size(); ++i) x.storage()[i] = static_cast(i); - for (std::size_t i = 0; i < column.size(); ++i) column.storage()[i] = static_cast(i * 3); + xarray x = xarray::from_shape({3, 4, 5}); + xarray row = {1.0, 2.0, 3.0, 4.0, 5.0}; + xarray column = xarray::from_shape( + {3, 1, 1} + ); + xarray result = xarray::from_shape( + x.shape() + ); + for (std::size_t i = 0; i < x.size(); ++i) + { + x.storage()[i] = static_cast(i); + } + for (std::size_t i = 0; i < column.size(); ++i) + { + column.storage()[i] = static_cast(i * 3); + } auto expression = (x + row) * 1.5 - column; EXPECT_TRUE(detail::try_runtime_plan(result, expression, false)); From 005bb0e0774ea2d9e7e0779199351963c2aba2d4 Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Mon, 5 Oct 2026 12:45:14 +0200 Subject: [PATCH 6/8] Restrict cast SIMD to floating-point batches --- include/xtensor/core/xoperation.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/xtensor/core/xoperation.hpp b/include/xtensor/core/xoperation.hpp index 74e77bbc4..e63403121 100644 --- a/include/xtensor/core/xoperation.hpp +++ b/include/xtensor/core/xoperation.hpp @@ -159,7 +159,10 @@ namespace xt template < class B, - std::enable_if_t::type>::value, int> = 0> + std::enable_if_t< + std::is_floating_point::value + && std::is_same::type>::value, + int> = 0> constexpr B simd_apply(const B& arg) const { return arg; From feb92e87bb294b0f8799b82f570c7130c956162b Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Mon, 5 Oct 2026 13:23:43 +0200 Subject: [PATCH 7/8] Advance scalar tails by logical stride --- include/xtensor/core/xassign.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/include/xtensor/core/xassign.hpp b/include/xtensor/core/xassign.hpp index 9c3ebde0f..e06a9fc5d 100644 --- a/include/xtensor/core/xassign.hpp +++ b/include/xtensor/core/xassign.hpp @@ -1377,6 +1377,7 @@ namespace xt { step_dim = cut; } + const std::size_t leading_dim = is_row_major ? loop_sizes.dimension - 1 : 0; #if defined(XTENSOR_USE_OPENMP) && defined(strided_parallel_assign) if (outer_loop_size >= XTENSOR_OPENMP_TRESHOLD / inner_loop_size) { @@ -1405,8 +1406,8 @@ namespace xt for (std::size_t i = 0; i < simd_rest; ++i) { *(res_stepper) = conditional_cast(*(fct_stepper)); - res_stepper.step_leading(); - fct_stepper.step_leading(); + res_stepper.step(leading_dim); + fct_stepper.step(leading_dim); } // next unaligned index @@ -1480,8 +1481,8 @@ namespace xt for (std::size_t i = 0; i < simd_rest; ++i) { *(res_stepper) = conditional_cast(*(fct_stepper)); - res_stepper.step_leading(); - fct_stepper.step_leading(); + res_stepper.step(leading_dim); + fct_stepper.step(leading_dim); } // next unaligned index @@ -1526,8 +1527,8 @@ namespace xt for (std::size_t i = 0; i < simd_rest; ++i) { *(res_stepper) = conditional_cast(*(fct_stepper)); - res_stepper.step_leading(); - fct_stepper.step_leading(); + res_stepper.step(leading_dim); + fct_stepper.step(leading_dim); } is_row_major From 68f55416adb889291ae9ae50358eff2ee707f393 Mon Sep 17 00:00:00 2001 From: Wolf Vollprecht Date: Mon, 5 Oct 2026 13:53:01 +0200 Subject: [PATCH 8/8] Limit optimized assignment to planned expressions --- include/xtensor/core/xassign.hpp | 56 +++-------------- include/xtensor/core/xiterator.hpp | 25 -------- include/xtensor/core/xoperation.hpp | 11 ---- include/xtensor/core/xstrides.hpp | 6 +- include/xtensor/views/xstrided_view.hpp | 13 ---- test/test_strided_assign.cpp | 82 ------------------------- 6 files changed, 11 insertions(+), 182 deletions(-) diff --git a/include/xtensor/core/xassign.hpp b/include/xtensor/core/xassign.hpp index e06a9fc5d..3ddd80ba3 100644 --- a/include/xtensor/core/xassign.hpp +++ b/include/xtensor/core/xassign.hpp @@ -328,20 +328,6 @@ namespace xt static constexpr bool value = std::conjunction>...>::value; }; - template - struct has_strided_simd_interface : has_simd_interface - { - }; - - template - struct has_strided_simd_interface, T> - : std::conjunction< - has_simd_type, - has_simd_apply>, - has_strided_simd_interface, T>...> - { - }; - template ::value, class = void> struct runtime_plan_expression : std::false_type { @@ -523,17 +509,6 @@ namespace xt return false; } - template - struct supports_row_major_gather : std::bool_constant - { - }; - - template - struct supports_row_major_gather> - : std::conjunction>...> - { - }; - /** * Considering the assignment LHS = RHS, if the requested value type used for * loading simd from RHS is not complex while LHS value_type is complex, @@ -637,9 +612,7 @@ namespace xt static constexpr bool simd_strided_assign() { - return strided_assign() && convertible_types() && simd_size() - && detail::has_strided_simd_interface::value - && detail::has_strided_simd_interface::value; + return strided_assign() && simd_assign(); } static constexpr bool simd_linear_assign(const E1& e1, const E2& e2) @@ -1266,18 +1239,8 @@ namespace xt if (is_row_major) { - if constexpr (possible && E1::contiguous_layout && E1::static_layout == layout_type::row_major && xt::detail::supports_row_major_gather::value) - { - // SIMD steppers can broadcast and gather arbitrary row-major input strides. - cut = e1.strides().size() - 1; - } - else - { - auto csf = check_strides_functor( - e1.strides() - ); - cut = csf(e2); - } + auto csf = check_strides_functor(e1.strides()); + cut = csf(e2); // This makes that only one dimension will be treated in the inner loop. if (cut < e1.strides().size() - 1) @@ -1377,7 +1340,6 @@ namespace xt { step_dim = cut; } - const std::size_t leading_dim = is_row_major ? loop_sizes.dimension - 1 : 0; #if defined(XTENSOR_USE_OPENMP) && defined(strided_parallel_assign) if (outer_loop_size >= XTENSOR_OPENMP_TRESHOLD / inner_loop_size) { @@ -1406,8 +1368,8 @@ namespace xt for (std::size_t i = 0; i < simd_rest; ++i) { *(res_stepper) = conditional_cast(*(fct_stepper)); - res_stepper.step(leading_dim); - fct_stepper.step(leading_dim); + res_stepper.step_leading(); + fct_stepper.step_leading(); } // next unaligned index @@ -1481,8 +1443,8 @@ namespace xt for (std::size_t i = 0; i < simd_rest; ++i) { *(res_stepper) = conditional_cast(*(fct_stepper)); - res_stepper.step(leading_dim); - fct_stepper.step(leading_dim); + res_stepper.step_leading(); + fct_stepper.step_leading(); } // next unaligned index @@ -1527,8 +1489,8 @@ namespace xt for (std::size_t i = 0; i < simd_rest; ++i) { *(res_stepper) = conditional_cast(*(fct_stepper)); - res_stepper.step(leading_dim); - fct_stepper.step(leading_dim); + res_stepper.step_leading(); + fct_stepper.step_leading(); } is_row_major diff --git a/include/xtensor/core/xiterator.hpp b/include/xtensor/core/xiterator.hpp index 2805afb14..448f6093d 100644 --- a/include/xtensor/core/xiterator.hpp +++ b/include/xtensor/core/xiterator.hpp @@ -566,31 +566,6 @@ namespace xt inline auto xstepper::step_simd() -> simd_return_type { using simd_type = simd_return_type; - const auto& strides = p_c->strides(); - const bool broadcast = !strides.empty() - && (p_c->layout() == layout_type::column_major ? strides.front() : strides.back() - ) == 0; - if (broadcast) - { - return simd_type(*m_it); - } -#ifdef XTENSOR_USE_XSIMD - if constexpr (!xt_simd::is_batch_bool::value && !xt_simd::is_batch_complex::value) - { - const auto stride = p_c->layout() == layout_type::column_major ? strides.front() : strides.back(); - if (stride != 1) - { - using index_value_type = xsimd::as_integer_t; - using arch_type = typename simd_type::arch_type; - xsimd::batch - index = xsimd::make_iota_batch_constant(); - index *= static_cast(stride); - simd_type reg = simd_type::gather(&(*m_it), index); - m_it += static_cast(xt_simd::revert_simd_traits::size * stride); - return reg; - } - } -#endif simd_type reg = detail::step_simd_invoker::template apply(m_it); m_it += xt_simd::revert_simd_traits::size; return reg; diff --git a/include/xtensor/core/xoperation.hpp b/include/xtensor/core/xoperation.hpp index e63403121..35d05da75 100644 --- a/include/xtensor/core/xoperation.hpp +++ b/include/xtensor/core/xoperation.hpp @@ -156,17 +156,6 @@ namespace xt { return static_cast(arg); } - - template < - class B, - std::enable_if_t< - std::is_floating_point::value - && std::is_same::type>::value, - int> = 0> - constexpr B simd_apply(const B& arg) const - { - return arg; - } }; }; diff --git a/include/xtensor/core/xstrides.hpp b/include/xtensor/core/xstrides.hpp index 278a93649..998ffcd6b 100644 --- a/include/xtensor/core/xstrides.hpp +++ b/include/xtensor/core/xstrides.hpp @@ -810,8 +810,7 @@ namespace xt for (; s2_index != 0; --s1_index, --s2_index) { - if (static_cast(s2[s2_index - 1]) != value_type(0) - && static_cast(s1[s1_index - 1]) != static_cast(s2[s2_index - 1])) + if (static_cast(s1[s1_index - 1]) != static_cast(s2[s2_index - 1])) { break; } @@ -842,8 +841,7 @@ namespace xt for (; index < size; ++index) { - if (static_cast(s2[index]) != value_type(0) - && static_cast(s1[index]) != static_cast(s2[index])) + if (static_cast(s1[index]) != static_cast(s2[index])) { break; } diff --git a/include/xtensor/views/xstrided_view.hpp b/include/xtensor/views/xstrided_view.hpp index 586a67133..b0fe2c6e0 100644 --- a/include/xtensor/views/xstrided_view.hpp +++ b/include/xtensor/views/xstrided_view.hpp @@ -313,19 +313,6 @@ namespace xt friend class xaxis_slice_iterator; }; - namespace detail - { - template - struct has_strided_simd_interface, T> - : std::conjunction< - has_simd_type, - std::bool_constant::value_type, T>::value>, - std::is_lvalue_reference::const_reference>, - std::bool_constant::value>> - { - }; - } - /************************** * xstrided_view builders * **************************/ diff --git a/test/test_strided_assign.cpp b/test/test_strided_assign.cpp index 99452e9e0..b4e9df304 100644 --- a/test/test_strided_assign.cpp +++ b/test/test_strided_assign.cpp @@ -215,10 +215,6 @@ namespace xt noalias(result) = expression; EXPECT_TRUE(try_permutation_assign(result, expression)); -#ifdef XTENSOR_USE_XSIMD - EXPECT_TRUE((xassign_traits::simd_strided_assign())); - EXPECT_TRUE(strided_assign_detail::get_loop_sizes(result, expression).can_do_strided_assign); -#endif for (std::size_t c = 0; c < 3; ++c) { for (std::size_t i = 0; i < 4; ++i) @@ -253,84 +249,6 @@ namespace xt } } -#ifdef XTENSOR_USE_XSIMD - TEST(xassign_strided, simd_broadcast) - { - xarray x = xarray::from_shape({4, 3, 8}); - xarray row = xarray::from_shape({8}); - xarray column = xarray::from_shape({4, 1, 1}); - xarray result = xarray::from_shape(x.shape()); - for (std::size_t i = 0; i < x.size(); ++i) - { - x.storage()[i] = static_cast(i); - } - for (std::size_t i = 0; i < row.size(); ++i) - { - row.storage()[i] = static_cast(i * 2); - } - for (std::size_t i = 0; i < column.size(); ++i) - { - column.storage()[i] = static_cast(i * 3); - } - - noalias(result) = (x + row) * 1.5 - column; - - EXPECT_TRUE(strided_assign_detail::get_loop_sizes(result, (x + row) * 1.5 - column).can_do_strided_assign - ); - for (std::size_t i = 0; i < 4; ++i) - { - for (std::size_t j = 0; j < 3; ++j) - { - for (std::size_t k = 0; k < 8; ++k) - { - EXPECT_EQ(result(i, j, k), (x(i, j, k) + row(k)) * 1.5 - column(i, 0, 0)); - } - } - } - - xarray cx = xarray::from_shape( - {8, 3, 4} - ); - xarray leading = xarray::from_shape( - {8, 1, 1} - ); - xarray rest = xarray::from_shape( - {1, 3, 4} - ); - xarray cresult = xarray::from_shape( - cx.shape() - ); - for (std::size_t i = 0; i < cx.size(); ++i) - { - cx.storage()[i] = static_cast(i); - } - for (std::size_t i = 0; i < leading.size(); ++i) - { - leading.storage()[i] = static_cast(i * 2); - } - for (std::size_t i = 0; i < rest.size(); ++i) - { - rest.storage()[i] = static_cast(i * 3); - } - - noalias(cresult) = (cx + leading) * 1.5 - rest; - - EXPECT_TRUE( - strided_assign_detail::get_loop_sizes(cresult, (cx + leading) * 1.5 - rest).can_do_strided_assign - ); - for (std::size_t i = 0; i < 8; ++i) - { - for (std::size_t j = 0; j < 3; ++j) - { - for (std::size_t k = 0; k < 4; ++k) - { - EXPECT_EQ(cresult(i, j, k), (cx(i, j, k) + leading(i, 0, 0)) * 1.5 - rest(0, j, k)); - } - } - } - } -#endif - TEST(xassign_strided, runtime_broadcast_plan) { xarray x = xarray::from_shape({3, 4, 5});