diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md index 59f87d0c95..2dd314c424 100644 --- a/.github/CHANGELOG.md +++ b/.github/CHANGELOG.md @@ -2,6 +2,12 @@ ### New features since last release +* Add `SparseHamiltonian` support for Lightning-GPU. + [(#526)] (https://github.com/PennyLaneAI/pennylane-lightning/pull/526) + +* Add `SparseHamiltonian` support for Lightning-Kokkos. + [(#527)] (https://github.com/PennyLaneAI/pennylane-lightning/pull/527) + * Integrate python/pybind layer of distributed Lightning-GPU into the Lightning monorepo with python unit tests. [(#518)] (https://github.com/PennyLaneAI/pennylane-lightning/pull/518) diff --git a/.github/workflows/format.yml b/.github/workflows/format.yml index 6cbde05023..350312d253 100644 --- a/.github/workflows/format.yml +++ b/.github/workflows/format.yml @@ -84,7 +84,7 @@ jobs: cp -rf ${{ github.workspace}}/Kokkos_install/${{ matrix.exec_model }}/* Kokkos/ - name: Install dependencies - run: sudo apt update && sudo apt -y install clang-tidy-14 cmake g++-10 ninja-build libomp-14-dev + run: sudo apt update && sudo apt -y install clang-tidy-14 cmake gcc-11 g++-11 ninja-build libomp-14-dev env: DEBIAN_FRONTEND: noninteractive @@ -96,5 +96,6 @@ jobs: -DBUILD_TESTS=ON \ -DENABLE_WARNINGS=ON \ -DPL_BACKEND=${{ matrix.pl_backend }} \ - -DCMAKE_CXX_COMPILER="$(which g++-10)" + -DCMAKE_CXX_COMPILER="$(which g++-11)" \ + -DCMAKE_C_COMPILER="$(which gcc-11)" cmake --build ./Build \ No newline at end of file diff --git a/.github/workflows/tests_gpu.yml b/.github/workflows/tests_gpu_kokkos.yml similarity index 99% rename from .github/workflows/tests_gpu.yml rename to .github/workflows/tests_gpu_kokkos.yml index bbf240d91c..40b9462f03 100644 --- a/.github/workflows/tests_gpu.yml +++ b/.github/workflows/tests_gpu_kokkos.yml @@ -1,4 +1,4 @@ -name: Testing (GPU) +name: Testing::LKokkos::GPU on: pull_request: push: diff --git a/.github/workflows/tests_linux_x86_mpi_gpu.yml b/.github/workflows/tests_linux_x86_mpi_gpu.yml index e43e9ca241..03d88cba06 100644 --- a/.github/workflows/tests_linux_x86_mpi_gpu.yml +++ b/.github/workflows/tests_linux_x86_mpi_gpu.yml @@ -246,15 +246,16 @@ jobs: run: | source /etc/profile.d/modules.sh && module use /opt/modules/ && module load ${{ matrix.mpilib }} PL_DEVICE=lightning.gpu python -m pytest ./tests/ $COVERAGE_FLAGS + mv coverage.xml coverage-${{ github.job }}-lightning_gpu_${{ matrix.mpilib }}-0.xml PL_DEVICE=lightning.gpu /opt/mpi/${{ matrix.mpilib }}/bin/mpirun -np 2 python -m pytest ./mpitests $COVERAGE_FLAGS + mv coverage.xml coverage-${{ github.job }}-lightning_gpu_${{ matrix.mpilib }}-1.xml # PL_DEVICE=lightning.gpu /opt/mpi/${{ matrix.mpilib }}/bin/mpirun --oversubscribe -n 4 pytest -s -x mpitests/test_device.py -k test_create_device $COVERAGE_FLAGS - mv coverage.xml coverage-${{ github.job }}-lightning_gpu_${{ matrix.mpilib }}.xml - name: Upload code coverage results uses: actions/upload-artifact@v3 with: name: ubuntu-codecov-results-python - path: coverage-${{ github.job }}-lightning_gpu_${{ matrix.mpilib }}.xml + path: coverage-${{ github.job }}-lightning_gpu_${{ matrix.mpilib }}-*.xml if-no-files-found: error - name: Cleanup diff --git a/pennylane_lightning/core/_serialize.py b/pennylane_lightning/core/_serialize.py index a848f0a3e3..12a17a0cc5 100644 --- a/pennylane_lightning/core/_serialize.py +++ b/pennylane_lightning/core/_serialize.py @@ -25,6 +25,8 @@ Identity, StatePrep, Rot, + Hamiltonian, + SparseHamiltonian, ) from pennylane.operation import Tensor from pennylane.tape import QuantumTape @@ -52,6 +54,7 @@ class QuantumScriptSerializer: # pylint: disable=import-outside-toplevel, too-many-instance-attributes, c-extension-no-member def __init__(self, device_name, use_csingle: bool = False, use_mpi: bool = False): self.use_csingle = use_csingle + self.device_name = device_name if device_name == "lightning.qubit": try: import pennylane_lightning.lightning_qubit_ops as lightning_ops @@ -75,6 +78,7 @@ def __init__(self, device_name, use_csingle: bool = False, use_mpi: bool = False ) from exception else: raise DeviceError(f'The device name "{device_name}" is not a valid option.') + self.statevector_c64 = lightning_ops.StateVectorC64 self.statevector_c128 = lightning_ops.StateVectorC128 self.named_obs_c64 = lightning_ops.observables.NamedObsC64 self.named_obs_c128 = lightning_ops.observables.NamedObsC128 @@ -84,20 +88,25 @@ def __init__(self, device_name, use_csingle: bool = False, use_mpi: bool = False self.tensor_prod_obs_c128 = lightning_ops.observables.TensorProdObsC128 self.hamiltonian_c64 = lightning_ops.observables.HamiltonianC64 self.hamiltonian_c128 = lightning_ops.observables.HamiltonianC128 + self.sparse_hamiltonian_c64 = lightning_ops.observables.SparseHamiltonianC64 + self.sparse_hamiltonian_c128 = lightning_ops.observables.SparseHamiltonianC128 self._use_mpi = False if use_mpi: self._use_mpi = use_mpi - self.statevectormpi_c128 = lightning_ops.StateVectorMPIC128 - self.named_obsmpi_c64 = lightning_ops.observablesMPI.NamedObsMPIC64 - self.named_obsmpi_c128 = lightning_ops.observablesMPI.NamedObsMPIC128 - self.hermitian_obsmpi_c64 = lightning_ops.observablesMPI.HermitianObsMPIC64 - self.hermitian_obsmpi_c128 = lightning_ops.observablesMPI.HermitianObsMPIC128 - self.tensor_prod_obsmpi_c64 = lightning_ops.observablesMPI.TensorProdObsMPIC64 - self.tensor_prod_obsmpi_c128 = lightning_ops.observablesMPI.TensorProdObsMPIC128 - self.hamiltonianmpi_c64 = lightning_ops.observablesMPI.HamiltonianMPIC64 - self.hamiltonianmpi_c128 = lightning_ops.observablesMPI.HamiltonianMPIC128 + self.statevector_mpi_c64 = lightning_ops.StateVectorMPIC64 + self.statevector_mpi_c128 = lightning_ops.StateVectorMPIC128 + self.named_obs_mpi_c64 = lightning_ops.observablesMPI.NamedObsMPIC64 + self.named_obs_mpi_c128 = lightning_ops.observablesMPI.NamedObsMPIC128 + self.hermitian_obs_mpi_c64 = lightning_ops.observablesMPI.HermitianObsMPIC64 + self.hermitian_obs_mpi_c128 = lightning_ops.observablesMPI.HermitianObsMPIC128 + self.tensor_prod_obs_mpi_c64 = lightning_ops.observablesMPI.TensorProdObsMPIC64 + self.tensor_prod_obs_mpi_c128 = lightning_ops.observablesMPI.TensorProdObsMPIC128 + self.hamiltonian_mpi_c64 = lightning_ops.observablesMPI.HamiltonianMPIC64 + self.hamiltonian_mpi_c128 = lightning_ops.observablesMPI.HamiltonianMPIC128 + + self._mpi_manager = lightning_ops.MPIManager @property def ctype(self): @@ -113,37 +122,44 @@ def rtype(self): def sv_type(self): """State vector matching ``use_csingle`` precision (and MPI if it is supported).""" if self._use_mpi: - return self.statevectormpi_c128 - return self.statevector_c128 + return self.statevector_mpi_c64 if self.use_csingle else self.statevector_mpi_c128 + return self.statevector_c64 if self.use_csingle else self.statevector_c128 @property def named_obs(self): """Named observable matching ``use_csingle`` precision.""" if self._use_mpi: - return self.named_obsmpi_c64 if self.use_csingle else self.named_obsmpi_c128 + return self.named_obs_mpi_c64 if self.use_csingle else self.named_obs_mpi_c128 return self.named_obs_c64 if self.use_csingle else self.named_obs_c128 @property def hermitian_obs(self): """Hermitian observable matching ``use_csingle`` precision.""" if self._use_mpi: - return self.hermitian_obsmpi_c64 if self.use_csingle else self.hermitian_obsmpi_c128 + return self.hermitian_obs_mpi_c64 if self.use_csingle else self.hermitian_obs_mpi_c128 return self.hermitian_obs_c64 if self.use_csingle else self.hermitian_obs_c128 @property def tensor_obs(self): """Tensor product observable matching ``use_csingle`` precision.""" if self._use_mpi: - return self.tensor_prod_obsmpi_c64 if self.use_csingle else self.tensor_prod_obsmpi_c128 + return ( + self.tensor_prod_obs_mpi_c64 if self.use_csingle else self.tensor_prod_obs_mpi_c128 + ) return self.tensor_prod_obs_c64 if self.use_csingle else self.tensor_prod_obs_c128 @property def hamiltonian_obs(self): """Hamiltonian observable matching ``use_csingle`` precision.""" if self._use_mpi: - return self.hamiltonianmpi_c64 if self.use_csingle else self.hamiltonianmpi_c128 + return self.hamiltonian_mpi_c64 if self.use_csingle else self.hamiltonian_mpi_c128 return self.hamiltonian_c64 if self.use_csingle else self.hamiltonian_c128 + @property + def sparse_hamiltonian_obs(self): + """SparseHamiltonian observable matching ``use_csingle`` precision.""" + return self.sparse_hamiltonian_c64 if self.use_csingle else self.sparse_hamiltonian_c128 + def _named_obs(self, observable, wires_map: dict): """Serializes a Named observable""" wires = [wires_map[w] for w in observable.wires] @@ -168,6 +184,37 @@ def _hamiltonian(self, observable, wires_map: dict): terms = [self._ob(t, wires_map) for t in observable.ops] return self.hamiltonian_obs(coeffs, terms) + def _sparse_hamiltonian(self, observable, wires_map: dict): + """Serialize an observable (Sparse Hamiltonian) + + Args: + observable (Observable): the input observable (Sparse Hamiltonian) + wire_map (dict): a dictionary mapping input wires to the device's backend wires + + Returns: + sparse_hamiltonian_obs (SparseHamiltonianC64 or SparseHamiltonianC128): A Sparse Hamiltonian observable object compatible with the C++ backend + """ + + if self._use_mpi: + Hmat = Hamiltonian([1.0], [Identity(0)]).sparse_matrix() + H_sparse = SparseHamiltonian(Hmat, wires=range(1)) + spm = H_sparse.sparse_matrix() + # Only root 0 needs the overall sparsematrix data + if self._mpi_manager().getRank() == 0: + spm = observable.sparse_matrix() + self._mpi_manager().Barrier() + else: + spm = observable.sparse_matrix() + data = np.array(spm.data).astype(self.ctype) + indices = np.array(spm.indices).astype(np.int64) + offsets = np.array(spm.indptr).astype(np.int64) + + wires = [] + wires_list = observable.wires.tolist() + wires.extend([wires_map[w] for w in wires_list]) + + return self.sparse_hamiltonian_obs(data, indices, offsets, wires) + def _pauli_word(self, observable, wires_map: dict): """Serialize a :class:`pennylane.pauli.PauliWord` into a Named or Tensor observable.""" if len(observable) == 1: @@ -195,6 +242,8 @@ def _ob(self, observable, wires_map): return self._tensor_ob(observable, wires_map) if observable.name == "Hamiltonian": return self._hamiltonian(observable, wires_map) + if observable.name == "SparseHamiltonian": + return self._sparse_hamiltonian(observable, wires_map) if isinstance(observable, (PauliX, PauliY, PauliZ, Identity, Hadamard)): return self._named_obs(observable, wires_map) if observable._pauli_rep is not None: diff --git a/pennylane_lightning/core/_version.py b/pennylane_lightning/core/_version.py index b31d7d35ba..d276e92d29 100644 --- a/pennylane_lightning/core/_version.py +++ b/pennylane_lightning/core/_version.py @@ -16,4 +16,4 @@ Version number (major.minor.patch[-label]) """ -__version__ = "0.33.0-dev23" +__version__ = "0.33.0-dev24" diff --git a/pennylane_lightning/core/src/bindings/Bindings.hpp b/pennylane_lightning/core/src/bindings/Bindings.hpp index 30c994a719..1368ee3a94 100644 --- a/pennylane_lightning/core/src/bindings/Bindings.hpp +++ b/pennylane_lightning/core/src/bindings/Bindings.hpp @@ -287,7 +287,8 @@ void registerInfo(py::module_ &m) { * @tparam StateVectorT * @param m Pybind module */ -template void registerObservables(py::module_ &m) { +template +void registerBackendAgnosticObservables(py::module_ &m) { using PrecisionT = typename StateVectorT::PrecisionT; // Statevector's precision. using ComplexT = @@ -627,7 +628,8 @@ template void lightningClassBindings(py::module_ &m) { /* Observables submodule */ py::module_ obs_submodule = m.def_submodule("observables", "Submodule for observables classes."); - registerObservables(obs_submodule); + registerBackendAgnosticObservables(obs_submodule); + registerBackendSpecificObservables(obs_submodule); //***********************************************************************// // Measurements diff --git a/pennylane_lightning/core/src/bindings/BindingsMPI.hpp b/pennylane_lightning/core/src/bindings/BindingsMPI.hpp index 6c828b0ff3..41276afe5d 100644 --- a/pennylane_lightning/core/src/bindings/BindingsMPI.hpp +++ b/pennylane_lightning/core/src/bindings/BindingsMPI.hpp @@ -82,6 +82,10 @@ template void registerObservablesMPI(py::module_ &m) { using np_arr_c = py::array_t, py::array::c_style>; using np_arr_r = py::array_t; + using np_arr_sparse_ind = typename std::conditional< + std::is_same::value, + py::array_t, + py::array_t>::type; std::string class_name; @@ -191,6 +195,48 @@ template void registerObservablesMPI(py::module_ &m) { return self == other_cast; }, "Compare two observables"); +#ifdef _ENABLE_PLGPU + class_name = "SparseHamiltonianMPIC" + bitsize; + using SpIDX = typename SparseHamiltonianMPI::IdxT; + py::class_, + std::shared_ptr>, + Observable>(m, class_name.c_str(), + py::module_local()) + .def(py::init([](const np_arr_c &data, const np_arr_sparse_ind &indices, + const np_arr_sparse_ind &offsets, + const std::vector &wires) { + const py::buffer_info buffer_data = data.request(); + const auto *data_ptr = static_cast(buffer_data.ptr); + + const py::buffer_info buffer_indices = indices.request(); + const auto *indices_ptr = static_cast(buffer_indices.ptr); + + const py::buffer_info buffer_offsets = offsets.request(); + const auto *offsets_ptr = static_cast(buffer_offsets.ptr); + + return SparseHamiltonianMPI{ + std::vector({data_ptr, data_ptr + data.size()}), + std::vector({indices_ptr, indices_ptr + indices.size()}), + std::vector({offsets_ptr, offsets_ptr + offsets.size()}), + wires}; + })) + .def("__repr__", &SparseHamiltonianMPI::getObsName) + .def("get_wires", &SparseHamiltonianMPI::getWires, + "Get wires of observables") + .def( + "__eq__", + [](const SparseHamiltonianMPI &self, + py::handle other) -> bool { + if (!py::isinstance>( + other)) { + return false; + } + auto other_cast = + other.cast>(); + return self == other_cast; + }, + "Compare two observables"); +#endif } /** diff --git a/pennylane_lightning/core/src/observables/Observables.hpp b/pennylane_lightning/core/src/observables/Observables.hpp index 42227f62e5..183a9fcd82 100644 --- a/pennylane_lightning/core/src/observables/Observables.hpp +++ b/pennylane_lightning/core/src/observables/Observables.hpp @@ -414,4 +414,112 @@ class HamiltonianBase : public Observable { } }; +/** + * @brief Sparse representation of SparseHamiltonian + * + * @tparam T Floating-point precision. + */ +template +class SparseHamiltonianBase : public Observable { + public: + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; +#ifdef _ENABLE_PLGPU + using IdxT = + typename std::conditional::value, + int32_t, int64_t>::type; +#else + using IdxT = std::size_t; +#endif + + protected: + std::vector data_; + std::vector indices_; + std::vector offsets_; + std::vector wires_; + + private: + [[nodiscard]] bool + isEqual(const Observable &other) const override { + const auto &other_cast = + static_cast &>(other); + return data_ == other_cast.data_ && indices_ == other_cast.indices_ && + offsets_ == other_cast.offsets_ && (wires_ == other_cast.wires_); + } + + public: + /** + * @brief Create a SparseHamiltonianBase from data, indices and offsets in + * CSR format. + * + * @param data Arguments to construct data + * @param indices Arguments to construct indices + * @param offsets Arguments to construct offsets + * @param wires Arguments to construct wires + */ + template > + SparseHamiltonianBase(T1 &&data, T2 &&indices, T3 &&offsets, T4 &&wires) + : data_{std::forward(data)}, indices_{std::forward(indices)}, + offsets_{std::forward(offsets)}, wires_{std::forward(wires)} { + PL_ASSERT(data_.size() == indices_.size()); + } + + /** + * @brief Convenient wrapper for the constructor as the constructor does not + * convert the std::shared_ptr with a derived class correctly. + * + * This function is useful as std::make_shared does not handle + * brace-enclosed initializer list correctly. + * + * @param data Argument to construct data + * @param indices Argument to construct indices + * @param offsets Argument to construct offsets + * @param wires Argument to construct wires + */ + static auto create(std::initializer_list data, + std::initializer_list indices, + std::initializer_list offsets, + std::initializer_list wires) + -> std::shared_ptr> { + // NOLINTBEGIN(*-move-const-arg) + return std::shared_ptr>( + new SparseHamiltonianBase{ + std::move(data), std::move(indices), std::move(offsets), + std::move(wires)}); + // NOLINTEND(*-move-const-arg) + } + + void applyInPlace([[maybe_unused]] StateVectorT &sv) const override { + PL_ABORT("For SparseHamiltonian Observables, the applyInPlace method " + "must be " + "defined at the backend level."); + } + + [[nodiscard]] auto getObsName() const -> std::string override { + using Pennylane::Util::operator<<; + std::ostringstream ss; + ss << "SparseHamiltonian: {\n'data' : \n"; + for (const auto &d : data_) { + ss << "{" << d.real() << ", " << d.imag() << "}, "; + } + ss << ",\n'indices' : \n"; + for (const auto &i : indices_) { + ss << i << ", "; + } + ss << ",\n'offsets' : \n"; + for (const auto &o : offsets_) { + ss << o << ", "; + } + ss << "\n}"; + return ss.str(); + } + /** + * @brief Get the wires the observable applies to. + */ + [[nodiscard]] auto getWires() const -> std::vector override { + return wires_; + }; +}; + } // namespace Pennylane::Observables \ No newline at end of file diff --git a/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp b/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp index 2c4e3b60da..b0541bd3e1 100644 --- a/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp +++ b/pennylane_lightning/core/src/observables/tests/Test_Observables.cpp @@ -30,8 +30,8 @@ /// @cond DEV namespace { using namespace Pennylane::Observables; - using Pennylane::Util::createProductState; +using Pennylane::Util::createRandomStateVectorData; using Pennylane::Util::createZeroState; using Pennylane::Util::isApproxEqual; using Pennylane::Util::LightningException; @@ -464,3 +464,85 @@ TEST_CASE("Methods implemented in the HamiltonianBase class", testHamiltonianBase(); } } + +template void testSparseHamiltonianBase() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + + const std::size_t num_qubits = 3; + std::mt19937 re{1337}; + + auto sparseH = SparseHamiltonianBase::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {7, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, {0, 1, 2}); + + DYNAMIC_SECTION("SparseHamiltonianBase - isEqual - " + << StateVectorToName::name) { + auto sparseH0 = SparseHamiltonianBase::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {7, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, + {0, 1, 2}); + auto sparseH1 = SparseHamiltonianBase::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {7, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, + {0, 1, 2}); + auto sparseH2 = SparseHamiltonianBase::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {8, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, + {0, 1, 2}); + + REQUIRE(*sparseH0 == *sparseH1); + REQUIRE(*sparseH0 != *sparseH2); + } + + DYNAMIC_SECTION("SparseHamiltonianBase - getWires - " + << StateVectorToName::name) { + REQUIRE(sparseH->getWires() == std::vector{0, 1, 2}); + } + + DYNAMIC_SECTION("SparseHamiltonianBase - getObsName - " + << StateVectorToName::name) { + REQUIRE(sparseH->getObsName() == + "SparseHamiltonian: {\n" + "'data' : \n" + "{1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 0}, " + "{1, 0}, ,\n" + "'indices' : \n" + "7, 6, 5, 4, 3, 2, 1, 0, ,\n" + "'offsets' : \n" + "0, 1, 2, 3, 4, 5, 6, 7, 8, \n" + "}"); + } + + DYNAMIC_SECTION("SparseHamiltonianBase - applyInPlace must fail - " + << StateVectorToName::name) { + + auto init_state = + createRandomStateVectorData(re, num_qubits); + + StateVectorT state_vector(init_state.data(), init_state.size()); + + REQUIRE_THROWS_AS(sparseH->applyInPlace(state_vector), + LightningException); + } + + testSparseHamiltonianBase(); + } +} + +TEST_CASE("Methods implemented in the SparseHamiltonianBase class", + "[SparseHamiltonianBase]") { + if constexpr (BACKEND_FOUND) { + testSparseHamiltonianBase(); + } +} \ No newline at end of file diff --git a/pennylane_lightning/core/src/observables/tests/mpi/Test_ObservablesMPI.cpp b/pennylane_lightning/core/src/observables/tests/mpi/Test_ObservablesMPI.cpp index fdcfa8c2ea..ca12e4272f 100644 --- a/pennylane_lightning/core/src/observables/tests/mpi/Test_ObservablesMPI.cpp +++ b/pennylane_lightning/core/src/observables/tests/mpi/Test_ObservablesMPI.cpp @@ -524,3 +524,67 @@ TEST_CASE("Methods implemented in the HamiltonianBase class", testHamiltonianBase(); } } + +template void testSparseHamiltonianBase() { + if constexpr (!std::is_same_v) { + using StateVectorT = typename TypeList::Type; + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + + const std::size_t num_qubits = 3; + std::mt19937 re{1337}; + + MPIManager mpi_manager(MPI_COMM_WORLD); + + size_t mpi_buffersize = 1; + size_t nGlobalIndexBits = + std::bit_width(static_cast(mpi_manager.getSize())) - 1; + size_t nLocalIndexBits = num_qubits - nGlobalIndexBits; + size_t subSvLength = 1 << nLocalIndexBits; + + int nDevices = 0; + cudaGetDeviceCount(&nDevices); + int deviceId = mpi_manager.getRank() % nDevices; + cudaSetDevice(deviceId); + DevTag dt_local(deviceId, 0); + mpi_manager.Barrier(); + + std::vector expected_sv(subSvLength); + std::vector local_state(subSvLength); + + auto init_state = + createRandomStateVectorData(re, num_qubits); + + mpi_manager.Scatter(init_state.data(), local_state.data(), subSvLength, + 0); + mpi_manager.Barrier(); + + DYNAMIC_SECTION("applyInPlace must fail - " + << StateVectorMPIToName::name) { + + auto sparseH = SparseHamiltonianBase::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {7, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, + {0, 1, 2}); + + StateVectorT sv_mpi(mpi_manager, dt_local, mpi_buffersize, + nGlobalIndexBits, nLocalIndexBits); + + sv_mpi.CopyHostDataToGpu(local_state, false); + + REQUIRE_THROWS_AS(sparseH->applyInPlace(sv_mpi), + LightningException); + } + + testSparseHamiltonianBase(); + } +} + +TEST_CASE("Methods implemented in the SparseHamiltonianBase class", + "[SparseHamiltonianBase]") { + if constexpr (BACKEND_FOUND) { + testSparseHamiltonianBase(); + } +} diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindings.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindings.hpp index e713ea2eef..2ebf7d3f95 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindings.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindings.hpp @@ -188,10 +188,13 @@ void registerBackendSpecificMeasurements(PyClass &pyclass) { using np_arr_c = py::array_t, py::array::c_style | py::array::forcecast>; - using sparse_index_type = std::size_t; - using np_arr_sparse_ind = - py::array_t; + using sparse_index_type = + typename std::conditional::value, int32_t, + int64_t>::type; + using np_arr_sparse_ind = typename std::conditional< + std::is_same::value, + py::array_t, + py::array_t>::type; pyclass .def("expval", @@ -205,10 +208,14 @@ void registerBackendSpecificMeasurements(PyClass &pyclass) { const np_arr_sparse_ind &entries, const np_arr_c &values) { return M.expval( static_cast(row_map.request().ptr), - static_cast(row_map.request().size), + static_cast( + row_map.request() + .size), // int64_t is required by cusparse static_cast(entries.request().ptr), static_cast(values.request().ptr), - static_cast(values.request().size)); + static_cast( + values.request() + .size)); // int64_t is required by cusparse }, "Expected value of a sparse Hamiltonian.") .def( @@ -249,14 +256,81 @@ void registerBackendSpecificMeasurements(PyClass &pyclass) { const np_arr_sparse_ind &entries, const np_arr_c &values) { return M.var( static_cast(row_map.request().ptr), - static_cast(row_map.request().size), + static_cast(row_map.request().size), static_cast(entries.request().ptr), static_cast(values.request().ptr), - static_cast(values.request().size)); + static_cast(values.request().size)); }, "Variance of a sparse Hamiltonian."); } +/** + * @brief Register backend specific observables. + * + * @tparam StateVectorT + * @param m Pybind module + */ +template +void registerBackendSpecificObservables(py::module_ &m) { + using PrecisionT = + typename StateVectorT::PrecisionT; // Statevector's precision. + using ComplexT = + typename StateVectorT::ComplexT; // Statevector's complex type. + using ParamT = PrecisionT; // Parameter's data precision + + const std::string bitsize = + std::to_string(sizeof(std::complex) * 8); + + using np_arr_c = py::array_t, py::array::c_style>; + + std::string class_name; + + class_name = "SparseHamiltonianC" + bitsize; + using np_arr_sparse_ind = typename std::conditional< + std::is_same::value, + py::array_t, + py::array_t>::type; + using IdxT = typename SparseHamiltonian::IdxT; + py::class_, + std::shared_ptr>, + Observable>(m, class_name.c_str(), + py::module_local()) + .def(py::init([](const np_arr_c &data, const np_arr_sparse_ind &indices, + const np_arr_sparse_ind &offsets, + const std::vector &wires) { + const py::buffer_info buffer_data = data.request(); + const auto *data_ptr = static_cast(buffer_data.ptr); + + const py::buffer_info buffer_indices = indices.request(); + const auto *indices_ptr = + static_cast(buffer_indices.ptr); + + const py::buffer_info buffer_offsets = offsets.request(); + const auto *offsets_ptr = + static_cast(buffer_offsets.ptr); + + return SparseHamiltonian{ + std::vector({data_ptr, data_ptr + data.size()}), + std::vector({indices_ptr, indices_ptr + indices.size()}), + std::vector({offsets_ptr, offsets_ptr + offsets.size()}), + wires}; + })) + .def("__repr__", &SparseHamiltonian::getObsName) + .def("get_wires", &SparseHamiltonian::getWires, + "Get wires of observables") + .def( + "__eq__", + [](const SparseHamiltonian &self, + py::handle other) -> bool { + if (!py::isinstance>(other)) { + return false; + } + auto other_cast = other.cast>(); + return self == other_cast; + }, + "Compare two observables"); +} + /** * @brief Register backend specific adjoint Jacobian methods. * diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindingsMPI.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindingsMPI.hpp index c79ad88fc5..1ca4670fe7 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindingsMPI.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/bindings/LGPUBindingsMPI.hpp @@ -191,10 +191,13 @@ void registerBackendSpecificMeasurementsMPI(PyClass &pyclass) { using np_arr_c = py::array_t, py::array::c_style | py::array::forcecast>; - using sparse_index_type = std::size_t; - using np_arr_sparse_ind = - py::array_t; + using sparse_index_type = + typename std::conditional::value, int32_t, + int64_t>::type; + using np_arr_sparse_ind = typename std::conditional< + std::is_same::value, + py::array_t, + py::array_t>::type; pyclass .def("expval", @@ -209,10 +212,10 @@ void registerBackendSpecificMeasurementsMPI(PyClass &pyclass) { const np_arr_sparse_ind &entries, const np_arr_c &values) { return M.expval( static_cast(row_map.request().ptr), - static_cast(row_map.request().size), + static_cast(row_map.request().size), static_cast(entries.request().ptr), static_cast(values.request().ptr), - static_cast(values.request().size)); + static_cast(values.request().size)); }, "Expected value of a sparse Hamiltonian.") .def( @@ -254,10 +257,10 @@ void registerBackendSpecificMeasurementsMPI(PyClass &pyclass) { const np_arr_sparse_ind &entries, const np_arr_c &values) { return M.var( static_cast(row_map.request().ptr), - static_cast(row_map.request().size), + static_cast(row_map.request().size), static_cast(entries.request().ptr), static_cast(values.request().ptr), - static_cast(values.request().size)); + static_cast(values.request().size)); }, "Variance of a sparse Hamiltonian."); } 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 af2ae38ec2..8ca6eacc69 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPU.hpp @@ -259,10 +259,10 @@ class Measurements final * @return auto Expectation value. */ template - auto expval(const index_type *csrOffsets_ptr, - const index_type csrOffsets_size, const index_type *columns_ptr, + auto expval(const index_type *csrOffsets_ptr, const int64_t csrOffsets_size, + const index_type *columns_ptr, const std::complex *values_ptr, - const index_type numNNZ) -> PrecisionT { + const int64_t numNNZ) -> PrecisionT { const std::size_t nIndexBits = this->_statevector.getNumQubits(); const std::size_t length = std::size_t{1} << nIndexBits; @@ -580,10 +580,10 @@ class Measurements final * @return Floating point with the variance of the sparse Hamiltonian. */ template - PrecisionT - var(const index_type *csrOffsets_ptr, const index_type csrOffsets_size, - const index_type *columns_ptr, - const std::complex *values_ptr, const index_type numNNZ) { + PrecisionT var(const index_type *csrOffsets_ptr, + const int64_t csrOffsets_size, const index_type *columns_ptr, + const std::complex *values_ptr, + const int64_t numNNZ) { PL_ABORT_IF( (this->_statevector.getLength() != (size_t(csrOffsets_size) - 1)), "Statevector and Hamiltonian have incompatible sizes."); diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp index f96e2bc217..ff101654df 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/MeasurementsGPUMPI.hpp @@ -47,6 +47,7 @@ using namespace Pennylane; using namespace Pennylane::Measures; using namespace Pennylane::Observables; using namespace Pennylane::LightningGPU::Observables; +using namespace Pennylane::LightningGPU::MPI; namespace cuUtil = Pennylane::LightningGPU::Util; using Pennylane::LightningGPU::StateVectorCudaManaged; using namespace Pennylane::Util; @@ -366,10 +367,10 @@ class MeasurementsMPI final * @return auto Expectation value. */ template - auto expval(const index_type *csrOffsets_ptr, - const index_type csrOffsets_size, const index_type *columns_ptr, + auto expval(const index_type *csrOffsets_ptr, const int64_t csrOffsets_size, + const index_type *columns_ptr, const std::complex *values_ptr, - const index_type numNNZ) -> PrecisionT { + const int64_t numNNZ) -> PrecisionT { if (mpi_manager_.getRank() == 0) { PL_ABORT_IF_NOT( static_cast(csrOffsets_size - 1) == @@ -657,10 +658,10 @@ class MeasurementsMPI final * @return Floating point with the variance of the sparse Hamiltonian. */ template - PrecisionT - var(const index_type *csrOffsets_ptr, const index_type csrOffsets_size, - const index_type *columns_ptr, - const std::complex *values_ptr, const index_type numNNZ) { + PrecisionT var(const index_type *csrOffsets_ptr, + const int64_t csrOffsets_size, const index_type *columns_ptr, + const std::complex *values_ptr, + const int64_t numNNZ) { if (mpi_manager_.getRank() == 0) { PL_ABORT_IF_NOT( static_cast(csrOffsets_size - 1) == diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/Test_StateVectorCudaManaged_Expval.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/Test_StateVectorCudaManaged_Expval.cpp index 82b4d41e93..36f1f1f128 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/Test_StateVectorCudaManaged_Expval.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/Test_StateVectorCudaManaged_Expval.cpp @@ -336,6 +336,8 @@ TEMPLATE_TEST_CASE("StateVectorCudaManaged::Hamiltonian_expval_Sparse", "[StateVectorCudaManaged_Expval]", float, double) { using StateVectorT = StateVectorCudaManaged; using ComplexT = StateVectorT::ComplexT; + using IdxT = typename std::conditional::value, + int32_t, int64_t>::type; SECTION("Sparse expval") { std::vector init_state{{0.0, 0.0}, {0.0, 0.1}, {0.1, 0.1}, @@ -344,17 +346,18 @@ TEMPLATE_TEST_CASE("StateVectorCudaManaged::Hamiltonian_expval_Sparse", StateVectorT sv{init_state.data(), init_state.size()}; auto m = Measurements(sv); - std::vector index_ptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; - std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, - 4, 7, 5, 6, 5, 6, 4, 7}; + std::vector index_ptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; + std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, + 4, 7, 5, 6, 5, 6, 4, 7}; std::vector values = { {3.1415, 0.0}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {3.1415, 0.0}, {3.1415, 0.0}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {3.1415, 0.0}}; - auto result = m.expval(index_ptr.data(), index_ptr.size(), - indices.data(), values.data(), values.size()); + auto result = m.expval( + index_ptr.data(), static_cast(index_ptr.size()), + indices.data(), values.data(), static_cast(values.size())); auto expected = TestType(3.1415); CHECK(expected == Approx(result).epsilon(1e-7)); } @@ -372,22 +375,23 @@ TEMPLATE_TEST_CASE("StateVectorCudaManaged::Hamiltonian_expval_Sparse", // measurements. Measurements Measurer(sv); const size_t num_qubits = 3; - const size_t data_size = Pennylane::Util::exp2(num_qubits); + size_t data_size = Pennylane::Util::exp2(num_qubits); - std::vector row_map; - std::vector entries; + std::vector row_map; + std::vector entries; std::vector values; - write_CSR_vectors(row_map, entries, values, data_size); + write_CSR_vectors(row_map, entries, values, + static_cast(data_size)); - PrecisionT exp_values = - Measurer.expval(row_map.data(), row_map.size(), entries.data(), - values.data(), values.size()); + PrecisionT exp_values = Measurer.expval( + row_map.data(), static_cast(row_map.size()), + entries.data(), values.data(), static_cast(values.size())); PrecisionT exp_values_ref = 0.5930885; REQUIRE(exp_values == Approx(exp_values_ref).margin(1e-6)); - PrecisionT var_values = - Measurer.var(row_map.data(), row_map.size(), entries.data(), - values.data(), values.size()); + PrecisionT var_values = Measurer.var( + row_map.data(), static_cast(row_map.size()), + entries.data(), values.data(), static_cast(values.size())); PrecisionT var_values_ref = 2.4624654; REQUIRE(var_values == Approx(var_values_ref).margin(1e-6)); } diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/mpi/Test_StateVectorCudaMPI_Expval.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/mpi/Test_StateVectorCudaMPI_Expval.cpp index 47505e73fe..b6fdab8737 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/mpi/Test_StateVectorCudaMPI_Expval.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/measurements/tests/mpi/Test_StateVectorCudaMPI_Expval.cpp @@ -401,6 +401,8 @@ TEMPLATE_TEST_CASE("StateVectorCudaMPI::Hamiltonian_expval_Sparse", "[StateVectorCudaMPI_Expval]", double) { using StateVectorT = StateVectorCudaMPI; using ComplexT = StateVectorT::ComplexT; + using IdxT = typename std::conditional::value, + int32_t, int64_t>::type; MPIManager mpi_manager(MPI_COMM_WORLD); REQUIRE(mpi_manager.getSize() == 2); @@ -431,17 +433,18 @@ TEMPLATE_TEST_CASE("StateVectorCudaMPI::Hamiltonian_expval_Sparse", sv.CopyHostDataToGpu(local_init_sv.data(), local_init_sv.size(), false); auto m = MeasurementsMPI(sv); - std::vector index_ptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; - std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, - 4, 7, 5, 6, 5, 6, 4, 7}; + std::vector index_ptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; + std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, + 4, 7, 5, 6, 5, 6, 4, 7}; std::vector values = { {3.1415, 0.0}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {3.1415, 0.0}, {3.1415, 0.0}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {0.0, -3.1415}, {3.1415, 0.0}, {0.0, 3.1415}, {3.1415, 0.0}}; - auto result = m.expval(index_ptr.data(), index_ptr.size(), - indices.data(), values.data(), values.size()); + auto result = m.expval( + index_ptr.data(), static_cast(index_ptr.size()), + indices.data(), values.data(), static_cast(values.size())); auto expected = TestType(3.1415); CHECK(expected == Approx(result).epsilon(1e-7)); } @@ -464,24 +467,25 @@ TEMPLATE_TEST_CASE("StateVectorCudaMPI::Hamiltonian_expval_Sparse", // This object attaches to the statevector allowing several // measurements. MeasurementsMPI Measurer(sv); - const size_t data_size = Pennylane::Util::exp2(num_qubits); + size_t data_size = Pennylane::Util::exp2(num_qubits); - std::vector row_map; - std::vector entries; + std::vector row_map; + std::vector entries; std::vector values; - write_CSR_vectors(row_map, entries, values, data_size); + write_CSR_vectors(row_map, entries, values, + static_cast(data_size)); - PrecisionT exp_values = - Measurer.expval(row_map.data(), row_map.size(), entries.data(), - values.data(), values.size()); + PrecisionT exp_values = Measurer.expval( + row_map.data(), static_cast(row_map.size()), + entries.data(), values.data(), static_cast(values.size())); PrecisionT exp_values_ref = 0.5930885; REQUIRE(exp_values == Approx(exp_values_ref).margin(1e-6)); mpi_manager.Barrier(); - PrecisionT var_values = - Measurer.var(row_map.data(), row_map.size(), entries.data(), - values.data(), values.size()); + PrecisionT var_values = Measurer.var( + row_map.data(), static_cast(row_map.size()), + entries.data(), values.data(), static_cast(values.size())); PrecisionT var_values_ref = 2.4624654; REQUIRE(var_values == Approx(var_values_ref).margin(1e-6)); } diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.cpp index b186caa39e..c4f2ca82d4 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.cpp @@ -28,3 +28,6 @@ template class Observables::TensorProdObs>; template class Observables::Hamiltonian>; template class Observables::Hamiltonian>; + +template class Observables::SparseHamiltonian>; +template class Observables::SparseHamiltonian>; diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp index c92e7dac5d..6d0d0a94e7 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPU.hpp @@ -209,4 +209,87 @@ class Hamiltonian final : public HamiltonianBase { } }; +/** + * @brief Sparse representation of Hamiltonian + * + */ +template +class SparseHamiltonian final : public SparseHamiltonianBase { + private: + using BaseType = SparseHamiltonianBase; + + public: + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + // cuSparse required index type + using IdxT = typename BaseType::IdxT; + + /** + * @brief Create a SparseHamiltonian from data, indices and offsets in CSR + * format. + * + * @param data Arguments to construct data + * @param indices Arguments to construct indices + * @param offsets Arguments to construct offsets + * @param wires Arguments to construct wires + */ + template + explicit SparseHamiltonian(T1 &&data, T2 &&indices, T3 &&offsets, + T4 &&wires) + : BaseType{data, indices, offsets, wires} {} + + /** + * @brief Convenient wrapper for the constructor as the constructor does not + * convert the std::shared_ptr with a derived class correctly. + * + * This function is useful as std::make_shared does not handle + * brace-enclosed initializer list correctly. + * + * @param data Argument to construct data + * @param indices Argument to construct indices + * @param offsets Argument to construct ofsets + * @param wires Argument to construct wires + */ + static auto create(std::initializer_list data, + std::initializer_list indices, + std::initializer_list offsets, + std::initializer_list wires) + -> std::shared_ptr> { + return std::shared_ptr>( + new SparseHamiltonian{ + std::move(data), std::move(indices), std::move(offsets), + std::move(wires)}); + } + + /** + * @brief Updates the statevector SV:->SV', where SV' = a*H*SV, and where H + * is a sparse Hamiltonian. + * + */ + void applyInPlace(StateVectorT &sv) const override { + PL_ABORT_IF_NOT(this->wires_.size() == sv.getNumQubits(), + "SparseH wire count does not match state-vector size"); + using CFP_t = typename StateVectorT::CFP_t; + + const std::size_t nIndexBits = sv.getNumQubits(); + const std::size_t length = std::size_t{1} << nIndexBits; + + auto device_id = sv.getDataBuffer().getDevTag().getDeviceID(); + auto stream_id = sv.getDataBuffer().getDevTag().getStreamID(); + + cusparseHandle_t handle = sv.getCusparseHandle(); + + std::unique_ptr> d_sv_prime = + std::make_unique>(length, device_id, stream_id, + true); + + SparseMV_cuSparse( + this->offsets_.data(), static_cast(this->offsets_.size()), + this->indices_.data(), this->data_.data(), + static_cast(this->data_.size()), sv.getData(), + d_sv_prime->getData(), device_id, stream_id, handle); + sv.updateData(std::move(d_sv_prime)); + } +}; + } // namespace Pennylane::LightningGPU::Observables diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.cpp index 9b1f776e0f..ae9ac9100a 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.cpp @@ -28,3 +28,6 @@ template class Observables::TensorProdObsMPI>; template class Observables::HamiltonianMPI>; template class Observables::HamiltonianMPI>; + +template class Observables::SparseHamiltonianMPI>; +template class Observables::SparseHamiltonianMPI>; diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp index d15df18207..94f5e45739 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/ObservablesGPUMPI.hpp @@ -21,6 +21,7 @@ #include "Constant.hpp" #include "ConstantUtil.hpp" // lookup #include "LinearAlg.hpp" +#include "MPILinearAlg.hpp" #include "Observables.hpp" #include "StateVectorCudaMPI.hpp" #include "Util.hpp" @@ -213,4 +214,97 @@ class HamiltonianMPI final : public HamiltonianBase { } }; +/** + * @brief Sparse representation of Hamiltonian + * + */ +template +class SparseHamiltonianMPI final : public SparseHamiltonianBase { + public: + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + // cuSparse required index type + using IdxT = + typename std::conditional::value, + int32_t, int64_t>::type; + + private: + using BaseType = SparseHamiltonianBase; + + public: + /** + * @brief Create a SparseHamiltonianMPI from data, indices and offsets in + * CSR format. + * + * @param data Arguments to construct data + * @param indices Arguments to construct indices + * @param offsets Arguments to construct offsets + * @param wires Arguments to construct wires + */ + template + explicit SparseHamiltonianMPI(T1 &&data, T2 &&indices, T3 &&offsets, + T4 &&wires) + : BaseType{data, indices, offsets, wires} {} + + /** + * @brief Convenient wrapper for the constructor as the constructor does not + * convert the std::shared_ptr with a derived class correctly. + * + * This function is useful as std::make_shared does not handle + * brace-enclosed initializer list correctly. + * + * @param data Argument to construct data + * @param indices Argument to construct indices + * @param offsets Argument to construct ofsets + * @param wires Argument to construct wires + */ + static auto create(std::initializer_list data, + std::initializer_list indices, + std::initializer_list offsets, + std::initializer_list wires) + -> std::shared_ptr> { + return std::shared_ptr>( + new SparseHamiltonianMPI{ + std::move(data), std::move(indices), std::move(offsets), + std::move(wires)}); + } + + /** + * @brief Updates the statevector SV:->SV', where SV' = a*H*SV, and where H + * is a sparse Hamiltonian. + * + */ + void applyInPlace(StateVectorT &sv) const override { + auto mpi_manager = sv.getMPIManager(); + if (mpi_manager.getRank() == 0) { + PL_ABORT_IF_NOT( + this->wires_.size() == sv.getTotalNumQubits(), + "SparseH wire count does not match state-vector size"); + } + using CFP_t = typename StateVectorT::CFP_t; + + auto device_id = sv.getDataBuffer().getDevTag().getDeviceID(); + auto stream_id = sv.getDataBuffer().getDevTag().getStreamID(); + + const size_t length_local = size_t{1} << sv.getNumLocalQubits(); + + std::unique_ptr> d_sv_prime = + std::make_unique>(length_local, device_id, + stream_id, true); + d_sv_prime->zeroInit(); + PL_CUDA_IS_SUCCESS(cudaDeviceSynchronize()); + mpi_manager.Barrier(); + + cuUtil::SparseMV_cuSparseMPI( + mpi_manager, length_local, this->offsets_.data(), + static_cast(this->offsets_.size()), this->indices_.data(), + this->data_.data(), const_cast(sv.getData()), + d_sv_prime->getData(), device_id, stream_id, + sv.getCusparseHandle()); + + sv.CopyGpuDataToGpuIn(d_sv_prime->getData(), d_sv_prime->getLength()); + mpi_manager.Barrier(); + } +}; + } // namespace Pennylane::LightningGPU::Observables diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/Test_ObservablesGPU.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/Test_ObservablesGPU.cpp index 0d8bd7d388..398f664ffc 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/Test_ObservablesGPU.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/Test_ObservablesGPU.cpp @@ -154,6 +154,20 @@ TEMPLATE_PRODUCT_TEST_CASE("Hamiltonian", "[Observables]", } } +TEMPLATE_PRODUCT_TEST_CASE("SparseHamiltonian", "[Observables]", + (StateVectorCudaManaged), (float, double)) { + using StateVectorT = TestType; + using SparseHamiltonianT = SparseHamiltonian; + + SECTION("Copy constructibility") { + REQUIRE(std::is_copy_constructible_v); + } + + SECTION("Move constructibility") { + REQUIRE(std::is_move_constructible_v); + } +} + TEMPLATE_PRODUCT_TEST_CASE("Observables::HermitianHasher", "[Observables]", (StateVectorCudaManaged), (float, double)) { using StateVectorT = TestType; @@ -257,3 +271,31 @@ TEMPLATE_PRODUCT_TEST_CASE("Hamiltonian::ApplyInPlace", "[Observables]", } } } + +TEMPLATE_PRODUCT_TEST_CASE("SparseHamiltonian::ApplyInPlace", "[Observables]", + (StateVectorCudaManaged), (float, double)) { + using StateVectorT = TestType; + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + + const std::size_t num_qubits = 3; + std::mt19937 re{1337}; + + auto sparseH = SparseHamiltonian::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {7, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, {0, 1, 2}); + + auto init_state = createRandomStateVectorData(re, num_qubits); + + StateVectorT state_vector(init_state.data(), init_state.size()); + + sparseH->applyInPlace(state_vector); + + std::reverse(init_state.begin(), init_state.end()); + + REQUIRE(isApproxEqual(state_vector.getDataVector().data(), + state_vector.getDataVector().size(), + init_state.data(), init_state.size())); +} \ No newline at end of file diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/mpi/Test_ObservablesGPUMPI.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/mpi/Test_ObservablesGPUMPI.cpp index 6abde861e6..bc4d0e517f 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/mpi/Test_ObservablesGPUMPI.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/observables/tests/mpi/Test_ObservablesGPUMPI.cpp @@ -289,4 +289,58 @@ TEMPLATE_PRODUCT_TEST_CASE("Observables::HermitianHasherMPI", "[Observables]", CHECK(ham_1->getObsName() == res1.str()); CHECK(ham_2->getObsName() == res2.str()); } +} + +TEMPLATE_PRODUCT_TEST_CASE("SparseHamiltonian::ApplyInPlace", "[Observables]", + (StateVectorCudaMPI), (float, double)) { + using StateVectorT = TestType; + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + MPIManager mpi_manager(MPI_COMM_WORLD); + + const std::size_t num_qubits = 3; + std::mt19937 re{1337}; + + auto sparseH = SparseHamiltonianMPI::create( + {ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}, + ComplexT{1.0, 0.0}, ComplexT{1.0, 0.0}}, + {7, 6, 5, 4, 3, 2, 1, 0}, {0, 1, 2, 3, 4, 5, 6, 7, 8}, {0, 1, 2}); + + size_t mpi_buffersize = 1; + size_t nGlobalIndexBits = + std::bit_width(static_cast(mpi_manager.getSize())) - 1; + size_t nLocalIndexBits = num_qubits - nGlobalIndexBits; + size_t subSvLength = 1 << nLocalIndexBits; + + mpi_manager.Barrier(); + std::vector expected_sv(subSvLength); + std::vector local_state(subSvLength); + + auto init_state = createRandomStateVectorData(re, num_qubits); + + mpi_manager.Scatter(init_state.data(), local_state.data(), subSvLength, 0); + mpi_manager.Barrier(); + + int nDevices = 0; + cudaGetDeviceCount(&nDevices); + int deviceId = mpi_manager.getRank() % nDevices; + cudaSetDevice(deviceId); + DevTag dt_local(deviceId, 0); + mpi_manager.Barrier(); + + StateVectorT sv_mpi(mpi_manager, dt_local, mpi_buffersize, nGlobalIndexBits, + nLocalIndexBits); + + sv_mpi.CopyHostDataToGpu(local_state, false); + + sparseH->applyInPlace(sv_mpi); + + std::reverse(init_state.begin(), init_state.end()); + mpi_manager.Scatter(init_state.data(), expected_sv.data(), subSvLength, 0); + mpi_manager.Barrier(); + + REQUIRE(isApproxEqual(sv_mpi.getDataVector().data(), + sv_mpi.getDataVector().size(), expected_sv.data(), + expected_sv.size())); } \ No newline at end of file diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/LinearAlg.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/LinearAlg.hpp index f70d4ea9f2..1827bfffa8 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/LinearAlg.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/LinearAlg.hpp @@ -307,20 +307,20 @@ inline SharedCusparseHandle make_shared_cusparse_handle() { * @param handle cuSparse handle. */ template -inline void SparseMV_cuSparse(const index_type *csrOffsets_ptr, - const index_type csrOffsets_size, - const index_type *columns_ptr, - const std::complex *values_ptr, - const index_type numNNZ, CFP_t *X, CFP_t *Y, - DevTypeID device_id, cudaStream_t stream_id, - cusparseHandle_t handle) { - const int64_t num_rows = static_cast( +inline void +SparseMV_cuSparse(const index_type *csrOffsets_ptr, + const int64_t csrOffsets_size, const index_type *columns_ptr, + const std::complex *values_ptr, + const int64_t numNNZ, CFP_t *X, CFP_t *Y, DevTypeID device_id, + cudaStream_t stream_id, cusparseHandle_t handle) { + + const int64_t num_rows = csrOffsets_size - - 1); // int64_t is required for num_rows by cusparseCreateCsr - const int64_t num_cols = static_cast( - num_rows); // int64_t is required for num_cols by cusparseCreateCsr - const int64_t nnz = static_cast( - numNNZ); // int64_t is required for nnz by cusparseCreateCsr + 1; // int64_t is required for num_rows by cusparseCreateCsr + const int64_t num_cols = + num_rows; // int64_t is required for num_cols by cusparseCreateCsr + const int64_t nnz = + numNNZ; // int64_t is required for nnz by cusparseCreateCsr const CFP_t alpha = {1.0, 0.0}; const CFP_t beta = {0.0, 0.0}; @@ -338,13 +338,15 @@ inline void SparseMV_cuSparse(const index_type *csrOffsets_ptr, d_values.CopyHostDataToGpu(values_ptr, d_values.getLength(), false); cudaDataType_t data_type; - cusparseIndexType_t compute_type = CUSPARSE_INDEX_64I; + cusparseIndexType_t compute_type; if constexpr (std::is_same_v || std::is_same_v) { data_type = CUDA_C_64F; + compute_type = CUSPARSE_INDEX_64I; } else { data_type = CUDA_C_32F; + compute_type = CUSPARSE_INDEX_32I; } // CUSPARSE APIs @@ -394,8 +396,7 @@ inline void SparseMV_cuSparse(const index_type *csrOffsets_ptr, /* cusparseSpMVAlg_t */ CUSPARSE_SPMV_ALG_DEFAULT, /* size_t* */ &bufferSize)); - DataBuffer dBuffer{bufferSize, device_id, stream_id, - true}; + DataBuffer dBuffer{bufferSize, device_id, stream_id, true}; // execute SpMV PL_CUSPARSE_IS_SUCCESS(cusparseSpMV( @@ -439,19 +440,19 @@ inline void SparseMV_cuSparse(const index_type *csrOffsets_ptr, */ template inline void SparseMV_cuSparse(const index_type *csrOffsets_ptr, - const index_type csrOffsets_size, + const int64_t csrOffsets_size, const index_type *columns_ptr, const std::complex *values_ptr, - const index_type numNNZ, const CFP_t *X, CFP_t *Y, + const int64_t numNNZ, const CFP_t *X, CFP_t *Y, DevTypeID device_id, cudaStream_t stream_id, cusparseHandle_t handle) { - const int64_t num_rows = static_cast( + const int64_t num_rows = csrOffsets_size - - 1); // int64_t is required for num_rows by cusparseCreateCsr - const int64_t num_cols = static_cast( - num_rows); // int64_t is required for num_cols by cusparseCreateCsr - const int64_t nnz = static_cast( - numNNZ); // int64_t is required for nnz by cusparseCreateCsr + 1; // int64_t is required for num_rows by cusparseCreateCsr + const int64_t num_cols = + num_rows; // int64_t is required for num_cols by cusparseCreateCsr + const int64_t nnz = + numNNZ; // int64_t is required for nnz by cusparseCreateCsr const CFP_t alpha = {1.0, 0.0}; const CFP_t beta = {0.0, 0.0}; diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/MPILinearAlg.hpp b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/MPILinearAlg.hpp index 4b2e905b8e..cd2afd426b 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/MPILinearAlg.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/MPILinearAlg.hpp @@ -45,15 +45,15 @@ namespace Pennylane::LightningGPU::Util { template inline void SparseMV_cuSparseMPI( MPIManager &mpi_manager, const size_t &length_local, - const index_type *csrOffsets_ptr, const index_type csrOffsets_size, + const index_type *csrOffsets_ptr, const int64_t csrOffsets_size, const index_type *columns_ptr, const std::complex *values_ptr, CFP_t *X, CFP_t *Y, DevTypeID device_id, cudaStream_t stream_id, cusparseHandle_t handle) { std::vector>> csrmatrix_blocks; if (mpi_manager.getRank() == 0) { csrmatrix_blocks = splitCSRMatrix( - mpi_manager, csrOffsets_size - 1, csrOffsets_ptr, columns_ptr, - values_ptr); + mpi_manager, static_cast(csrOffsets_size - 1), + csrOffsets_ptr, columns_ptr, values_ptr); } mpi_manager.Barrier(); @@ -79,11 +79,11 @@ inline void SparseMV_cuSparseMPI( color = 1; SparseMV_cuSparse( localCSRMatrix.getCsrOffsets().data(), - localCSRMatrix.getCsrOffsets().size(), + static_cast(localCSRMatrix.getCsrOffsets().size()), localCSRMatrix.getColumns().data(), localCSRMatrix.getValues().data(), - localCSRMatrix.getValues().size(), X, d_res_per_block.getData(), - device_id, stream_id, handle); + static_cast(localCSRMatrix.getValues().size()), X, + d_res_per_block.getData(), device_id, stream_id, handle); } PL_CUDA_IS_SUCCESS(cudaDeviceSynchronize()); diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/Test_LinearAlgebra.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/Test_LinearAlgebra.cpp index cd1dd3937a..a2b35d0742 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/Test_LinearAlgebra.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/Test_LinearAlgebra.cpp @@ -40,6 +40,8 @@ TEMPLATE_TEST_CASE("Linear Algebra::SparseMV", "[Linear Algebra]", float, using StateVectorT = StateVectorCudaManaged; using ComplexT = StateVectorT::ComplexT; using CFP_t = StateVectorT::CFP_t; + using IdxT = typename std::conditional::value, + int32_t, int64_t>::type; std::size_t num_qubits = 3; std::size_t data_size = exp2(num_qubits); @@ -52,9 +54,9 @@ TEMPLATE_TEST_CASE("Linear Algebra::SparseMV", "[Linear Algebra]", float, {0.2, -0.1}, {-0.1, 0.2}, {0.2, 0.1}, {0.1, 0.2}, {0.7, -0.2}, {-0.1, 0.6}, {0.6, 0.1}, {0.2, 0.7}}; - std::vector indptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; - std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, - 4, 7, 5, 6, 5, 6, 4, 7}; + std::vector indptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; + std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, + 4, 7, 5, 6, 5, 6, 4, 7}; std::vector values = { {1.0, 0.0}, {0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}, {0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}, {1.0, 0.0}, @@ -69,10 +71,10 @@ TEMPLATE_TEST_CASE("Linear Algebra::SparseMV", "[Linear Algebra]", float, SECTION("Testing sparse matrix vector product:") { std::vector result(data_size); - cuUtil::SparseMV_cuSparse( - indptr.data(), indptr.size(), indices.data(), values.data(), - values.size(), sv_x.getData(), sv_y.getData(), - sv_x.getDataBuffer().getDevTag().getDeviceID(), + cuUtil::SparseMV_cuSparse( + indptr.data(), static_cast(indptr.size()), indices.data(), + values.data(), static_cast(values.size()), sv_x.getData(), + sv_y.getData(), sv_x.getDataBuffer().getDevTag().getDeviceID(), sv_x.getDataBuffer().getDevTag().getStreamID(), sv_x.getCusparseHandle()); diff --git a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/mpi/Test_LinearAlgebraMPI.cpp b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/mpi/Test_LinearAlgebraMPI.cpp index 64df9d77e8..371bc66f75 100644 --- a/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/mpi/Test_LinearAlgebraMPI.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_gpu/utils/tests/mpi/Test_LinearAlgebraMPI.cpp @@ -40,6 +40,8 @@ TEMPLATE_TEST_CASE("Linear Algebra::SparseMV", "[Linear Algebra]", float, using StateVectorT = StateVectorCudaMPI; using ComplexT = StateVectorT::ComplexT; using CFP_t = StateVectorT::CFP_t; + using IdxT = typename std::conditional::value, + int32_t, int64_t>::type; MPIManager mpi_manager(MPI_COMM_WORLD); REQUIRE(mpi_manager.getSize() == 2); @@ -54,9 +56,9 @@ TEMPLATE_TEST_CASE("Linear Algebra::SparseMV", "[Linear Algebra]", float, {0.1, 0.2}, {0.7, -0.2}, {-0.1, 0.6}, {0.6, 0.1}, {0.2, 0.7}}; - std::vector indptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; - std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, - 4, 7, 5, 6, 5, 6, 4, 7}; + std::vector indptr = {0, 2, 4, 6, 8, 10, 12, 14, 16}; + std::vector indices = {0, 3, 1, 2, 1, 2, 0, 3, + 4, 7, 5, 6, 5, 6, 4, 7}; std::vector values = { {1.0, 0.0}, {0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}, {0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}, {1.0, 0.0}, @@ -95,9 +97,10 @@ TEMPLATE_TEST_CASE("Linear Algebra::SparseMV", "[Linear Algebra]", float, nGlobalIndexBits, nLocalIndexBits); sv_x.CopyHostDataToGpu(local_state, false); - cuUtil::SparseMV_cuSparseMPI( - mpi_manager, sv_x.getLength(), indptr.data(), indptr.size(), - indices.data(), values.data(), sv_x.getData(), sv_y.getData(), + cuUtil::SparseMV_cuSparseMPI( + mpi_manager, sv_x.getLength(), indptr.data(), + static_cast(indptr.size()), indices.data(), values.data(), + sv_x.getData(), sv_y.getData(), sv_x.getDataBuffer().getDevTag().getDeviceID(), sv_x.getDataBuffer().getDevTag().getStreamID(), sv_x.getCusparseHandle()); diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp index a3f0951b24..02064e811f 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/StateVectorKokkos.hpp @@ -17,6 +17,7 @@ */ #pragma once +#include #include #include #include @@ -170,7 +171,7 @@ class StateVectorKokkos final * * @param num_qubits Number of qubits */ - StateVectorKokkos(ComplexT *hostdata_, size_t length, + StateVectorKokkos(ComplexT *hostdata_, std::size_t length, const Kokkos::InitializationSettings &kokkos_args = {}) : StateVectorKokkos(log2(length), kokkos_args) { PL_ABORT_IF_NOT(isPerfectPowerOf2(length), @@ -178,12 +179,20 @@ class StateVectorKokkos final HostToDevice(hostdata_, length); } + StateVectorKokkos(std::complex *hostdata_, std::size_t length, + const Kokkos::InitializationSettings &kokkos_args = {}) + : StateVectorKokkos(log2(length), kokkos_args) { + PL_ABORT_IF_NOT(isPerfectPowerOf2(length), + "The size of provided data must be a power of 2."); + HostToDevice(reinterpret_cast(hostdata_), length); + } + /** * @brief Create a new state vector from data on the host. * * @param num_qubits Number of qubits */ - StateVectorKokkos(const ComplexT *hostdata_, size_t length, + StateVectorKokkos(const ComplexT *hostdata_, std::size_t length, const Kokkos::InitializationSettings &kokkos_args = {}) : StateVectorKokkos(log2(length), kokkos_args) { PL_ABORT_IF_NOT(isPerfectPowerOf2(length), @@ -692,7 +701,7 @@ class StateVectorKokkos final * @param new_data data pointer to new data. * @param new_size size of underlying data storage. */ - void updateData(ComplexT *new_data, size_t new_size) { + void updateData(ComplexT *new_data, std::size_t new_size) { updateData(KokkosVector(new_data, new_size)); } @@ -744,7 +753,7 @@ class StateVectorKokkos final * @brief Copy data from the host space to the device space. * */ - inline void HostToDevice(ComplexT *sv, size_t length) { + inline void HostToDevice(ComplexT *sv, std::size_t length) { Kokkos::deep_copy(*data_, UnmanagedComplexHostView(sv, length)); } @@ -752,7 +761,7 @@ class StateVectorKokkos final * @brief Copy data from the device space to the host space. * */ - inline void DeviceToHost(ComplexT *sv, size_t length) const { + inline void DeviceToHost(ComplexT *sv, std::size_t length) const { Kokkos::deep_copy(UnmanagedComplexHostView(sv, length), *data_); } diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/bindings/LKokkosBindings.hpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/bindings/LKokkosBindings.hpp index bd9d89d72f..6432864c4e 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/bindings/LKokkosBindings.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/bindings/LKokkosBindings.hpp @@ -24,6 +24,7 @@ #include "ConstantUtil.hpp" // lookup #include "GateOperation.hpp" #include "MeasurementsKokkos.hpp" +#include "ObservablesKokkos.hpp" #include "StateVectorKokkos.hpp" #include "TypeList.hpp" #include "Util.hpp" // exp2 @@ -33,6 +34,7 @@ namespace { using namespace Pennylane::Bindings; using namespace Pennylane::LightningKokkos::Algorithms; using namespace Pennylane::LightningKokkos::Measures; +using namespace Pennylane::LightningKokkos::Observables; using Kokkos::InitializationSettings; using Pennylane::LightningKokkos::StateVectorKokkos; using Pennylane::Util::exp2; @@ -214,6 +216,58 @@ void registerBackendSpecificMeasurements(PyClass &pyclass) { "Variance of a sparse Hamiltonian."); } +/** + * @brief Register observable classes. + * + * @tparam StateVectorT + * @param m Pybind module + */ +template +void registerBackendSpecificObservables(py::module_ &m) { + using PrecisionT = + typename StateVectorT::PrecisionT; // Statevector's precision. + using ParamT = PrecisionT; // Parameter's data precision + + const std::string bitsize = + std::to_string(sizeof(std::complex) * 8); + + using np_arr_c = py::array_t, py::array::c_style>; + + std::string class_name; + + class_name = "SparseHamiltonianC" + bitsize; + py::class_, + std::shared_ptr>, + Observable>(m, class_name.c_str(), + py::module_local()) + .def(py::init([](const np_arr_c &data, + const std::vector &indices, + const std::vector &indptr, + const std::vector &wires) { + using ComplexT = typename StateVectorT::ComplexT; + const py::buffer_info buffer_data = data.request(); + const auto *data_ptr = static_cast(buffer_data.ptr); + + return SparseHamiltonian{ + std::vector({data_ptr, data_ptr + data.size()}), + indices, indptr, wires}; + })) + .def("__repr__", &SparseHamiltonian::getObsName) + .def("get_wires", &SparseHamiltonian::getWires, + "Get wires of observables") + .def( + "__eq__", + [](const SparseHamiltonian &self, + py::handle other) -> bool { + if (!py::isinstance>(other)) { + return false; + } + auto other_cast = other.cast>(); + return self == other_cast; + }, + "Compare two observables"); +} + /** * @brief Register backend specific adjoint Jacobian methods. * diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.cpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.cpp index d90f3e6019..66192b934a 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.cpp @@ -28,3 +28,6 @@ template class Observables::TensorProdObs>; template class Observables::Hamiltonian>; template class Observables::Hamiltonian>; + +template class Observables::SparseHamiltonian>; +template class Observables::SparseHamiltonian>; diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp index a0371df7f2..c3fae6b3ea 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/ObservablesKokkos.hpp @@ -31,6 +31,7 @@ namespace { using namespace Pennylane::Util; using namespace Pennylane::Observables; using Pennylane::LightningKokkos::StateVectorKokkos; +using Pennylane::LightningKokkos::Util::SparseMV_Kokkos; } // namespace /// @endcond @@ -199,6 +200,76 @@ class Hamiltonian final : public HamiltonianBase { } }; +/** + * @brief Sparse representation of Hamiltonian + * + */ +template +class SparseHamiltonian final : public SparseHamiltonianBase { + private: + using BaseType = SparseHamiltonianBase; + + public: + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + using IdxT = typename BaseType::IdxT; + + /** + * @brief Create a SparseHamiltonian from data, indices and offsets in CSR + * format. + * + * @param data Arguments to construct data + * @param indices Arguments to construct indices + * @param offsets Arguments to construct offsets + * @param wires Arguments to construct wires + */ + template + explicit SparseHamiltonian(T1 &&data, T2 &&indices, T3 &&offsets, + T4 &&wires) + : BaseType{data, indices, offsets, wires} {} + + /** + * @brief Convenient wrapper for the constructor as the constructor does not + * convert the std::shared_ptr with a derived class correctly. + * + * This function is useful as std::make_shared does not handle + * brace-enclosed initializer list correctly. + * + * @param data Argument to construct data + * @param indices Argument to construct indices + * @param offsets Argument to construct ofsets + * @param wires Argument to construct wires + */ + static auto create(std::initializer_list data, + std::initializer_list indices, + std::initializer_list offsets, + std::initializer_list wires) + -> std::shared_ptr> { + return std::shared_ptr>( + new SparseHamiltonian{ + std::move(data), std::move(indices), std::move(offsets), + std::move(wires)}); + } + + /** + * @brief Updates the statevector SV:->SV', where SV' = a*H*SV, and where H + * is a sparse Hamiltonian. + * + */ + void applyInPlace(StateVectorT &sv) const override { + PL_ABORT_IF_NOT(this->wires_.size() == sv.getNumQubits(), + "SparseH wire count does not match state-vector size"); + StateVectorT d_sv_prime(sv.getNumQubits()); + + SparseMV_Kokkos( + sv.getView(), d_sv_prime.getView(), this->offsets_.data(), + this->offsets_.size(), this->indices_.data(), this->data_.data(), + this->data_.size()); + + sv.updateData(d_sv_prime); + } +}; + /// @cond DEV namespace detail { using Pennylane::LightningKokkos::Util::axpy_Kokkos; diff --git a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/tests/Test_ObservablesKokkos.cpp b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/tests/Test_ObservablesKokkos.cpp index aa616d2e33..5d402b1e2d 100644 --- a/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/tests/Test_ObservablesKokkos.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_kokkos/observables/tests/Test_ObservablesKokkos.cpp @@ -153,6 +153,20 @@ TEMPLATE_PRODUCT_TEST_CASE("Hamiltonian", "[Observables]", (StateVectorKokkos), } } +TEMPLATE_PRODUCT_TEST_CASE("SparseHamiltonian", "[Observables]", + (StateVectorKokkos), (float, double)) { + using StateVectorT = TestType; + using SparseHamiltonianT = SparseHamiltonian; + + SECTION("Copy constructibility") { + REQUIRE(std::is_copy_constructible_v); + } + + SECTION("Move constructibility") { + REQUIRE(std::is_move_constructible_v); + } +} + TEMPLATE_PRODUCT_TEST_CASE("Hamiltonian::ApplyInPlace", "[Observables]", (StateVectorKokkos), (float, double)) { using StateVectorT = TestType; diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/bindings/LQubitBindings.hpp b/pennylane_lightning/core/src/simulators/lightning_qubit/bindings/LQubitBindings.hpp index 190a9c6525..2f83f2f39d 100644 --- a/pennylane_lightning/core/src/simulators/lightning_qubit/bindings/LQubitBindings.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_qubit/bindings/LQubitBindings.hpp @@ -25,15 +25,17 @@ #include "DynamicDispatcher.hpp" #include "GateOperation.hpp" #include "MeasurementsLQubit.hpp" +#include "ObservablesLQubit.hpp" #include "StateVectorLQubitRaw.hpp" #include "TypeList.hpp" #include "VectorJacobianProduct.hpp" /// @cond DEV namespace { -using namespace Pennylane::LightningQubit::Measures; -using namespace Pennylane::LightningQubit::Algorithms; using namespace Pennylane::Bindings; +using namespace Pennylane::LightningQubit::Algorithms; +using namespace Pennylane::LightningQubit::Measures; +using namespace Pennylane::LightningQubit::Observables; using Pennylane::LightningQubit::StateVectorLQubitRaw; } // namespace /// @endcond @@ -180,6 +182,58 @@ void registerBackendSpecificMeasurements(PyClass &pyclass) { }); } +/** + * @brief Register backend specific observables. + * + * @tparam StateVectorT + * @param m Pybind module + */ +template +void registerBackendSpecificObservables([[maybe_unused]] py::module_ &m) { + using PrecisionT = + typename StateVectorT::PrecisionT; // Statevector's precision. + using ParamT = PrecisionT; // Parameter's data precision + + const std::string bitsize = + std::to_string(sizeof(std::complex) * 8); + + using np_arr_c = py::array_t, py::array::c_style>; + + std::string class_name; + + class_name = "SparseHamiltonianC" + bitsize; + py::class_, + std::shared_ptr>, + Observable>(m, class_name.c_str(), + py::module_local()) + .def(py::init([](const np_arr_c &data, + const std::vector &indices, + const std::vector &indptr, + const std::vector &wires) { + using ComplexT = typename StateVectorT::ComplexT; + const py::buffer_info buffer_data = data.request(); + const auto *data_ptr = static_cast(buffer_data.ptr); + + return SparseHamiltonian{ + std::vector({data_ptr, data_ptr + data.size()}), + indices, indptr, wires}; + })) + .def("__repr__", &SparseHamiltonian::getObsName) + .def("get_wires", &SparseHamiltonian::getWires, + "Get wires of observables") + .def( + "__eq__", + [](const SparseHamiltonian &self, + py::handle other) -> bool { + if (!py::isinstance>(other)) { + return false; + } + auto other_cast = other.cast>(); + return self == other_cast; + }, + "Compare two observables"); +} + /** * @brief Register Vector Jacobian Product. */ diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.cpp b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.cpp index 0a1fb54c48..e45e2c1572 100644 --- a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.cpp @@ -41,3 +41,9 @@ template class Observables::Hamiltonian>; template class Observables::Hamiltonian>; template class Observables::Hamiltonian>; + +template class Observables::SparseHamiltonian>; +template class Observables::SparseHamiltonian>; + +template class Observables::SparseHamiltonian>; +template class Observables::SparseHamiltonian>; diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp index 3433a3fcc3..b659969037 100644 --- a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp +++ b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/ObservablesLQubit.hpp @@ -27,6 +27,7 @@ #include "LinearAlgebra.hpp" // scaleAndAdd #include "Macros.hpp" // use_openmp #include "Observables.hpp" +#include "SparseLinAlg.hpp" #include "StateVectorLQubitManaged.hpp" #include "StateVectorLQubitRaw.hpp" #include "Util.hpp" @@ -359,4 +360,74 @@ class Hamiltonian final : public HamiltonianBase { } }; +/** + * @brief Sparse representation of Hamiltonian + * + */ +template +class SparseHamiltonian final : public SparseHamiltonianBase { + private: + using BaseType = SparseHamiltonianBase; + + public: + using PrecisionT = typename StateVectorT::PrecisionT; + using ComplexT = typename StateVectorT::ComplexT; + using IdxT = typename BaseType::IdxT; + + /** + * @brief Create a SparseHamiltonian from data, indices and offsets in CSR + * format. + * + * @param data Arguments to construct data + * @param indices Arguments to construct indices + * @param offsets Arguments to construct offsets + * @param wires Arguments to construct wires + */ + template + explicit SparseHamiltonian(T1 &&data, T2 &&indices, T3 &&offsets, + T4 &&wires) + : BaseType{data, indices, offsets, wires} {} + + /** + * @brief Convenient wrapper for the constructor as the constructor does not + * convert the std::shared_ptr with a derived class correctly. + * + * This function is useful as std::make_shared does not handle + * brace-enclosed initializer list correctly. + * + * @param data Argument to construct data + * @param indices Argument to construct indices + * @param offsets Argument to construct ofsets + * @param wires Argument to construct wires + */ + static auto create(std::initializer_list data, + std::initializer_list indices, + std::initializer_list offsets, + std::initializer_list wires) + -> std::shared_ptr> { + // NOLINTBEGIN(*-move-const-arg) + return std::shared_ptr>( + new SparseHamiltonian{ + std::move(data), std::move(indices), std::move(offsets), + std::move(wires)}); + // NOLINTEND(*-move-const-arg) + } + + /** + * @brief Updates the statevector SV:->SV', where SV' = a*H*SV, and where H + * is a sparse Hamiltonian. + * + */ + void applyInPlace(StateVectorT &sv) const override { + PL_ABORT_IF_NOT(this->wires_.size() == sv.getNumQubits(), + "SparseH wire count does not match state-vector size"); + auto operator_vector = Util::apply_Sparse_Matrix( + sv.getData(), sv.getLength(), this->offsets_.data(), + this->offsets_.size(), this->indices_.data(), this->data_.data(), + this->data_.size()); + + sv.updateData(operator_vector); + } +}; + } // namespace Pennylane::LightningQubit::Observables \ No newline at end of file diff --git a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/tests/Test_ObservablesLQubit.cpp b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/tests/Test_ObservablesLQubit.cpp index 624248da20..4ef59b04c5 100644 --- a/pennylane_lightning/core/src/simulators/lightning_qubit/observables/tests/Test_ObservablesLQubit.cpp +++ b/pennylane_lightning/core/src/simulators/lightning_qubit/observables/tests/Test_ObservablesLQubit.cpp @@ -159,6 +159,21 @@ TEMPLATE_PRODUCT_TEST_CASE("Hamiltonian", "[Observables]", } } +TEMPLATE_PRODUCT_TEST_CASE("SparseHamiltonian", "[Observables]", + (StateVectorLQubitManaged, StateVectorLQubitRaw), + (float, double)) { + using StateVectorT = TestType; + using SparseHamiltonianT = SparseHamiltonian; + + SECTION("Copy constructibility") { + REQUIRE(std::is_copy_constructible_v); + } + + SECTION("Move constructibility") { + REQUIRE(std::is_move_constructible_v); + } +} + TEMPLATE_PRODUCT_TEST_CASE("Hamiltonian::ApplyInPlace", "[Observables]", (StateVectorLQubitManaged, StateVectorLQubitRaw), (float, double)) { diff --git a/pennylane_lightning/core/src/utils/TestHelpers.hpp b/pennylane_lightning/core/src/utils/TestHelpers.hpp index 605e584bcb..af1a46618a 100644 --- a/pennylane_lightning/core/src/utils/TestHelpers.hpp +++ b/pennylane_lightning/core/src/utils/TestHelpers.hpp @@ -429,7 +429,8 @@ void write_CSR_vectors(std::vector &row_map, const ComplexT SC_ONE = 1.0; row_map.resize(numRows + 1); - for (IndexT rowIdx = 1; rowIdx < (IndexT)row_map.size(); ++rowIdx) { + for (IndexT rowIdx = 1; rowIdx < static_cast(row_map.size()); + ++rowIdx) { row_map[rowIdx] = row_map[rowIdx - 1] + 3; }; const IndexT numNNZ = row_map[numRows]; @@ -437,6 +438,7 @@ void write_CSR_vectors(std::vector &row_map, entries.resize(numNNZ); values.resize(numNNZ); for (IndexT rowIdx = 0; rowIdx < numRows; ++rowIdx) { + size_t idx = row_map[rowIdx]; if (rowIdx == 0) { entries[0] = rowIdx; entries[1] = rowIdx + 1; @@ -446,21 +448,21 @@ void write_CSR_vectors(std::vector &row_map, values[1] = -SC_ONE; values[2] = -SC_ONE; } else if (rowIdx == numRows - 1) { - entries[row_map[rowIdx]] = 0; - entries[row_map[rowIdx] + 1] = rowIdx - 1; - entries[row_map[rowIdx] + 2] = rowIdx; + entries[idx] = 0; + entries[idx + 1] = rowIdx - 1; + entries[idx + 2] = rowIdx; - values[row_map[rowIdx]] = -SC_ONE; - values[row_map[rowIdx] + 1] = -SC_ONE; - values[row_map[rowIdx] + 2] = SC_ONE; + values[idx] = -SC_ONE; + values[idx + 1] = -SC_ONE; + values[idx + 2] = SC_ONE; } else { - entries[row_map[rowIdx]] = rowIdx - 1; - entries[row_map[rowIdx] + 1] = rowIdx; - entries[row_map[rowIdx] + 2] = rowIdx + 1; + entries[idx] = rowIdx - 1; + entries[idx + 1] = rowIdx; + entries[idx + 2] = rowIdx + 1; - values[row_map[rowIdx]] = -SC_ONE; - values[row_map[rowIdx] + 1] = SC_ONE; - values[row_map[rowIdx] + 2] = -SC_ONE; + values[idx] = -SC_ONE; + values[idx + 1] = SC_ONE; + values[idx + 2] = -SC_ONE; } } }; diff --git a/tests/test_adjoint_jacobian.py b/tests/test_adjoint_jacobian.py index c8d65dd093..5e685d8d74 100644 --- a/tests/test_adjoint_jacobian.py +++ b/tests/test_adjoint_jacobian.py @@ -25,6 +25,7 @@ from pennylane import QNode, qnode from pennylane import qchem + I, X, Y, Z = ( np.eye(2), qml.PauliX.compute_matrix(), @@ -892,7 +893,10 @@ def circuit_ansatz(params, wires): qml.RX(params[29], wires=wires[1]) -@pytest.mark.skipif(not ld._CPP_BINARY_AVAILABLE, reason="Lightning binary required") +@pytest.mark.skipif( + device_name != "lightning.gpu" or not ld._CPP_BINARY_AVAILABLE, + reason="Lightning binary required", +) def test_tape_qchem(tol): """Tests the circuit Ansatz with a QChem Hamiltonian produces correct results""" @@ -944,6 +948,60 @@ def circuit(params): assert np.allclose(qml.grad(circuit_ld)(params), qml.grad(circuit_dq)(params), tol) +custom_wires = ["alice", 3.14, -1, 0] + + +@pytest.mark.skipif(not ld._CPP_BINARY_AVAILABLE, reason="Lightning binary required") +@pytest.mark.parametrize( + "returns", + [ + qml.SparseHamiltonian( + qml.Hamiltonian( + [0.1], + [qml.PauliX(wires=custom_wires[0]) @ qml.PauliZ(wires=custom_wires[1])], + ).sparse_matrix(custom_wires), + wires=custom_wires, + ), + qml.SparseHamiltonian( + qml.Hamiltonian( + [2.0], + [qml.PauliX(wires=custom_wires[2]) @ qml.PauliZ(wires=custom_wires[0])], + ).sparse_matrix(custom_wires), + wires=custom_wires, + ), + qml.SparseHamiltonian( + qml.Hamiltonian( + [1.1], + [qml.PauliX(wires=custom_wires[0]) @ qml.PauliZ(wires=custom_wires[2])], + ).sparse_matrix(custom_wires), + wires=custom_wires, + ), + ], +) +def test_adjoint_SparseHamiltonian(returns): + """Integration tests that compare to default.qubit for a large circuit containing parametrized + operations and when using custom wire labels""" + + dev = qml.device(device_name, wires=custom_wires) + dev_default = qml.device("default.qubit", wires=custom_wires) + + def circuit(params): + circuit_ansatz(params, wires=custom_wires) + return qml.expval(returns) + + n_params = 30 + np.random.seed(1337) + params = np.random.rand(n_params) + + qnode = qml.QNode(circuit, dev, diff_method="adjoint") + qnode_default = qml.QNode(circuit, dev_default, diff_method="parameter-shift") + + j_device = qml.jacobian(qnode)(params) + j_default = qml.jacobian(qnode_default)(params) + + assert np.allclose(j_device, j_default) + + @pytest.mark.parametrize( "returns", [ @@ -1024,9 +1082,6 @@ def casted_to_array_batched(params): assert np.allclose(j_def, j_lightning_batched) -custom_wires = ["alice", 3.14, -1, 0] - - @pytest.mark.parametrize( "returns", [ diff --git a/tests/test_serialize.py b/tests/test_serialize.py index ab6df3c26e..360ac0e71b 100644 --- a/tests/test_serialize.py +++ b/tests/test_serialize.py @@ -34,6 +34,8 @@ TensorProdObsC128, HamiltonianC64, HamiltonianC128, + SparseHamiltonianC64, + SparseHamiltonianC128, ) elif device_name == "lightning.gpu": from pennylane_lightning.lightning_gpu_ops.observables import ( @@ -45,6 +47,8 @@ TensorProdObsC128, HamiltonianC64, HamiltonianC128, + SparseHamiltonianC64, + SparseHamiltonianC128, ) else: from pennylane_lightning.lightning_qubit_ops.observables import ( @@ -56,6 +60,8 @@ TensorProdObsC128, HamiltonianC64, HamiltonianC128, + SparseHamiltonianC64, + SparseHamiltonianC128, ) @@ -92,6 +98,10 @@ def test_wrong_device_name(): (qml.Projector([0], wires=0), HermitianObsC128), (qml.Hamiltonian([1], [qml.PauliZ(0)]), HamiltonianC128), (qml.sum(qml.Hadamard(0), qml.PauliX(1)), HermitianObsC128), + ( + qml.SparseHamiltonian(qml.Hamiltonian([1], [qml.PauliZ(0)]).sparse_matrix(), wires=[0]), + SparseHamiltonianC128, + ), ], ) def test_obs_returns_expected_type(obs, obs_type):