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..a73902685
--- /dev/null
+++ b/benchmark/benchmark_compare.cpp
@@ -0,0 +1,309 @@
+/***************************************************************************
+ * 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 <algorithm>
+#include <cstddef>
+#include <vector>
+
+#include <benchmark/benchmark.h>
+
+#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 <Eigen/Core>
+#include <unsupported/Eigen/CXX11/Tensor>
+#endif
+
+#ifdef XTENSOR_BENCHMARK_USE_ARMADILLO
+#define ARMA_DONT_PRINT_FAST_MATH_WARNING
+#include <armadillo>
+#endif
+
+namespace xt::compare
+{
+    constexpr std::size_t dynamic_size = 64 * 64 * 16;
+    constexpr std::size_t fixed_size = 16 * 16;
+
+    template <class X, class Y, class Z>
+    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<std::size_t> shape = {64, 64, 16};
+        xarray<double> 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<double> 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<std::size_t> shape = {64, 64, 16};
+        xarray<double> x, row, column, out;
+        x.resize(shape);
+        row.resize({16});
+        column.resize({64, 1, 1});
+        out.resize(shape);
+        xarray<double> 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<double> 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<unsigned char> input(height * width * channels);
+        std::vector<float> out(input.size());
+        for (std::size_t i = 0; i < input.size(); ++i)
+        {
+            input[i] = static_cast<unsigned char>(i % 251);
+        }
+        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<unsigned char> input = xarray<unsigned char>::from_shape({128, 256, 3});
+        xarray<float> out = xarray<float>::from_shape({3, 128, 256});
+        for (std::size_t i = 0; i < input.size(); ++i)
+        {
+            input.storage()[i] = static_cast<unsigned char>(i % 251);
+        }
+        for (auto _ : state)
+        {
+            noalias(out) = cast<float>(transpose(input, {2, 0, 1})) / 255.0f;
+            benchmark::DoNotOptimize(out.data());
+            benchmark::ClobberMemory();
+        }
+    }
+
+    void linear_fixed_xtensor(benchmark::State& state)
+    {
+        xtensor_fixed<double, xshape<16, 16>> 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<double, fixed_size, 1> 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();
+        }
+    }
+
+    void broadcast_dynamic_eigen(benchmark::State& state)
+    {
+        using tensor = Eigen::Tensor<double, 3, Eigen::RowMajor>;
+        tensor x(64, 64, 16), out(64, 64, 16);
+        Eigen::Tensor<double, 1, Eigen::RowMajor> 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<Eigen::Index, 3> row_shape = {1, 1, 16};
+        const Eigen::array<Eigen::Index, 3> row_broadcast = {64, 64, 1};
+        const Eigen::array<Eigen::Index, 3> column_shape = {64, 1, 1};
+        const Eigen::array<Eigen::Index, 3> 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<unsigned char, 3, Eigen::RowMajor> input(128, 256, 3);
+        Eigen::Tensor<float, 3, Eigen::RowMajor> out(3, 128, 256);
+        for (Eigen::Index i = 0; i < input.size(); ++i)
+        {
+            input.data()[i] = static_cast<unsigned char>(i % 251);
+        }
+        const Eigen::array<Eigen::Index, 3> permutation = {2, 0, 1};
+        for (auto _ : state)
+        {
+            out = input.shuffle(permutation).cast<float>() / 255.0f;
+            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<fixed_size> 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);
+    BENCHMARK(broadcast_dynamic_eigen);
+    BENCHMARK(transpose_cast_eigen);
+#endif
+#ifdef XTENSOR_BENCHMARK_USE_ARMADILLO
+    BENCHMARK(linear_dynamic_armadillo);
+    BENCHMARK(linear_fixed_armadillo);
+#endif
+}
diff --git a/include/xtensor/core/xassign.hpp b/include/xtensor/core/xassign.hpp
index 174ddcdc0..3ddd80ba3 100644
--- a/include/xtensor/core/xassign.hpp
+++ b/include/xtensor/core/xassign.hpp
@@ -11,7 +11,9 @@
 #define XTENSOR_ASSIGN_HPP
 
 #include <algorithm>
+#include <array>
 #include <functional>
+#include <tuple>
 #include <type_traits>
 #include <utility>
 
@@ -58,6 +60,12 @@ namespace xt
     template <class E1, class E2>
     void strided_assign(E1& e1, const E2& e2, std::true_type /*enable*/);
 
+    template <class E1, class E2>
+    bool try_permutation_assign(E1&, const E2&)
+    {
+        return false;
+    }
+
     /************************
      * xexpression_assigner *
      ************************/
@@ -320,6 +328,187 @@ namespace xt
             static constexpr bool value = std::conjunction<use_strided_loop<std::decay_t<CT>>...>::value;
         };
 
+        template <class E, bool = !std::is_same<typename E::value_type, bool>::value, class = void>
+        struct runtime_plan_expression : std::false_type
+        {
+        };
+
+        template <class E>
+        struct runtime_plan_expression<E, true, std::void_t<decltype(std::declval<const E&>().storage().data())>>
+            : std::bool_constant<
+                  data_interface_expression<E> && std::is_arithmetic<typename E::value_type>::value
+                  && std::is_pointer<decltype(std::declval<const E&>().storage().data())>::value
+                  && std::is_same<
+                      std::remove_cv_t<std::remove_pointer_t<decltype(std::declval<const E&>().storage().data())>>,
+                      typename E::value_type>::value>
+        {
+        };
+
+        template <class T>
+        struct runtime_plan_expression<xscalar<T>, true, void>
+            : std::bool_constant<std::is_arithmetic<T>::value>
+        {
+        };
+
+        template <class F, class... CT>
+        struct runtime_plan_expression<xfunction<F, CT...>, true, void>
+            : std::conjunction<runtime_plan_expression<std::decay_t<CT>>...>
+        {
+        };
+
+        template <class E, bool = !std::is_same<typename E::value_type, bool>::value, class = void>
+        struct runtime_plan_output : std::false_type
+        {
+        };
+
+        template <class E>
+        struct runtime_plan_output<
+            E,
+            true,
+            std::void_t<decltype(std::declval<E&>().storage().data()), decltype(std::declval<E&>().data_offset())>>
+            : std::bool_constant<
+                  E::contiguous_layout && E::static_layout == layout_type::row_major
+                  && std::is_pointer<decltype(std::declval<E&>().storage().data())>::value
+                  && std::is_same<
+                      std::remove_cv_t<std::remove_pointer_t<decltype(std::declval<E&>().storage().data())>>,
+                      typename E::value_type>::value>
+        {
+        };
+
+        template <class E>
+        struct runtime_leaf_plan
+        {
+            using value_type = typename E::value_type;
+            const value_type* data;
+            std::array<std::ptrdiff_t, 3> strides;
+
+            value_type operator()(std::size_t i, std::size_t j, std::size_t k) const
+            {
+                const auto offset = static_cast<std::ptrdiff_t>(i) * strides[0]
+                                    + static_cast<std::ptrdiff_t>(j) * strides[1];
+                if (strides[2] == 0)
+                {
+                    return data[offset];
+                }
+                if (strides[2] == 1)
+                {
+                    return data[offset + static_cast<std::ptrdiff_t>(k)];
+                }
+                return data[offset + static_cast<std::ptrdiff_t>(k) * strides[2]];
+            }
+        };
+
+        template <class T>
+        struct runtime_scalar_plan
+        {
+            T value;
+
+            T operator()(std::size_t, std::size_t, std::size_t) const
+            {
+                return value;
+            }
+        };
+
+        template <class F, class P>
+        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 <class E>
+        struct runtime_plan_builder
+        {
+            static auto make(const E& expression)
+            {
+                std::array<std::ptrdiff_t, 3> 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<E>{expression.storage().data() + expression.data_offset(), strides};
+            }
+        };
+
+        template <class T>
+        struct runtime_plan_builder<xscalar<T>>
+        {
+            static auto make(const xscalar<T>& expression)
+            {
+                return runtime_scalar_plan<T>{expression()};
+            }
+        };
+
+        template <class F, class... CT>
+        struct runtime_plan_builder<xfunction<F, CT...>>
+        {
+            using expression_type = xfunction<F, CT...>;
+
+            template <std::size_t... I>
+            static auto make(const expression_type& expression, std::index_sequence<I...>)
+            {
+                auto arguments = std::make_tuple(
+                    runtime_plan_builder<std::decay_t<decltype(std::get<I>(expression.arguments()))>>::make(
+                        std::get<I>(expression.arguments())
+                    )...
+                );
+                return runtime_function_plan<std::decay_t<F>, decltype(arguments)>{
+                    expression.functor(),
+                    std::move(arguments)
+                };
+            }
+
+            static auto make(const expression_type& expression)
+            {
+                return make(expression, std::index_sequence_for<CT...>{});
+            }
+        };
+
+        template <class E1, class E2>
+        bool try_runtime_plan(E1& dst, const E2& expression, bool trivial)
+        {
+            if (trivial)
+            {
+                return false;
+            }
+            if constexpr (runtime_plan_output<E1>::value && runtime_plan_expression<E2>::value)
+            {
+                if (dst.dimension() == 3 && expression.dimension() <= 3)
+                {
+                    auto plan = runtime_plan_builder<E2>::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;
+        }
+
         /**
          * 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,
@@ -447,6 +636,11 @@ namespace xt
         const E2& de2 = e2.derived_cast();
         using traits = xassign_traits<E1, E2>;
 
+        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();
@@ -1047,6 +1241,7 @@ namespace xt
             {
                 auto csf = check_strides_functor<layout_type::row_major, decltype(e1.strides())>(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 +1249,7 @@ namespace xt
                     cut = e1.strides().size() - 1;
                 }
             }
-            else if (!is_row_major)
+            else
             {
                 auto csf = check_strides_functor<layout_type::column_major, decltype(e1.strides())>(e1.strides()
                 );
@@ -1064,7 +1259,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::size_t>(std::accumulate(
                 e1.shape().begin(),
diff --git a/include/xtensor/core/xoperation.hpp b/include/xtensor/core/xoperation.hpp
index 76bfc7c50..35d05da75 100644
--- a/include/xtensor/core/xoperation.hpp
+++ b/include/xtensor/core/xoperation.hpp
@@ -11,6 +11,7 @@
 #define XTENSOR_OPERATION_HPP
 
 #include <algorithm>
+#include <array>
 #include <type_traits>
 
 #include <xtl/xsequence.hpp>
@@ -148,22 +149,133 @@ namespace xt
             struct functor
             {
                 using result_type = R;
+                using cast_result_type = R;
 
                 template <class A1>
                 constexpr result_type operator()(const A1& arg) const
                 {
                     return static_cast<R>(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
-                {
-                    return static_cast<R>(arg);
-                }*/
             };
         };
 
+        template <class E, class = void>
+        struct cast_expression : std::false_type
+        {
+        };
+
+        template <class F, class CT>
+        struct cast_expression<xfunction<F, CT>, std::void_t<typename F::cast_result_type>> : std::true_type
+        {
+            using result_type = typename F::cast_result_type;
+        };
+
+        template <class E>
+        struct permutation_view : std::false_type
+        {
+        };
+
+        template <class CT, class S, layout_type L, class FST>
+        struct permutation_view<xstrided_view<CT, S, L, FST>>
+            : std::bool_constant<L == layout_type::dynamic && is_specialization_of<inner_storage_getter, FST>::value>
+        {
+        };
+
+        template <std::size_t A, std::size_t B, std::size_t C, class E1, class V, class R, class S>
+        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<std::size_t>(src.strides()[A]);
+            const std::size_t src_stride_b = static_cast<std::size_t>(src.strides()[B]);
+            const std::size_t src_stride_c = static_cast<std::size_t>(src.strides()[C]);
+            const std::size_t dst_stride_a = static_cast<std::size_t>(dst.strides()[A]);
+            const std::size_t dst_stride_b = static_cast<std::size_t>(dst.strides()[B]);
+            const std::size_t dst_stride_c = static_cast<std::size_t>(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<R>(src_data[3 * i]) / static_cast<R>(scale);
+                    dst_1[i] = static_cast<R>(src_data[3 * i + 1]) / static_cast<R>(scale);
+                    dst_2[i] = static_cast<R>(src_data[3 * i + 2]) / static_cast<R>(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<R>(src_data[src_offset])
+                                                       / static_cast<R>(scale);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        template <class E1, class V, class R, class S>
+        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<std::size_t, 3> 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<std::size_t>(src.strides()[axes[1]]) != src.shape()[axes[2]]
+                || static_cast<std::size_t>(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<std::size_t, 3>{A, B, C})                \
+    {                                                               \
+        run_permutation_assign<A, B, C, E1, V, R>(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 <class Tag, class F, class... E>
         struct select_xfunction_expression;
 
@@ -208,6 +320,28 @@ namespace xt
             enable_if_t<has_xexpression<std::decay_t<E>...>::value, xfunction_type<F, E...>>::type;
     }
 
+    template <class E1, class CT1, class CT2>
+    inline bool try_permutation_assign(E1& dst, const xfunction<detail::divides, CT1, CT2>& expression)
+    {
+        using cast_expression_type = std::decay_t<decltype(std::get<0>(expression.arguments()))>;
+        if constexpr (detail::cast_expression<cast_expression_type>::value && is_xscalar<std::decay_t<CT2>>::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<decltype(view)>;
+            using cast_result_type = typename detail::cast_expression<cast_expression_type>::result_type;
+            if constexpr (detail::permutation_view<view_type>::value && std::is_same<cast_result_type, typename E1::value_type>::value && std::is_floating_point<cast_result_type>::value)
+            {
+                return detail::try_3d_permutation_assign<E1, view_type, cast_result_type>(
+                    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 a0131b27d..b4e9df304 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"
@@ -196,4 +198,87 @@ namespace xt
             }
         }
     }
+
+    TEST(xassign_strided, transpose_cast_permutation_plan)
+    {
+        xarray<unsigned char, layout_type::row_major>
+            input = xarray<unsigned char, layout_type::row_major>::from_shape({4, 16, 3});
+        xarray<float, layout_type::row_major> result = xarray<float, layout_type::row_major>::from_shape(
+            {3, 4, 16}
+        );
+        for (std::size_t i = 0; i < input.size(); ++i)
+        {
+            input.storage()[i] = static_cast<unsigned char>(i % 251);
+        }
+
+        auto expression = cast<float>(transpose(input, {2, 0, 1})) / 255.0f;
+        noalias(result) = expression;
+
+        EXPECT_TRUE(try_permutation_assign(result, expression));
+        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<float>(input(i, j, c)) / 255.0f);
+                }
+            }
+        }
+
+        xarray<unsigned char, layout_type::row_major>
+            tiled_input = xarray<unsigned char, layout_type::row_major>::from_shape({4, 5, 2});
+        xarray<float, layout_type::row_major> tiled_result = xarray<float, layout_type::row_major>::from_shape(
+            {5, 4, 2}
+        );
+        for (std::size_t i = 0; i < tiled_input.size(); ++i)
+        {
+            tiled_input.storage()[i] = static_cast<unsigned char>(i);
+        }
+        auto tiled_expression = cast<float>(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<float>(tiled_input(j, i, c)) / 7.0f);
+                }
+            }
+        }
+    }
+
+    TEST(xassign_strided, runtime_broadcast_plan)
+    {
+        xarray<double, layout_type::row_major> x = xarray<double, layout_type::row_major>::from_shape({3, 4, 5});
+        xarray<double, layout_type::row_major> row = {1.0, 2.0, 3.0, 4.0, 5.0};
+        xarray<double, layout_type::row_major> column = xarray<double, layout_type::row_major>::from_shape(
+            {3, 1, 1}
+        );
+        xarray<double, layout_type::row_major> result = xarray<double, layout_type::row_major>::from_shape(
+            x.shape()
+        );
+        for (std::size_t i = 0; i < x.size(); ++i)
+        {
+            x.storage()[i] = static_cast<double>(i);
+        }
+        for (std::size_t i = 0; i < column.size(); ++i)
+        {
+            column.storage()[i] = static_cast<double>(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));
+                }
+            }
+        }
+    }
 }