Bugfix/lm multi qubit op

.github/CHANGELOG.md

@@ -56,6 +56,9 @@
 
 ### Bug fixes
 
+* Switch most L-Qubit default kernels to `LM`. Add `LM::multiQubitOp` tests, failing when targeting out-of-order wires clustered close to `num_qubits-1`. Fix the `LM::multiQubitOp` kernel implementation by introducing a generic `revWireParity` routine and replacing the `bitswap`-based implementation. Mimic the changes fixing the corresponding `multiQubitOp` and `expval` functors in L-Kokkos.
+  [(#511)](https://github.com/PennyLaneAI/pennylane-lightning/pull/511)
+
 * Fix RTD builds by removing unsupported `sytem_packages` configuration option.
   [(#491)](https://github.com/PennyLaneAI/pennylane-lightning/pull/491)
 

pennylane_lightning/core/_version.py

@@ -16,4 +16,4 @@
    Version number (major.minor.patch[-label])
 """
 
-__version__ = "0.33.0-dev16"
+__version__ = "0.33.0-dev17"

pennylane_lightning/core/src/simulators/lightning_kokkos/gates/tests/CMakeLists.txt

@@ -20,6 +20,7 @@ target_link_libraries(lightning_kokkos_gates_tests INTERFACE    Catch2::Catch2
                                                                 lightning_kokkos_measurements
                                                                 lightning_kokkos_observables
                                                                 lightning_kokkos
+                                                                lightning_kokkos_utils
                                                                 )
 
 ProcessTestOptions(lightning_kokkos_gates_tests)

pennylane_lightning/core/src/simulators/lightning_kokkos/utils/BitUtilKokkos.hpp

@@ -0,0 +1,72 @@
+// 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.
+
+/**
+ * @file
+ * Defines utility functions for Bitwise operations.
+ */
+#pragma once
+#include <Kokkos_Core.hpp>
+
+#include "BitUtil.hpp"
+
+/// @cond DEV
+namespace {
+using namespace Pennylane::Util;
+using KokkosIntVector = Kokkos::View<std::size_t *>;
+} // namespace
+/// @endcond
+
+namespace Pennylane::LightningKokkos::Util {
+
+constexpr std::size_t one{1};
+
+/**
+ * @brief Compute the parities and shifts for multi-qubit operations.
+ *
+ * @param num_qubits Number of qubits in the state vector.
+ * @param wires List of target wires.
+ * @return std::pair<KokkosIntVector, KokkosIntVector> Parities and shifts for
+ * multi-qubit operations.
+ */
+inline auto wires2Parity(const std::size_t num_qubits,
+                         const std::vector<std::size_t> &wires)
+    -> std::pair<KokkosIntVector, KokkosIntVector> {
+    KokkosIntVector parity;
+    KokkosIntVector rev_wire_shifts;
+
+    std::vector<std::size_t> rev_wires_(wires.size());
+    std::vector<std::size_t> rev_wire_shifts_(wires.size());
+    for (std::size_t k = 0; k < wires.size(); k++) {
+        rev_wires_[k] = (num_qubits - 1) - wires[(wires.size() - 1) - k];
+        rev_wire_shifts_[k] = (one << rev_wires_[k]);
+    }
+    const std::vector<std::size_t> parity_ = revWireParity(rev_wires_);
+
+    Kokkos::View<const size_t *, Kokkos::HostSpace,
+                 Kokkos::MemoryTraits<Kokkos::Unmanaged>>
+        rev_wire_shifts_host(rev_wire_shifts_.data(), rev_wire_shifts_.size());
+    Kokkos::resize(rev_wire_shifts, rev_wire_shifts_host.size());
+    Kokkos::deep_copy(rev_wire_shifts, rev_wire_shifts_host);
+
+    Kokkos::View<const size_t *, Kokkos::HostSpace,
+                 Kokkos::MemoryTraits<Kokkos::Unmanaged>>
+        parity_host(parity_.data(), parity_.size());
+    Kokkos::resize(parity, parity_host.size());
+    Kokkos::deep_copy(parity, parity_host);
+
+    return {parity, rev_wire_shifts};
+}
+
+} // namespace Pennylane::LightningKokkos::Util

pennylane_lightning/core/src/simulators/lightning_qubit/gates/AssignKernelMap_Default.cpp

@@ -121,10 +121,10 @@ void assignKernelsForGateOp_Default() {
     /* Three-qubit gates */
     instance.assignKernelForOp(GateOperation::Toffoli, all_threading,
                                all_memory_model, all_qubit_numbers,
-                               KernelType::PI);
+                               KernelType::LM);
     instance.assignKernelForOp(GateOperation::CSWAP, all_threading,
                                all_memory_model, all_qubit_numbers,
-                               KernelType::PI);
+                               KernelType::LM);
 
     /* QChem gates */
     instance.assignKernelForOp(GateOperation::SingleExcitation, all_threading,
@@ -138,13 +138,13 @@ void assignKernelsForGateOp_Default() {
                                all_qubit_numbers, KernelType::LM);
     instance.assignKernelForOp(GateOperation::DoubleExcitation, all_threading,
                                all_memory_model, all_qubit_numbers,
-                               KernelType::PI);
+                               KernelType::LM);
     instance.assignKernelForOp(GateOperation::DoubleExcitationPlus,
                                all_threading, all_memory_model,
-                               all_qubit_numbers, KernelType::PI);
+                               all_qubit_numbers, KernelType::LM);
     instance.assignKernelForOp(GateOperation::DoubleExcitationMinus,
                                all_threading, all_memory_model,
-                               all_qubit_numbers, KernelType::PI);
+                               all_qubit_numbers, KernelType::LM);
 
     /* Multi-qubit gates */
     instance.assignKernelForOp(GateOperation::MultiRZ, all_threading,
@@ -226,6 +226,6 @@ void assignKernelsForMatrixOp_Default() {
                                KernelType::LM);
     instance.assignKernelForOp(MatrixOperation::MultiQubitOp, all_threading,
                                all_memory_model, all_qubit_numbers,
-                               KernelType::PI);
+                               KernelType::LM);
 }
 } // namespace Pennylane::LightningQubit::KernelMap::Internal

pennylane_lightning/core/src/simulators/lightning_qubit/gates/cpu_kernels/GateImplementationsLM.hpp

@@ -51,79 +51,31 @@
   private:
     /* Alias utility functions */
     static std::pair<size_t, size_t> revWireParity(size_t rev_wire) {
-        using Pennylane::Util::fillLeadingOnes;
-        using Pennylane::Util::fillTrailingOnes;
-
-        const size_t parity_low = fillTrailingOnes(rev_wire);
-        const size_t parity_high = fillLeadingOnes(rev_wire + 1);
-        return {parity_high, parity_low};
+        const auto parity = Pennylane::Util::revWireParity(
+            std::array<std::size_t, 1>{rev_wire});
+        return {parity[1], parity[0]};
     }
-
     static std::tuple<size_t, size_t, size_t> revWireParity(size_t rev_wire0,
                                                             size_t rev_wire1) {
-        using Pennylane::Util::fillLeadingOnes;
-        using Pennylane::Util::fillTrailingOnes;
-
-        const size_t rev_wire_min = std::min(rev_wire0, rev_wire1);
-        const size_t rev_wire_max = std::max(rev_wire0, rev_wire1);
-
-        const size_t parity_low = fillTrailingOnes(rev_wire_min);
-        const size_t parity_high = fillLeadingOnes(rev_wire_max + 1);
-        const size_t parity_middle =
-            fillLeadingOnes(rev_wire_min + 1) & fillTrailingOnes(rev_wire_max);
-
-        return {parity_high, parity_middle, parity_low};
+        const auto parity = Pennylane::Util::revWireParity(
+            std::array<std::size_t, 2>{rev_wire0, rev_wire1});
+        return {parity[2], parity[1], parity[0]};
     }
-
     template <const size_t wire_size = 3>
     static constexpr auto revWireParity(size_t rev_wire0, size_t rev_wire1,
                                         size_t rev_wire2)
         -> std::array<size_t, wire_size + 1> {
-        using Pennylane::Util::fillLeadingOnes;
-        using Pennylane::Util::fillTrailingOnes;
-
-        std::array<size_t, wire_size> rev_wire{rev_wire0, rev_wire1, rev_wire2};
-
-        std::sort(rev_wire.begin(), rev_wire.end());
-
-        const size_t parity_size = rev_wire.size() + 1;
-
-        std::array<size_t, parity_size> parity{};
-
-        parity[0] = fillTrailingOnes(rev_wire[0]);
-        parity[1] =
-            fillLeadingOnes(rev_wire[0] + 1) & fillTrailingOnes(rev_wire[1]);
-        parity[2] =
-            fillLeadingOnes(rev_wire[1] + 1) & fillTrailingOnes(rev_wire[2]);
-        parity[3] = fillLeadingOnes(rev_wire[2] + 1);
-
-        return parity;
+        return Pennylane::Util::revWireParity(
+            std::array<std::size_t, wire_size>{rev_wire0, rev_wire1,
+                                               rev_wire2});
     }
     template <const size_t wire_size = 4>
     static constexpr auto revWireParity(size_t rev_wire0, size_t rev_wire1,
                                         size_t rev_wire2, size_t rev_wire3)
         -> std::array<size_t, wire_size + 1> {
-        using Pennylane::Util::fillLeadingOnes;
-        using Pennylane::Util::fillTrailingOnes;
-
-        std::array<size_t, wire_size> rev_wire{rev_wire0, rev_wire1, rev_wire2,
-                                               rev_wire3};
-
-        std::sort(rev_wire.begin(), rev_wire.end());
-
-        const size_t parity_size = rev_wire.size() + 1;
-        std::array<size_t, parity_size> parity{};
-
-        parity[0] = fillTrailingOnes(rev_wire[0]);
-        parity[1] =
-            fillLeadingOnes(rev_wire[0] + 1) & fillTrailingOnes(rev_wire[1]);
-        parity[2] =
-            fillLeadingOnes(rev_wire[1] + 1) & fillTrailingOnes(rev_wire[2]);
-        parity[3] =
-            fillLeadingOnes(rev_wire[2] + 1) & fillTrailingOnes(rev_wire[3]);
-        parity[4] = fillLeadingOnes(rev_wire[3] + 1);
-
-        return parity;
+        return Pennylane::Util::revWireParity(
+            std::array<std::size_t, wire_size>{rev_wire0, rev_wire1, rev_wire2,
+                                               rev_wire3});
     }
 
   public:
@@ -241,6 +193,7 @@
             }
         }
     }
+
     /**
      * @brief Apply a two qubit gate to the statevector.
      *
@@ -328,29 +281,45 @@
 
     template <class PrecisionT>
     static void
-    applyMultiQubitOp(std::complex<PrecisionT> *arr, size_t num_qubits,
+    applyMultiQubitOp(std::complex<PrecisionT> *arr, std::size_t num_qubits,
                       const std::complex<PrecisionT> *matrix,
-                      const std::vector<size_t> &wires, bool inverse) {
+                      const std::vector<std::size_t> &wires, bool inverse) {
         using Pennylane::Util::bitswap;
+        constexpr std::size_t one{1};
         PL_ASSERT(num_qubits >= wires.size());
 
-        const size_t dim = static_cast<size_t>(1U) << wires.size();
-        std::vector<size_t> indices(dim);
+        const std::size_t dim = one << wires.size();
+        std::vector<std::size_t> indices(dim);
         std::vector<std::complex<PrecisionT>> coeffs_in(dim, 0.0);
 
+        std::vector<std::size_t> rev_wires(wires.size());
+        std::vector<std::size_t> rev_wire_shifts(wires.size());
+        for (std::size_t k = 0; k < wires.size(); k++) {
+            rev_wires[k] = (num_qubits - 1) - wires[(wires.size() - 1) - k];
+            rev_wire_shifts[k] = (one << rev_wires[k]);
+        }
+        const std::vector<std::size_t> parity =
+            Pennylane::Util::revWireParity(rev_wires);
+        PL_ASSERT(wires.size() == parity.size() - 1);
+
         if (inverse) {
-            for (size_t k = 0; k < exp2(num_qubits); k += dim) {
-                for (size_t inner_idx = 0; inner_idx < dim; inner_idx++) {
-                    size_t idx = k | inner_idx;
-                    const size_t n_wires = wires.size();
-                    for (size_t pos = 0; pos < n_wires; pos++) {
-                        idx = bitswap(idx, n_wires - pos - 1,
-                                      num_qubits - wires[pos] - 1);
+            for (std::size_t k = 0; k < exp2(num_qubits - wires.size()); k++) {
+                std::size_t idx = (k & parity[0]);
+                for (std::size_t i = 1; i < parity.size(); i++) {
+                    idx |= ((k << i) & parity[i]);
+                }
+                indices[0] = idx;
+                coeffs_in[0] = arr[idx];
+                for (std::size_t inner_idx = 1; inner_idx < dim; inner_idx++) {
+                    idx = indices[0];
+                    for (std::size_t i = 0; i < wires.size(); i++) {
+                        if ((inner_idx & (one << i)) != 0) {
+                            idx |= rev_wire_shifts[i];
+                        }
                     }
                     indices[inner_idx] = idx;
                     coeffs_in[inner_idx] = arr[idx];
                 }
-
                 for (size_t i = 0; i < dim; i++) {
                     const auto idx = indices[i];
                     arr[idx] = 0.0;
@@ -363,25 +332,29 @@
                 }
             }
         } else {
-            for (size_t k = 0; k < exp2(num_qubits); k += dim) {
-                for (size_t inner_idx = 0; inner_idx < dim; inner_idx++) {
-                    size_t idx = k | inner_idx;
-                    const size_t n_wires = wires.size();
-                    for (size_t pos = 0; pos < n_wires; pos++) {
-                        idx = bitswap(idx, n_wires - pos - 1,
-                                      num_qubits - wires[pos] - 1);
+            for (std::size_t k = 0; k < exp2(num_qubits - wires.size()); k++) {
+                std::size_t idx = (k & parity[0]);
+                for (std::size_t i = 1; i < parity.size(); i++) {
+                    idx |= ((k << i) & parity[i]);
+                }
+                indices[0] = idx;
+                coeffs_in[0] = arr[idx];
+                for (std::size_t inner_idx = 1; inner_idx < dim; inner_idx++) {
+                    idx = indices[0];
+                    for (std::size_t i = 0; i < wires.size(); i++) {
+                        if ((inner_idx & (one << i)) != 0) {
+                            idx |= rev_wire_shifts[i];
+                        }
                     }
                     indices[inner_idx] = idx;
                     coeffs_in[inner_idx] = arr[idx];
                 }
-
-                for (size_t i = 0; i < dim; i++) {
-                    const auto idx = indices[i];
-                    arr[idx] = 0.0;
-                    const size_t base_idx = i * dim;
-
-                    for (size_t j = 0; j < dim; j++) {
-                        arr[idx] += matrix[base_idx + j] * coeffs_in[j];
+                for (std::size_t i = 0; i < dim; i++) {
+                    const auto index = indices[i];
+                    arr[index] = 0.0;
+                    const std::size_t base_idx = i * dim;
+                    for (std::size_t j = 0; j < dim; j++) {
+                        arr[index] += matrix[base_idx + j] * coeffs_in[j];
                     }
                 }
             }
@@ -390,86 +363,86 @@
 
     template <class PrecisionT>
     static void applyIdentity(std::complex<PrecisionT> *arr,
                               const size_t num_qubits,
                               const std::vector<size_t> &wires,
                               [[maybe_unused]] bool inverse) {
         PL_ASSERT(wires.size() == 1);
         static_cast<void>(arr);        // No-op
         static_cast<void>(num_qubits); // No-op
         static_cast<void>(wires);      // No-op
     }
 
     template <class PrecisionT>
     static void applyPauliX(std::complex<PrecisionT> *arr,
                             const size_t num_qubits,
                             const std::vector<size_t> &wires,
                             [[maybe_unused]] bool inverse) {
         PL_ASSERT(wires.size() == 1);
 
         const size_t rev_wire = num_qubits - wires[0] - 1;
         const size_t rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
 
         const auto [parity_high, parity_low] = revWireParity(rev_wire);
 
         for (size_t k = 0; k < exp2(num_qubits - 1); k++) {
             const size_t i0 = ((k << 1U) & parity_high) | (parity_low & k);
             const size_t i1 = i0 | rev_wire_shift;
             std::swap(arr[i0], arr[i1]);
         }
     }
 
     template <class PrecisionT>
     static void applyPauliY(std::complex<PrecisionT> *arr,
                             const size_t num_qubits,
                             const std::vector<size_t> &wires,
                             [[maybe_unused]] bool inverse) {
         PL_ASSERT(wires.size() == 1);
         const size_t rev_wire = num_qubits - wires[0] - 1;
         const size_t rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
 
         const auto [parity_high, parity_low] = revWireParity(rev_wire);
 
         for (size_t k = 0; k < exp2(num_qubits - 1); k++) {
             const size_t i0 = ((k << 1U) & parity_high) | (parity_low & k);
             const size_t i1 = i0 | rev_wire_shift;
             const auto v0 = arr[i0];
             const auto v1 = arr[i1];
             arr[i0] = {std::imag(v1), -std::real(v1)};
             arr[i1] = {-std::imag(v0), std::real(v0)};
         }
     }
 
     template <class PrecisionT>
     static void applyPauliZ(std::complex<PrecisionT> *arr,
                             const size_t num_qubits,
                             const std::vector<size_t> &wires,
                             [[maybe_unused]] bool inverse) {
         PL_ASSERT(wires.size() == 1);
         const size_t rev_wire = num_qubits - wires[0] - 1;
         const size_t rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
         const auto [parity_high, parity_low] = revWireParity(rev_wire);
 
         for (size_t k = 0; k < exp2(num_qubits - 1); k++) {
             const size_t i0 = ((k << 1U) & parity_high) | (parity_low & k);
             const size_t i1 = i0 | rev_wire_shift;
             arr[i1] *= -1;
         }
     }
 
     template <class PrecisionT>
     static void applyHadamard(std::complex<PrecisionT> *arr,
                               const size_t num_qubits,
                               const std::vector<size_t> &wires,
                               [[maybe_unused]] bool inverse) {
         PL_ASSERT(wires.size() == 1);
         constexpr static auto isqrt2 = INVSQRT2<PrecisionT>();
         const size_t rev_wire = num_qubits - wires[0] - 1;
         const size_t rev_wire_shift = (static_cast<size_t>(1U) << rev_wire);
         const auto [parity_high, parity_low] = revWireParity(rev_wire);
 
         for (size_t k = 0; k < exp2(num_qubits - 1); k++) {
             const size_t i0 = ((k << 1U) & parity_high) | (parity_low & k);
             const size_t i1 = i0 | rev_wire_shift;
             const std::complex<PrecisionT> v0 = arr[i0];
             const std::complex<PrecisionT> v1 = arr[i1];
             arr[i0] = isqrt2 * v0 + isqrt2 * v1;

pennylane_lightning/core/src/simulators/lightning_qubit/gates/tests/Test_KernelMap.cpp

@@ -110,7 +110,7 @@ TEST_CASE("Test several limiting cases of default kernels", "[KernelMap]") {
         OperationKernelMap<Pennylane::Gates::GateOperation>::getInstance();
     SECTION("Single thread, large number of qubits") {
         // For large N, single thread calls "LM" for all single- and two-qubit
-        // gates. For k-qubit gates with k >= 3, we use PI.
+        // gates.
         auto gate_map = instance.getKernelMap(28, Threading::SingleThread,
                                               CPUMemoryModel::Unaligned);
         for_each_enum<Pennylane::Gates::GateOperation>(
@@ -123,9 +123,6 @@ TEST_CASE("Test several limiting cases of default kernels", "[KernelMap]") {
                                   gate_op) <= 2) {
                     REQUIRE(gate_map[gate_op] ==
                             Pennylane::Gates::KernelType::LM);
-                } else {
-                    REQUIRE(gate_map[gate_op] ==
-                            Pennylane::Gates::KernelType::PI);
                 }
             });
     }