diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index 24c391fb36..b6d83609a3 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -108,6 +108,9 @@
 
 ### Bug fixes
 
+* Fix wire order permutations when using `qml.probs` with out-of-order wires.
+  [(#707)](https://github.com/PennyLaneAI/pennylane-lightning/pull/707)
+
 * Lightning Qubit once again respects the wire order specified on device instantiation.
   [(#705)](https://github.com/PennyLaneAI/pennylane-lightning/pull/705)
 
diff --git a/.github/workflows/tests_linux_cpp.yml b/.github/workflows/tests_linux_cpp.yml
index 2a4d2c7058..521c9ab0e7 100644
--- a/.github/workflows/tests_linux_cpp.yml
+++ b/.github/workflows/tests_linux_cpp.yml
@@ -142,7 +142,7 @@ jobs:
 
       - name: Install dependencies
         run: |
-          sudo apt-get update && sudo apt-get -y -q install cmake gcc-$GCC_VERSION  g++-$GCC_VERSION libopenblas-dev ninja-build gcovr lcov
+          sudo apt-get update && sudo apt-get -y -q install cmake gcc-$GCC_VERSION  g++-$GCC_VERSION libopenblas-base libopenblas-dev ninja-build gcovr lcov
           python -m pip install scipy
 
       - name: Build and run unit tests
diff --git a/.github/workflows/tests_lqcpu_python.yml b/.github/workflows/tests_lqcpu_python.yml
index a9d9cdb125..5c48006936 100644
--- a/.github/workflows/tests_lqcpu_python.yml
+++ b/.github/workflows/tests_lqcpu_python.yml
@@ -73,7 +73,7 @@ jobs:
 
       - name: Install dependencies
         run: |
-          sudo apt-get update && sudo apt-get -y -q install cmake gcc-$GCC_VERSION g++-$GCC_VERSION libopenblas-dev
+          sudo apt-get update && sudo apt-get -y -q install cmake gcc-$GCC_VERSION g++-$GCC_VERSION libopenblas-base libopenblas-dev
           python -m pip install scipy wheel
 
       - name: Get required Python packages
diff --git a/README.rst b/README.rst
index 63980f2b8f..d2df483c03 100644
--- a/README.rst
+++ b/README.rst
@@ -455,5 +455,6 @@ PennyLane Lightning makes use of the following libraries and tools, which are un
 - **pybind11:** https://github.com/pybind/pybind11
 - **Kokkos Core:** https://github.com/kokkos/kokkos
 - **NVIDIA cuQuantum:** https://developer.nvidia.com/cuquantum-sdk
+- **Xanadu JET:** https://github.com/XanaduAI/jet
 
 .. acknowledgements-end-inclusion-marker-do-not-remove
diff --git a/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp b/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp
index ff28ba9c62..a4ca85d556 100644
--- a/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp
+++ b/pennylane_lightning/core/src/measurements/tests/Test_MeasurementsBase.cpp
@@ -84,7 +84,7 @@ template <typename TypeList> void testProbabilities() {
         // Expected results calculated with Pennylane default.qubit:
         std::vector<std::pair<std::vector<size_t>, std::vector<PrecisionT>>>
             input = {
-#ifdef _ENABLE_PLGPU
+#if defined(_ENABLE_PLGPU)
                 // Bit index reodering conducted in the python layer
                 // for L-GPU. Also L-GPU backend doesn't support
                 // out of order wires for probability calculation
@@ -92,9 +92,9 @@ template <typename TypeList> void testProbabilities() {
                  {0.67078706, 0.03062806, 0.0870997, 0.00397696, 0.17564072,
                   0.00801973, 0.02280642, 0.00104134}}
 #else
-                {{0, 1, 2},
-                 {0.67078706, 0.03062806, 0.0870997, 0.00397696, 0.17564072,
-                  0.00801973, 0.02280642, 0.00104134}},
+#if defined(_ENABLE_PLQUBIT)
+                // LightningQubit currently supports arbitrary wire index
+                // ordering.
                 {{0, 2, 1},
                  {0.67078706, 0.0870997, 0.03062806, 0.00397696, 0.17564072,
                   0.02280642, 0.00801973, 0.00104134}},
@@ -102,18 +102,23 @@ template <typename TypeList> void testProbabilities() {
                  {0.67078706, 0.03062806, 0.17564072, 0.00801973, 0.0870997,
                   0.00397696, 0.02280642, 0.00104134}},
                 {{1, 2, 0},
-                 {0.67078706, 0.0870997, 0.17564072, 0.02280642, 0.03062806,
-                  0.00397696, 0.00801973, 0.00104134}},
-                {{2, 0, 1},
                  {0.67078706, 0.17564072, 0.03062806, 0.00801973, 0.0870997,
                   0.02280642, 0.00397696, 0.00104134}},
+                {{2, 0, 1},
+                 {0.67078706, 0.0870997, 0.17564072, 0.02280642, 0.03062806,
+                  0.00397696, 0.00801973, 0.00104134}},
                 {{2, 1, 0},
                  {0.67078706, 0.17564072, 0.0870997, 0.02280642, 0.03062806,
                   0.00801973, 0.00397696, 0.00104134}},
+                {{2, 1}, {0.84642778, 0.10990612, 0.0386478, 0.0050183}},
+
+#endif
+                {{0, 1, 2},
+                 {0.67078706, 0.03062806, 0.0870997, 0.00397696, 0.17564072,
+                  0.00801973, 0.02280642, 0.00104134}},
                 {{0, 1}, {0.70141512, 0.09107666, 0.18366045, 0.02384776}},
                 {{0, 2}, {0.75788676, 0.03460502, 0.19844714, 0.00906107}},
                 {{1, 2}, {0.84642778, 0.0386478, 0.10990612, 0.0050183}},
-                {{2, 1}, {0.84642778, 0.10990612, 0.0386478, 0.0050183}},
                 {{0}, {0.79249179, 0.20750821}},
                 {{1}, {0.88507558, 0.11492442}},
                 {{2}, {0.9563339, 0.0436661}}
diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp
index dc213c126f..3a80ca8990 100644
--- a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp
+++ b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp
@@ -92,6 +92,11 @@ class Measurements final
      * @return std::vector<PrecisionT>
      */
     auto probs(const std::vector<size_t> &wires) -> std::vector<PrecisionT> {
+        PL_ABORT_IF_NOT(std::is_sorted(wires.cbegin(), wires.cend()) ||
+                            std::is_sorted(wires.rbegin(), wires.rend()),
+                        "LightningGPU does not currently support out-of-order "
+                        "wire indices with probability calculations");
+
         // Data return type fixed as double in custatevec function call
         std::vector<double> probabilities(Pennylane::Util::exp2(wires.size()));
         // this should be built upon by the wires not participating
@@ -193,6 +198,10 @@ class Measurements final
 
     std::vector<PrecisionT> probs(const std::vector<size_t> &wires,
                                   size_t num_shots) {
+        PL_ABORT_IF_NOT(std::is_sorted(wires.cbegin(), wires.cend()),
+                        "LightningGPU does not currently support out-of-order "
+                        "wire indices with probability calculations");
+
         return BaseType::probs(wires, num_shots);
     }
 
diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp
index d4f85aeead..0b5cd01012 100644
--- a/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp
+++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/MeasurementsKokkos.hpp
@@ -506,6 +506,10 @@ class Measurements final
      * The basis columns are rearranged according to wires.
      */
     std::vector<PrecisionT> probs(const std::vector<size_t> &wires) {
+        PL_ABORT_IF_NOT(
+            std::is_sorted(wires.cbegin(), wires.cend()),
+            "LightningKokkos does not currently support out-of-order wire "
+            "indices with probability calculations");
         using MDPolicyType_2D =
             Kokkos::MDRangePolicy<Kokkos::Rank<2, Kokkos::Iterate::Left>>;
 
@@ -645,6 +649,11 @@ class Measurements final
 
     std::vector<PrecisionT> probs(const std::vector<size_t> &wires,
                                   size_t num_shots) {
+        PL_ABORT_IF_NOT(
+            std::is_sorted(wires.cbegin(), wires.cend()),
+            "LightningKokkos does not currently support out-of-order wire "
+            "indices with probability calculations");
+
         return BaseType::probs(wires, num_shots);
     }
 
diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/tests/Test_StateVectorKokkos_Measure.cpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/tests/Test_StateVectorKokkos_Measure.cpp
index 6add1f2140..662bafc7be 100644
--- a/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/tests/Test_StateVectorKokkos_Measure.cpp
+++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/measurements/tests/Test_StateVectorKokkos_Measure.cpp
@@ -226,25 +226,10 @@ TEMPLATE_TEST_CASE("Probabilities", "[Measures]", float, double) {
         {{0, 1, 2},
          {0.67078706, 0.03062806, 0.0870997, 0.00397696, 0.17564072, 0.00801973,
           0.02280642, 0.00104134}},
-        {{0, 2, 1},
-         {0.67078706, 0.0870997, 0.03062806, 0.00397696, 0.17564072, 0.02280642,
-          0.00801973, 0.00104134}},
-        {{1, 0, 2},
-         {0.67078706, 0.03062806, 0.17564072, 0.00801973, 0.0870997, 0.00397696,
-          0.02280642, 0.00104134}},
-        {{1, 2, 0},
-         {0.67078706, 0.0870997, 0.17564072, 0.02280642, 0.03062806, 0.00397696,
-          0.00801973, 0.00104134}},
-        {{2, 0, 1},
-         {0.67078706, 0.17564072, 0.03062806, 0.00801973, 0.0870997, 0.02280642,
-          0.00397696, 0.00104134}},
-        {{2, 1, 0},
-         {0.67078706, 0.17564072, 0.0870997, 0.02280642, 0.03062806, 0.00801973,
-          0.00397696, 0.00104134}},
+        // TODO: Fix LK out-of-order permutations
         {{0, 1}, {0.70141512, 0.09107666, 0.18366045, 0.02384776}},
         {{0, 2}, {0.75788676, 0.03460502, 0.19844714, 0.00906107}},
         {{1, 2}, {0.84642778, 0.0386478, 0.10990612, 0.0050183}},
-        {{2, 1}, {0.84642778, 0.10990612, 0.0386478, 0.0050183}},
         {{0}, {0.79249179, 0.20750821}},
         {{1}, {0.88507558, 0.11492442}},
         {{2}, {0.9563339, 0.0436661}}};
diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp
index 0e74849964..580f223367 100644
--- a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp
+++ b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/MeasurementsLQubit.hpp
@@ -31,6 +31,7 @@
 
 #include "LinearAlgebra.hpp"
 #include "MeasurementsBase.hpp"
+#include "NDPermuter.hpp"
 #include "Observables.hpp"
 #include "SparseLinAlg.hpp"
 #include "StateVectorLQubitManaged.hpp"
@@ -44,6 +45,7 @@ using namespace Pennylane::Measures;
 using namespace Pennylane::Observables;
 using Pennylane::LightningQubit::StateVectorLQubitManaged;
 using Pennylane::LightningQubit::Util::innerProdC;
+namespace PUtil = Pennylane::Util;
 } // namespace
 /// @endcond
 
@@ -100,16 +102,17 @@ class Measurements final
         // Determining index that would sort the vector.
         // This information is needed later.
         const auto sorted_ind_wires = Pennylane::Util::sorting_indices(wires);
+
         // Sorting wires.
         std::vector<size_t> sorted_wires(wires.size());
         for (size_t pos = 0; pos < wires.size(); pos++) {
             sorted_wires[pos] = wires[sorted_ind_wires[pos]];
         }
+
         // Determining probabilities for the sorted wires.
         const ComplexT *arr_data = this->_statevector.getData();
 
         size_t num_qubits = this->_statevector.getNumQubits();
-
         const std::vector<size_t> all_indices =
             Gates::generateBitPatterns(sorted_wires, num_qubits);
         const std::vector<size_t> all_offsets = Gates::generateBitPatterns(
@@ -125,12 +128,29 @@ class Measurements final
             }
             ind_probs++;
         }
-        // Transposing the probabilities tensor with the indices determined
-        // at the beginning.
+
+        // Permute the data according to the required wire ordering
         if (wires != sorted_wires) {
-            probabilities = Pennylane::Util::transpose_state_tensor(
-                probabilities, sorted_ind_wires);
+            static constexpr std::size_t CACHE_SIZE = 8;
+            PUtil::Permuter<PUtil::DefaultPermuter<CACHE_SIZE>> p{};
+            std::vector<std::size_t> shape(wires.size(), 2);
+            std::vector<std::string> wire_labels_old(sorted_wires.size(), "");
+            std::vector<std::string> wire_labels_new(wires.size(), "");
+
+            std::transform(sorted_wires.begin(), sorted_wires.end(),
+                           wire_labels_old.begin(), [](std::size_t index) {
+                               return std::to_string(index);
+                           });
+            std::transform(
+                wires.begin(), wires.end(), wire_labels_new.begin(),
+                [](std::size_t index) { return std::to_string(index); });
+
+            auto probs_sorted = probabilities;
+            p.Transpose(probabilities, shape, probs_sorted, wire_labels_old,
+                        wire_labels_new);
+            return probs_sorted;
         }
+
         return probabilities;
     }
 
diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/CMakeLists.txt b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/CMakeLists.txt
index 1ea6861b8a..0bfbe65f90 100644
--- a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/CMakeLists.txt
+++ b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/CMakeLists.txt
@@ -17,6 +17,7 @@ FetchAndIncludeCatch()
 add_library(lightning_qubit_measurements_tests INTERFACE)
 target_link_libraries(lightning_qubit_measurements_tests INTERFACE  Catch2::Catch2
                                                                     lightning_measurements
+                                                                    lightning_qubit_observables
                                                                     lightning_qubit_measurements
                                                                     )
 
diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/Test_MeasurementsLQubit.cpp b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/Test_MeasurementsLQubit.cpp
index b4be45e99a..63de0d3a95 100644
--- a/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/Test_MeasurementsLQubit.cpp
+++ b/pennylane_lightning/core/src/simulators/lightning_qubit/measurements/tests/Test_MeasurementsLQubit.cpp
@@ -20,6 +20,7 @@
 #include <catch2/catch.hpp>
 
 #include "MeasurementsLQubit.hpp"
+#include "ObservablesLQubit.hpp"
 #include "StateVectorLQubitManaged.hpp"
 #include "StateVectorLQubitRaw.hpp"
 #include "Util.hpp"
@@ -33,6 +34,7 @@ namespace {
 using namespace Pennylane::Util;
 
 using namespace Pennylane::LightningQubit;
+using namespace Pennylane::LightningQubit::Observables;
 using namespace Pennylane::LightningQubit::Measures;
 
 }; // namespace
@@ -81,17 +83,17 @@ TEMPLATE_PRODUCT_TEST_CASE("Expected Values", "[Measurements]",
         operations_list = {PauliX, PauliX, PauliX};
         exp_values = Measurer.expval(operations_list, wires_list);
         exp_values_ref = {0.49272486, 0.42073549, 0.28232124};
-        REQUIRE_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
+        CHECK_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
 
         operations_list = {PauliY, PauliY, PauliY};
         exp_values = Measurer.expval(operations_list, wires_list);
         exp_values_ref = {-0.64421768, -0.47942553, -0.29552020};
-        REQUIRE_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
+        CHECK_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
 
         operations_list = {PauliZ, PauliZ, PauliZ};
         exp_values = Measurer.expval(operations_list, wires_list);
         exp_values_ref = {0.58498357, 0.77015115, 0.91266780};
-        REQUIRE_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
+        CHECK_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
     }
 
     SECTION("Testing list of operators defined by its name:") {
@@ -103,17 +105,17 @@ TEMPLATE_PRODUCT_TEST_CASE("Expected Values", "[Measurements]",
         operations_list = {"PauliX", "PauliX", "PauliX"};
         exp_values = Measurer.expval(operations_list, wires_list);
         exp_values_ref = {0.49272486, 0.42073549, 0.28232124};
-        REQUIRE_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
+        CHECK_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
 
         operations_list = {"PauliY", "PauliY", "PauliY"};
         exp_values = Measurer.expval(operations_list, wires_list);
         exp_values_ref = {-0.64421768, -0.47942553, -0.29552020};
-        REQUIRE_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
+        CHECK_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
 
         operations_list = {"PauliZ", "PauliZ", "PauliZ"};
         exp_values = Measurer.expval(operations_list, wires_list);
         exp_values_ref = {0.58498357, 0.77015115, 0.91266780};
-        REQUIRE_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
+        CHECK_THAT(exp_values, Catch::Approx(exp_values_ref).margin(1e-6));
     }
 }
 
@@ -160,17 +162,17 @@ TEMPLATE_PRODUCT_TEST_CASE("Variances", "[Measurements]",
         operations_list = {PauliX, PauliX, PauliX};
         variances = Measurer.var(operations_list, wires_list);
         variances_ref = {0.7572222, 0.8229816, 0.9202947};
-        REQUIRE_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
+        CHECK_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
 
         operations_list = {PauliY, PauliY, PauliY};
         variances = Measurer.var(operations_list, wires_list);
         variances_ref = {0.5849835, 0.7701511, 0.9126678};
-        REQUIRE_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
+        CHECK_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
 
         operations_list = {PauliZ, PauliZ, PauliZ};
         variances = Measurer.var(operations_list, wires_list);
         variances_ref = {0.6577942, 0.4068672, 0.1670374};
-        REQUIRE_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
+        CHECK_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
     }
 
     SECTION("Testing list of operators defined by its name:") {
@@ -182,17 +184,343 @@ TEMPLATE_PRODUCT_TEST_CASE("Variances", "[Measurements]",
         operations_list = {"PauliX", "PauliX", "PauliX"};
         variances = Measurer.var(operations_list, wires_list);
         variances_ref = {0.7572222, 0.8229816, 0.9202947};
-        REQUIRE_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
+        CHECK_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
 
         operations_list = {"PauliY", "PauliY", "PauliY"};
         variances = Measurer.var(operations_list, wires_list);
         variances_ref = {0.5849835, 0.7701511, 0.9126678};
-        REQUIRE_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
+        CHECK_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
 
         operations_list = {"PauliZ", "PauliZ", "PauliZ"};
         variances = Measurer.var(operations_list, wires_list);
         variances_ref = {0.6577942, 0.4068672, 0.1670374};
-        REQUIRE_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
+        CHECK_THAT(variances, Catch::Approx(variances_ref).margin(1e-6));
+    }
+}
+
+TEMPLATE_PRODUCT_TEST_CASE("Probabilities", "[Measurements]",
+                           (StateVectorLQubitManaged, StateVectorLQubitRaw),
+                           (float, double)) {
+    using StateVectorT = TestType;
+    using PrecisionT = typename StateVectorT::PrecisionT;
+    using ComplexT = typename StateVectorT::ComplexT;
+    SECTION("1 qubit") {
+        auto statevector_data = std::vector<ComplexT>{{1.0, 0.0}, {0.0, 0.0}};
+        StateVectorT statevector(statevector_data.data(),
+                                 statevector_data.size());
+
+        Measurements<StateVectorT> Measurer(statevector);
+
+        SECTION("Testing probs()") {
+            auto p0 = Measurer.probs();
+            statevector.applyOperation("Hadamard", {0}, false);
+            auto p1 = Measurer.probs();
+
+            CHECK_THAT(
+                p0,
+                Catch::Approx(std::vector<PrecisionT>{1.0, 0.0}).margin(1e-7));
+            CHECK_THAT(
+                p1,
+                Catch::Approx(std::vector<PrecisionT>{0.5, 0.5}).margin(1e-7));
+        }
+        SECTION("Testing probs(NamedObs)") {
+            const auto obs1 = Observables::NamedObs<StateVectorT>(
+                {"PauliX"}, std::vector<size_t>{0});
+            const auto obs2 = Observables::NamedObs<StateVectorT>(
+                {"PauliZ"}, std::vector<size_t>{0});
+            const auto obs3 = Observables::NamedObs<StateVectorT>(
+                {"Hadamard"}, std::vector<size_t>{0});
+            auto p0_obs1 = Measurer.probs(obs1);
+            auto p0_obs2 = Measurer.probs(obs2);
+            auto p0_obs3 = Measurer.probs(obs3);
+
+            CHECK_THAT(
+                p0_obs1,
+                Catch::Approx(std::vector<PrecisionT>{0.5, 0.5}).margin(1e-7));
+            CHECK_THAT(
+                p0_obs2,
+                Catch::Approx(std::vector<PrecisionT>{1.0, 0.0}).margin(1e-7));
+            CHECK_THAT(p0_obs3, Catch::Approx(std::vector<PrecisionT>{
+                                                  0.85355339, 0.14644661})
+                                    .margin(1e-7));
+
+            statevector.applyOperation("Hadamard", {0}, false);
+            auto p1_obs1 = Measurer.probs(obs1);
+            auto p1_obs2 = Measurer.probs(obs2);
+            auto p1_obs3 = Measurer.probs(obs3);
+
+            CHECK_THAT(
+                p1_obs1,
+                Catch::Approx(std::vector<PrecisionT>{1.0, 0.0}).margin(1e-7));
+            CHECK_THAT(
+                p1_obs2,
+                Catch::Approx(std::vector<PrecisionT>{0.5, 0.5}).margin(1e-7));
+            CHECK_THAT(p1_obs3, Catch::Approx(std::vector<PrecisionT>{
+                                                  0.85355339, 0.14644661})
+                                    .margin(1e-7));
+
+            statevector.applyOperation("Hadamard", {0}, false);
+            auto p2_obs1 = Measurer.probs(obs1);
+            auto p2_obs2 = Measurer.probs(obs2);
+            auto p2_obs3 = Measurer.probs(obs3);
+
+            CHECK_THAT(
+                p0_obs1,
+                Catch::Approx(std::vector<PrecisionT>{0.5, 0.5}).margin(1e-7));
+            CHECK_THAT(
+                p0_obs2,
+                Catch::Approx(std::vector<PrecisionT>{1.0, 0.0}).margin(1e-7));
+            CHECK_THAT(p0_obs3, Catch::Approx(std::vector<PrecisionT>{
+                                                  0.85355339, 0.14644661})
+                                    .margin(1e-7));
+        }
+    }
+    SECTION("n-qubit") {
+        SECTION("2 qubits") {
+            SECTION("|00> state") {
+                constexpr std::size_t num_qubits = 2;
+                auto statevector_data =
+                    std::vector<ComplexT>((1UL << num_qubits), {0.0, 0.0});
+                const std::vector<std::size_t> wires{0, 1};
+                statevector_data[0] = {1.0, 0.0};
+
+                StateVectorT statevector(statevector_data.data(),
+                                         statevector_data.size());
+                Measurements<StateVectorT> Measurer(statevector);
+
+                auto p0_full = Measurer.probs();
+                auto p0_0 = Measurer.probs({0}, wires);
+                auto p0_1 = Measurer.probs({1}, wires);
+                auto p0_perm0 = Measurer.probs({1, 0}, wires);
+
+                CHECK_THAT(p0_full, Catch::Approx(std::vector<PrecisionT>{
+                                                      1.0, 0.0, 0.0, 0.0})
+                                        .margin(1e-7));
+                CHECK_THAT(p0_0,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+                CHECK_THAT(p0_1,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+
+                CHECK_THAT(p0_perm0, Catch::Approx(std::vector<PrecisionT>{
+                                                       1.0, 0.0, 0.0, 0.0})
+                                         .margin(1e-7));
+
+                statevector.applyOperation("Hadamard", {0}, false);
+                auto p1_full = Measurer.probs();
+                auto p1_0 = Measurer.probs({0}, wires);
+                auto p1_1 = Measurer.probs({1}, wires);
+                auto p1_perm0 = Measurer.probs({1, 0}, wires);
+
+                CHECK_THAT(p1_full, Catch::Approx(std::vector<PrecisionT>{
+                                                      0.5, 0.0, 0.5, 0.0})
+                                        .margin(1e-7));
+                CHECK_THAT(p1_0,
+                           Catch::Approx(std::vector<PrecisionT>{0.5, 0.5})
+                               .margin(1e-7));
+                CHECK_THAT(p1_1,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+                CHECK_THAT(p1_perm0, Catch::Approx(std::vector<PrecisionT>{
+                                                       0.5, 0.5, 0.0, 0.0})
+                                         .margin(1e-7));
+
+                statevector.applyOperation("Hadamard", {1}, false);
+                auto p2_full = Measurer.probs();
+                auto p2_0 = Measurer.probs({0}, wires);
+                auto p2_1 = Measurer.probs({1}, wires);
+                auto p2_perm0 = Measurer.probs({1, 0}, wires);
+
+                CHECK_THAT(p2_full, Catch::Approx(std::vector<PrecisionT>{
+                                                      0.25, 0.25, 0.25, 0.25})
+                                        .margin(1e-7));
+                CHECK_THAT(p2_0,
+                           Catch::Approx(std::vector<PrecisionT>{0.5, 0.5})
+                               .margin(1e-7));
+                CHECK_THAT(p2_1,
+                           Catch::Approx(std::vector<PrecisionT>{0.5, 0.5})
+                               .margin(1e-7));
+                CHECK_THAT(p2_perm0, Catch::Approx(std::vector<PrecisionT>{
+                                                       0.25, 0.25, 0.25, 0.25})
+                                         .margin(1e-7));
+            }
+        }
+        SECTION("3 qubits") {
+            SECTION("|000> state") {
+                constexpr std::size_t num_qubits = 3;
+                auto statevector_data =
+                    std::vector<ComplexT>((1UL << num_qubits), {0.0, 0.0});
+                const std::vector<std::size_t> wires{0, 1, 2};
+                statevector_data[0] = {1.0, 0.0};
+
+                StateVectorT statevector(statevector_data.data(),
+                                         statevector_data.size());
+                Measurements<StateVectorT> Measurer(statevector);
+
+                auto p0_full = Measurer.probs();
+                auto p0_0 = Measurer.probs({0}, wires);
+                auto p0_1 = Measurer.probs({1}, wires);
+                auto p0_2 = Measurer.probs({2}, wires);
+
+                auto p0_01 = Measurer.probs({0, 1}, wires);
+                auto p0_02 = Measurer.probs({0, 2}, wires);
+                auto p0_12 = Measurer.probs({1, 2}, wires);
+
+                auto p0_10 = Measurer.probs({1, 0}, wires);
+                auto p0_20 = Measurer.probs({2, 0}, wires);
+                auto p0_21 = Measurer.probs({2, 1}, wires);
+
+                auto p0_012 = Measurer.probs({0, 1, 2}, wires);
+                auto p0_021 = Measurer.probs({0, 2, 1}, wires);
+                auto p0_102 = Measurer.probs({1, 0, 2}, wires);
+                auto p0_120 = Measurer.probs({1, 2, 0}, wires);
+                auto p0_201 = Measurer.probs({2, 0, 1}, wires);
+                auto p0_210 = Measurer.probs({2, 1, 0}, wires);
+
+                CHECK_THAT(p0_full, Catch::Approx(std::vector<PrecisionT>{
+                                                      1.0, 0.0, 0.0, 0.0, 0.0,
+                                                      0.0, 0.0, 0.0})
+                                        .margin(1e-7));
+                CHECK_THAT(p0_0,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+                CHECK_THAT(p0_1,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+                CHECK_THAT(p0_2,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+
+                CHECK_THAT(p0_01, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p0_02, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p0_12, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+
+                CHECK_THAT(p0_10, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p0_20, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p0_21, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+
+                CHECK_THAT(p0_012, Catch::Approx(std::vector<PrecisionT>{
+                                                     1.0, 0.0, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p0_021, Catch::Approx(std::vector<PrecisionT>{
+                                                     1.0, 0.0, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p0_102, Catch::Approx(std::vector<PrecisionT>{
+                                                     1.0, 0.0, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+
+                CHECK_THAT(p0_120, Catch::Approx(std::vector<PrecisionT>{
+                                                     1.0, 0.0, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p0_201, Catch::Approx(std::vector<PrecisionT>{
+                                                     1.0, 0.0, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p0_210, Catch::Approx(std::vector<PrecisionT>{
+                                                     1.0, 0.0, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+
+                statevector.applyOperation("Hadamard", {0}, false);
+
+                auto p1_full = Measurer.probs();
+                auto p1_0 = Measurer.probs({0}, wires);
+                auto p1_1 = Measurer.probs({1}, wires);
+                auto p1_2 = Measurer.probs({2}, wires);
+
+                auto p1_01 = Measurer.probs({0, 1}, wires);
+                auto p1_02 = Measurer.probs({0, 2}, wires);
+                auto p1_12 = Measurer.probs({1, 2}, wires);
+
+                auto p1_10 = Measurer.probs({1, 0}, wires);
+                auto p1_20 = Measurer.probs({2, 0}, wires);
+                auto p1_21 = Measurer.probs({2, 1}, wires);
+
+                auto p1_012 = Measurer.probs({0, 1, 2}, wires);
+                auto p1_021 = Measurer.probs({0, 2, 1}, wires);
+                auto p1_102 = Measurer.probs({1, 0, 2}, wires);
+                auto p1_120 = Measurer.probs({1, 2, 0}, wires);
+                auto p1_201 = Measurer.probs({2, 0, 1}, wires);
+                auto p1_210 = Measurer.probs({2, 1, 0}, wires);
+
+                CHECK_THAT(p1_full, Catch::Approx(std::vector<PrecisionT>{
+                                                      0.5, 0.0, 0.0, 0.0, 0.5,
+                                                      0.0, 0.0, 0.0})
+                                        .margin(1e-7));
+                CHECK_THAT(p1_0,
+                           Catch::Approx(std::vector<PrecisionT>{0.5, 0.5})
+                               .margin(1e-7));
+                CHECK_THAT(p1_1,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+                CHECK_THAT(p1_2,
+                           Catch::Approx(std::vector<PrecisionT>{1.0, 0.0})
+                               .margin(1e-7));
+
+                CHECK_THAT(p1_01, Catch::Approx(std::vector<PrecisionT>{
+                                                    0.5, 0.0, 0.5, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p1_02, Catch::Approx(std::vector<PrecisionT>{
+                                                    0.5, 0.0, 0.5, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p1_12, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+
+                CHECK_THAT(p1_10, Catch::Approx(std::vector<PrecisionT>{
+                                                    0.5, 0.5, 0.0, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p1_20, Catch::Approx(std::vector<PrecisionT>{
+                                                    0.5, 0.5, 0.0, 0.0})
+                                      .margin(1e-7));
+                CHECK_THAT(p1_21, Catch::Approx(std::vector<PrecisionT>{
+                                                    1.0, 0.0, 0.0, 0.0})
+                                      .margin(1e-7));
+
+                CHECK_THAT(p1_012, Catch::Approx(std::vector<PrecisionT>{
+                                                     0.5, 0.0, 0.0, 0.0, 0.5,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p1_021, Catch::Approx(std::vector<PrecisionT>{
+                                                     0.5, 0.0, 0.0, 0.0, 0.5,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p1_102, Catch::Approx(std::vector<PrecisionT>{
+                                                     0.5, 0.0, 0.5, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+
+                CHECK_THAT(p1_120, Catch::Approx(std::vector<PrecisionT>{
+                                                     0.5, 0.5, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p1_201, Catch::Approx(std::vector<PrecisionT>{
+                                                     0.5, 0.0, 0.5, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+                CHECK_THAT(p1_210, Catch::Approx(std::vector<PrecisionT>{
+                                                     0.5, 0.5, 0.0, 0.0, 0.0,
+                                                     0.0, 0.0, 0.0})
+                                       .margin(1e-7));
+            }
+        }
     }
 }
 
@@ -253,8 +581,8 @@ TEMPLATE_PRODUCT_TEST_CASE("Sample with Metropolis (Local Kernel)",
 
     // compare estimated probabilities to real probabilities
     SECTION("No wires provided:") {
-        REQUIRE_THAT(probabilities,
-                     Catch::Approx(expected_probabilities).margin(.05));
+        CHECK_THAT(probabilities,
+                   Catch::Approx(expected_probabilities).margin(.05));
     }
 }
 
@@ -315,7 +643,7 @@ TEMPLATE_PRODUCT_TEST_CASE("Sample with Metropolis (NonZeroRandom Kernel)",
 
     // compare estimated probabilities to real probabilities
     SECTION("No wires provided:") {
-        REQUIRE_THAT(probabilities,
-                     Catch::Approx(expected_probabilities).margin(.05));
+        CHECK_THAT(probabilities,
+                   Catch::Approx(expected_probabilities).margin(.05));
     }
 }
diff --git a/pennylane_lightning/core/src/utils/NDPermuter.hpp b/pennylane_lightning/core/src/utils/NDPermuter.hpp
new file mode 100644
index 0000000000..46e25c7c56
--- /dev/null
+++ b/pennylane_lightning/core/src/utils/NDPermuter.hpp
@@ -0,0 +1,294 @@
+// Copyright 2024 Xanadu Quantum Technologies Inc.
+
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+
+//     http://www.apache.org/licenses/LICENSE-2.0
+
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// This file is verbatim copied from JET's permuter module at:
+// https://github.com/XanaduAI/jet/tree/v0.2.2/include/jet/permute and reserves
+// all licensing and attributions to that repository's implementations,
+// including inspiration from QFlex https://github.com/ngnrsaa/qflex.
+
+#pragma once
+
+#include <complex>
+#include <numeric>
+#include <set>
+#include <string>
+#include <vector>
+
+#include "Error.hpp"
+
+namespace Pennylane::Util {
+
+/**
+ * @brief Interface for tensor permutation backend.
+ *
+ * The Permuter class represents the front-end interface for calling
+ * permutations, which are a generalization of transposition to high-rank
+ * tensors. The class follows a composition-based approach, where we instantiate
+ * with a given backend permuter, who makes available two `Transpose` methods,
+ * one which returns the transform result, and another which modifies a
+ * reference directly.
+ *
+ *   Example 1:
+ *   const std::vector<size_t> data_in {0,1,2,3,4,5};
+ *   std::vector<size_t> data_out(data_in.size(), 0);
+ *   Permuter<DefaultPermuter<size_t>> p;
+ *   p.Transpose(data_in, {2,3}, data_out, {"a","b"}, {"b","a"});
+ *
+ *   Example 2:
+ *   const std::vector<size_t> data_in {0,1,2,3,4,5};
+ *   Permuter<DefaultPermuter<size_t>> p;
+ *   auto data_out = p.Transpose(data_in, {2,3}, {"a","b"}, {"b","a"});
+ *
+ * @tparam PermuteBackend
+ */
+template <class PermuterBackend> class Permuter {
+  public:
+    /**
+     * @brief Reshape the given lexicographic data vector from old to new index
+     * ordering.
+     *
+     * @tparam T Data participating in the permutation.
+     * @param data_in Input data to be transposed.
+     * @param shape Current shape of the tensor data in each dimension.
+     * @param data_out Output data following the transpose.
+     * @param current_order Current index ordering of the tensor.
+     * @param new_order New index ordering of the tensor.
+     */
+    template <class T>
+    void Transpose(const std::vector<T> &data_in,
+                   const std::vector<size_t> &shape, std::vector<T> &data_out,
+                   const std::vector<std::string> &current_order,
+                   const std::vector<std::string> &new_order) {
+        const std::set<std::string> idx_old(current_order.begin(),
+                                            current_order.end());
+        const std::set<std::string> idx_new(new_order.begin(), new_order.end());
+        const std::size_t data_size =
+            std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<>());
+        PL_ABORT_IF_NOT(idx_old.size() == current_order.size(),
+                        "Duplicate existing indices found. Please ensure "
+                        "indices are unique.");
+        PL_ABORT_IF_NOT(idx_new.size() == new_order.size(),
+                        "Duplicate transpose indices found. Please ensure "
+                        "indices are unique.");
+        PL_ABORT_IF_NOT(shape.size() == new_order.size(),
+                        "Tensor shape does not match number of indices.");
+        PL_ABORT_IF_NOT(data_size == data_in.size(),
+                        "Tensor shape does not match given input tensor data.");
+        PL_ABORT_IF_NOT(
+            data_size == data_out.size(),
+            "Tensor shape does not match given output tensor data.");
+        PL_ABORT_IF_NOT(
+            idx_old == idx_new,
+            "New indices are an invalid permutation of the existing indices");
+
+        permuter_b_.Transpose(data_in, shape, data_out, current_order,
+                              new_order);
+    }
+
+    /**
+     * @brief Reshape the given lexicographic data vector from old to new index
+     * ordering.
+     *
+     * @tparam T Data participating in the permutation.
+     * @param data_in Input data to be transposed.
+     * @param shape Current shape of the tensor data in each dimension.
+     * @param current_order Current index ordering of the tensor.
+     * @param new_order New index ordering of the tensor.
+     * @return std::vector<T> Output data following the transpose.
+     */
+    template <class T>
+    std::vector<T> Transpose(const std::vector<T> &data_in,
+                             const std::vector<size_t> &shape,
+                             const std::vector<std::string> &current_order,
+                             const std::vector<std::string> &new_order) {
+        const std::set<std::string> idx_old(current_order.begin(),
+                                            current_order.end());
+        const std::set<std::string> idx_new(new_order.begin(), new_order.end());
+        const auto data_size =
+            std::accumulate(shape.begin(), shape.end(), 1, std::multiplies<>());
+        PL_ABORT_IF_NOT(idx_old.size() == current_order.size(),
+                        "Duplicate existing indices found. Please ensure "
+                        "indices are unique.");
+        PL_ABORT_IF_NOT(idx_new.size() == new_order.size(),
+                        "Duplicate transpose indices found. Please ensure "
+                        "indices are unique.");
+        PL_ABORT_IF_NOT(shape.size() == new_order.size(),
+                        "Tensor shape does not match number of indices.");
+        PL_ABORT_IF_NOT(data_size == data_in.size(),
+                        "Tensor shape does not match given tensor data.");
+        PL_ABORT_IF_NOT(
+            idx_old == idx_new,
+            "New indices are an invalid permutation of the existing indices");
+
+        PL_ABORT_IF(shape.empty(), "Tensor shape cannot be empty.");
+        PL_ABORT_IF(new_order.empty(), "Tensor indices cannot be empty.");
+        return permuter_b_.Transpose(data_in, shape, current_order, new_order);
+    }
+
+  protected:
+    friend PermuterBackend;
+
+  private:
+    PermuterBackend permuter_b_;
+};
+
+/**
+ * @brief Default Permuter backend class for generalised transforms. Adapted
+ * from QFlex.
+ *
+ * @tparam blocksize Controls the internal data chunk size for cache blocking.
+ */
+template <size_t BLOCKSIZE = 1024> class DefaultPermuter {
+
+  public:
+    /**
+     * @brief Reference-based transpose operation. See `Permuter` class for more
+     * details.
+     */
+    template <class T>
+    void Transpose(const std::vector<T> &data_,
+                   const std::vector<size_t> &shape, std::vector<T> &data_out,
+                   const std::vector<std::string> &old_indices,
+                   const std::vector<std::string> &new_indices) {
+        data_out = data_;
+
+        if (new_indices == old_indices) {
+            return;
+        }
+
+        const std::size_t num_indices = old_indices.size();
+        const std::size_t total_dim = data_.size();
+        std::size_t remaining_data = total_dim;
+
+        if (num_indices == 0) {
+            PL_ABORT("Number of indices cannot be zero.");
+        }
+
+        // Create map_old_to_new_idxpos from old to new indices, and
+        // new_dimensions.
+        std::vector<size_t> map_old_to_new_idxpos(num_indices);
+        std::vector<size_t> new_dimensions(num_indices);
+        for (size_t i = 0; i < num_indices; ++i) {
+            for (size_t j = 0; j < num_indices; ++j) {
+                if (old_indices[i] == new_indices[j]) {
+                    map_old_to_new_idxpos[i] = j;
+                    new_dimensions[j] = shape[i];
+                    break;
+                }
+            }
+        }
+
+        std::vector<size_t> old_super_dimensions(num_indices, 1);
+        std::vector<size_t> new_super_dimensions(num_indices, 1);
+
+        const std::size_t old_dimensions_size = shape.size();
+        for (size_t i = old_dimensions_size; --i;) {
+            old_super_dimensions[i - 1] = old_super_dimensions[i] * shape[i];
+            new_super_dimensions[i - 1] =
+                new_super_dimensions[i] * new_dimensions[i];
+        }
+
+        std::vector<std::size_t> small_map_old_to_new_position(blocksize_);
+
+        // Position old and new.
+        std::size_t po = 0;
+        std::size_t pn = 0;
+        // Counter of the values of each indices in the iteration (old
+        // ordering).
+        std::vector<size_t> old_counter(num_indices, 0);
+        // offset is important when doing this in blocks, as it's indeed
+        // implemented.
+        std::size_t offset = 0;
+        // internal_po keeps track of interations within a block.
+        // Blocks have size `blocksize`.
+        std::size_t internal_po = 0;
+
+        T *data = data_out.data();
+        const T *scratch =
+            data_.data(); // internal pointer offers better performance than
+                          // pointer from argument
+
+        std::size_t effective_max;
+
+        while (true) {
+            // If end of entire opration, break.
+            if (po == total_dim - 1) {
+                break;
+            }
+
+            internal_po = 0;
+            // Each iteration of the while block goes through a new position.
+            // Inside the while, j takes care of increasing indices properly.
+            while (true) {
+                po = 0;
+                pn = 0;
+                for (size_t i = 0; i < num_indices; i++) {
+                    po += old_super_dimensions[i] * old_counter[i];
+                    pn += new_super_dimensions[map_old_to_new_idxpos[i]] *
+                          old_counter[i];
+                }
+                small_map_old_to_new_position[po - offset] = pn;
+
+                bool complete{true};
+                // NOLINTBEGIN
+                for (size_t j = num_indices; j--;) {
+                    if (++old_counter[j] < shape[j]) {
+                        complete = false;
+                        break;
+                    } else {
+                        old_counter[j] = 0;
+                    }
+                }
+                // NOLINTEND
+                // If end of block or end of entire operation, break.
+                if ((++internal_po == blocksize_) || (po == total_dim - 1)) {
+                    break;
+                }
+                // If last index (0) was increased, then go back to fastest
+                // index.
+                if (complete) {
+                    break;
+                }
+            }
+            // Copy data for this block, taking into account offset of
+            // small_map...
+            effective_max = std::min(blocksize_, remaining_data);
+            for (size_t p = 0; p < effective_max; p++) {
+                data[small_map_old_to_new_position[p]] = scratch[offset + p];
+            }
+
+            offset += blocksize_;
+            remaining_data -= blocksize_;
+        }
+    }
+
+    /**
+     * @brief Return-based transpose operation. See `Permuter` class for more
+     * details.
+     */
+    template <class T>
+    std::vector<T> Transpose(std::vector<T> data_,
+                             const std::vector<size_t> &shape,
+                             const std::vector<std::string> &old_indices,
+                             const std::vector<std::string> &new_indices) {
+        std::vector<T> data_out(std::move(data_));
+        Transpose(data_, shape, data_out, old_indices, new_indices);
+        return data_out;
+    }
+
+  private:
+    static constexpr std::size_t blocksize_ = BLOCKSIZE;
+};
+
+} // namespace Pennylane::Util
diff --git a/pennylane_lightning/core/src/utils/tests/CMakeLists.txt b/pennylane_lightning/core/src/utils/tests/CMakeLists.txt
index ccb0c8d062..679dd01fc8 100644
--- a/pennylane_lightning/core/src/utils/tests/CMakeLists.txt
+++ b/pennylane_lightning/core/src/utils/tests/CMakeLists.txt
@@ -29,6 +29,7 @@ target_sources(utils_tests INTERFACE runner_utils.cpp)
 set(TEST_SOURCES    Test_BitUtil.cpp
                     Test_ConstantUtil.cpp
                     Test_Error.cpp
+                    Test_NDPermuter.cpp
                     Test_RuntimeInfo.cpp
                     Test_TypeTraits.cpp
                     Test_Util.cpp
diff --git a/pennylane_lightning/core/src/utils/tests/Test_NDPermuter.cpp b/pennylane_lightning/core/src/utils/tests/Test_NDPermuter.cpp
new file mode 100644
index 0000000000..b906834fe5
--- /dev/null
+++ b/pennylane_lightning/core/src/utils/tests/Test_NDPermuter.cpp
@@ -0,0 +1,35 @@
+
+// Copyright 2018-2023 Xanadu Quantum Technologies Inc.
+
+// Licensed under the Apache License, Version 2.0 (the License);
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+
+// http://www.apache.org/licenses/LICENSE-2.0
+
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an AS IS BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include <cmath>
+#include <complex>
+#include <vector>
+
+#include "NDPermuter.hpp"
+#include "TestHelpers.hpp"
+#include <catch2/catch.hpp>
+
+/// @cond DEV
+namespace {
+using namespace Pennylane;
+using namespace Pennylane::Util;
+} // namespace
+/// @endcond
+
+TEMPLATE_TEST_CASE("Util::DefaultPermuter::Constructibility",
+                   "[Default Constructibility]", DefaultPermuter<>,
+                   DefaultPermuter<8>) {
+    SECTION("DefaultPermuter") { REQUIRE(std::is_constructible_v<TestType>); }
+}