diff --git a/README.md b/README.md index 75da61d8..04c59a4d 100644 --- a/README.md +++ b/README.md @@ -110,10 +110,10 @@ n_qubits = 3 state = StateVector.Haar_random_state(n_qubits, 0) circuit = Circuit(n_qubits) -circuit.add_gate(gate::X(0)) -circuit.add_gate(gate::CNot(0, 1)) -circuit.add_gate(gate::Y(1)) -circuit.add_gate(gate::RX(1, math.pi / 2)) +circuit.add_gate(gate.X(0)) +circuit.add_gate(gate.CNot(0, 1)) +circuit.add_gate(gate.Y(1)) +circuit.add_gate(gate.RX(1, math.pi / 2)) circuit.update_quantum_state(state) observable = Operator(n_qubits) diff --git a/exe/main.cpp b/exe/main.cpp index 29bf8f06..5659677e 100644 --- a/exe/main.cpp +++ b/exe/main.cpp @@ -80,11 +80,6 @@ int main() { std::cout << Json(gate::Pauli(pauli)) << std::endl; std::cout << Json(gate::PauliRotation(pauli, 0.5)) << std::endl; - std::cout << Json(gate::OneTargetMatrix(2, {0, 1, 2, 3})) << std::endl; - std::cout << Json(gate::TwoTargetMatrix( - 2, 3, {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15})) - << std::endl; - auto probgate = gate::Probablistic({.1, .9}, {gate::X(0), gate::I()}); std::cout << Json(probgate) << std::endl; @@ -119,13 +114,6 @@ int main() { Gate gate = j; std::cout << gate << std::endl; } - { - auto x = gate::OneTargetMatrix(1, {2., 3., 0., 10.}, {0, 3}); - Json j = x; - std::cout << j << std::endl; - Gate gate = j; - std::cout << gate << std::endl; - } { PauliOperator pauli("X 2 Y 1"); auto x = gate::PauliRotation(pauli, 1.5, {0, 3}); diff --git a/include/scaluq/gate/gate.hpp b/include/scaluq/gate/gate.hpp index 967273c0..672bd98f 100644 --- a/include/scaluq/gate/gate.hpp +++ b/include/scaluq/gate/gate.hpp @@ -1,6 +1,7 @@ #pragma once #include "../state/state_vector.hpp" +#include "../state/state_vector_batched.hpp" #include "../types.hpp" #include "../util/utility.hpp" @@ -56,12 +57,8 @@ class U2GateImpl; template class U3GateImpl; template -class OneTargetMatrixGateImpl; -template class SwapGateImpl; template -class TwoTargetMatrixGateImpl; -template class PauliGateImpl; template class PauliRotationGateImpl; @@ -98,9 +95,7 @@ enum class GateType { U1, U2, U3, - OneTargetMatrix, Swap, - TwoTargetMatrix, Pauli, PauliRotation, SparseMatrix, @@ -157,12 +152,8 @@ constexpr GateType get_gate_type() { return GateType::U2; else if constexpr (std::is_same_v>) return GateType::U3; - else if constexpr (std::is_same_v>) - return GateType::OneTargetMatrix; else if constexpr (std::is_same_v>) return GateType::Swap; - else if constexpr (std::is_same_v>) - return GateType::TwoTargetMatrix; else if constexpr (std::is_same_v>) return GateType::Pauli; else if constexpr (std::is_same_v>) @@ -188,6 +179,7 @@ class GateBase : public std::enable_shared_from_this> { std::uint64_t _target_mask, _control_mask; void check_qubit_mask_within_bounds(const StateVector& state_vector) const; + void check_qubit_mask_within_bounds(const StateVectorBatched& states) const; std::string get_qubit_info_as_string(const std::string& indent) const; @@ -214,6 +206,7 @@ class GateBase : public std::enable_shared_from_this> { [[nodiscard]] virtual internal::ComplexMatrix get_matrix() const = 0; virtual void update_quantum_state(StateVector& state_vector) const = 0; + virtual void update_quantum_state(StateVectorBatched& states) const = 0; [[nodiscard]] virtual std::string to_string(const std::string& indent = "") const = 0; @@ -316,8 +309,6 @@ class GatePtr { else if (type == "U2") gate = get_from_json>(j); else if (type == "U3") gate = get_from_json>(j); else if (type == "Swap") gate = get_from_json>(j); - else if (type == "OneTargetMatrix") gate = get_from_json>(j); - else if (type == "TwoTargetMatrix") gate = get_from_json>(j); else if (type == "Pauli") gate = get_from_json>(j); else if (type == "PauliRotation") gate = get_from_json>(j); else if (type == "Probablistic") gate = get_from_json>(j); @@ -429,11 +420,11 @@ void bind_gate_gate_hpp_without_precision(nb::module_& m) { .value("U1", GateType::U1) .value("U2", GateType::U2) .value("U3", GateType::U3) - .value("OneTargetMatrix", GateType::OneTargetMatrix) .value("Swap", GateType::Swap) - .value("TwoTargetMatrix", GateType::TwoTargetMatrix) .value("Pauli", GateType::Pauli) - .value("PauliRotation", GateType::PauliRotation); + .value("PauliRotation", GateType::PauliRotation) + .value("SparseMatrix", GateType::SparseMatrix) + .value("DenseMatrix", GateType::DenseMatrix); } template diff --git a/include/scaluq/gate/gate_factory.hpp b/include/scaluq/gate/gate_factory.hpp index 3accdf8f..5fb79040 100644 --- a/include/scaluq/gate/gate_factory.hpp +++ b/include/scaluq/gate/gate_factory.hpp @@ -143,13 +143,6 @@ inline Gate U3(std::uint64_t target, internal::vector_to_mask({target}), internal::vector_to_mask(controls), theta, phi, lambda); } template -inline Gate OneTargetMatrix(std::uint64_t target, - const std::array, 2>, 2>& matrix, - const std::vector& controls = {}) { - return internal::GateFactory::create_gate>( - internal::vector_to_mask({target}), internal::vector_to_mask(controls), matrix); -} -template inline Gate CX(std::uint64_t control, std::uint64_t target) { return internal::GateFactory::create_gate>( internal::vector_to_mask({target}), internal::vector_to_mask({control})); @@ -178,14 +171,6 @@ inline Gate Swap(std::uint64_t target1, internal::vector_to_mask({target1, target2}), internal::vector_to_mask(controls)); } template -inline Gate TwoTargetMatrix(std::uint64_t target1, - std::uint64_t target2, - const std::array, 4>, 4>& matrix, - const std::vector& controls = {}) { - return internal::GateFactory::create_gate>( - internal::vector_to_mask({target1, target2}), internal::vector_to_mask(controls), matrix); -} -template inline Gate Pauli(const PauliOperator& pauli, const std::vector& controls = {}) { auto tar = pauli.target_qubit_list(); @@ -397,19 +382,6 @@ void bind_gate_gate_factory_hpp(nb::module_& mgate) { &gate::CCX, "Generate general Gate class instance of Toffoli.\n\n.. note:: Toffoli is an alias " "of CCX."); - mgate.def("OneTargetMatrix", - &gate::OneTargetMatrix, - "Generate general Gate class instance of OneTargetMatrix.", - "target"_a, - "matrix"_a, - "controls"_a = std::vector{}); - mgate.def("TwoTargetMatrix", - &gate::TwoTargetMatrix, - "Generate general Gate class instance of TwoTargetMatrix.", - "target1"_a, - "target2"_a, - "matrix"_a, - "controls"_a = std::vector{}); mgate.def("DenseMatrix", &gate::DenseMatrix, "Generate general Gate class instance of DenseMatrix.", diff --git a/include/scaluq/gate/gate_matrix.hpp b/include/scaluq/gate/gate_matrix.hpp index bdb4b63d..4de24bd6 100644 --- a/include/scaluq/gate/gate_matrix.hpp +++ b/include/scaluq/gate/gate_matrix.hpp @@ -10,105 +10,6 @@ namespace scaluq { namespace internal { - -template -class OneTargetMatrixGateImpl : public GateBase { - Matrix2x2 _matrix; - -public: - OneTargetMatrixGateImpl(std::uint64_t target_mask, - std::uint64_t control_mask, - const std::array, 2>, 2>& matrix) - : GateBase(target_mask, control_mask) { - _matrix[0][0] = matrix[0][0]; - _matrix[0][1] = matrix[0][1]; - _matrix[1][0] = matrix[1][0]; - _matrix[1][1] = matrix[1][1]; - } - - std::array, 2>, 2> matrix() const { - return {_matrix[0][0], _matrix[0][1], _matrix[1][0], _matrix[1][1]}; - } - - std::shared_ptr> get_inverse() const override { - return std::make_shared( - this->_target_mask, - this->_control_mask, - std::array, 2>, 2>{Kokkos::conj(_matrix[0][0]), - Kokkos::conj(_matrix[1][0]), - Kokkos::conj(_matrix[0][1]), - Kokkos::conj(_matrix[1][1])}); - } - internal::ComplexMatrix get_matrix() const override; - - void update_quantum_state(StateVector& state_vector) const override; - std::string to_string(const std::string& indent) const override; - - void get_as_json(Json& j) const override { - j = Json{{"type", "OneTargetMatrix"}, - {"target", this->target_qubit_list()}, - {"control", this->control_qubit_list()}, - {"matrix", - std::vector>>{{_matrix[0][0], _matrix[0][1]}, - {_matrix[1][0], _matrix[1][1]}}}}; - } -}; - -template -class TwoTargetMatrixGateImpl : public GateBase { - Matrix4x4 _matrix; - -public: - TwoTargetMatrixGateImpl(std::uint64_t target_mask, - std::uint64_t control_mask, - const std::array, 4>, 4>& matrix) - : GateBase(target_mask, control_mask) { - for (std::uint64_t i : std::views::iota(0, 4)) { - for (std::uint64_t j : std::views::iota(0, 4)) { - _matrix[i][j] = matrix[i][j]; - } - } - } - - std::array, 4>, 4> matrix() const { - std::array, 4>, 4> matrix; - for (std::uint64_t i : std::views::iota(0, 4)) { - for (std::uint64_t j : std::views::iota(0, 4)) { - matrix[i][j] = _matrix[i][j]; - } - } - return matrix; - } - - std::shared_ptr> get_inverse() const override { - std::array, 4>, 4> matrix_dag; - for (std::uint64_t i : std::views::iota(0, 4)) { - for (std::uint64_t j : std::views::iota(0, 4)) { - matrix_dag[i][j] = Kokkos::conj(_matrix[j][i]); - } - } - return std::make_shared( - this->_target_mask, this->_control_mask, matrix_dag); - } - internal::ComplexMatrix get_matrix() const override; - - void update_quantum_state(StateVector& state_vector) const override; - - std::string to_string(const std::string& indent) const override; - - void get_as_json(Json& j) const override { - j = Json{{"type", "TwoTargetMatrix"}, - {"target", this->target_qubit_list()}, - {"control", this->control_qubit_list()}, - {"matrix", - std::vector>>{ - {_matrix[0][0], _matrix[0][1], _matrix[0][2], _matrix[0][3]}, - {_matrix[1][0], _matrix[1][1], _matrix[1][2], _matrix[1][3]}, - {_matrix[2][0], _matrix[2][1], _matrix[2][2], _matrix[2][3]}, - {_matrix[3][0], _matrix[3][1], _matrix[3][2], _matrix[3][3]}}}}; - } -}; - template class DenseMatrixGateImpl : public GateBase { Matrix _matrix; @@ -127,6 +28,7 @@ class DenseMatrixGateImpl : public GateBase { ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -157,6 +59,7 @@ class SparseMatrixGateImpl : public GateBase { SparseComplexMatrix get_sparse_matrix() const { return get_matrix().sparseView(); } void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -170,74 +73,17 @@ class SparseMatrixGateImpl : public GateBase { } // namespace internal -template -using OneTargetMatrixGate = internal::GatePtr>; -template -using TwoTargetMatrixGate = internal::GatePtr>; template using SparseMatrixGate = internal::GatePtr>; template using DenseMatrixGate = internal::GatePtr>; -namespace internal { -#define DECLARE_GET_FROM_JSON_ONETARGETMATRIXGATE_WITH_TYPE(Type) \ - template <> \ - inline std::shared_ptr> get_from_json(const Json& j) { \ - auto targets = j.at("target").get>(); \ - auto controls = j.at("control").get>(); \ - auto matrix = j.at("matrix").get>>>(); \ - return std::make_shared>( \ - vector_to_mask(targets), \ - vector_to_mask(controls), \ - std::array, 2>, 2>{ \ - matrix[0][0], matrix[0][1], matrix[1][0], matrix[1][1]}); \ - } -DECLARE_GET_FROM_JSON_ONETARGETMATRIXGATE_WITH_TYPE(double) -DECLARE_GET_FROM_JSON_ONETARGETMATRIXGATE_WITH_TYPE(float) -#undef DECLARE_GET_FROM_JSON_ONETARGETMATRIXGATE_WITH_TYPE - -#define DECLARE_GET_FROM_JSON_TWOTARGETMATRIXGATE_WITH_TYPE(Type) \ - template <> \ - inline std::shared_ptr> get_from_json(const Json& j) { \ - auto targets = j.at("target").get>(); \ - auto controls = j.at("control").get>(); \ - auto matrix = j.at("matrix").get>>>(); \ - return std::make_shared>( \ - vector_to_mask(targets), \ - vector_to_mask(controls), \ - std::array, 4>, 4>{matrix[0][0], \ - matrix[0][1], \ - matrix[0][2], \ - matrix[0][3], \ - matrix[1][0], \ - matrix[1][1], \ - matrix[1][2], \ - matrix[1][3], \ - matrix[2][0], \ - matrix[2][1], \ - matrix[2][2], \ - matrix[2][3], \ - matrix[3][0], \ - matrix[3][1], \ - matrix[3][2], \ - matrix[3][3]}); \ - } -DECLARE_GET_FROM_JSON_TWOTARGETMATRIXGATE_WITH_TYPE(double) -DECLARE_GET_FROM_JSON_TWOTARGETMATRIXGATE_WITH_TYPE(float) -#undef DECLARE_GET_FROM_JSON_TWOTARGETMATRIXGATE_WITH_TYPE - -} // namespace internal - #ifdef SCALUQ_USE_NANOBIND namespace internal { template void bind_gate_gate_matrix_hpp(nb::module_& m) { - DEF_GATE(OneTargetMatrixGate, Fp, "Specific class of one-qubit dense matrix gate.") - .def("matrix", [](const OneTargetMatrixGate& gate) { return gate->matrix(); }); - DEF_GATE(TwoTargetMatrixGate, Fp, "Specific class of two-qubit dense matrix gate.") - .def("matrix", [](const TwoTargetMatrixGate& gate) { return gate->matrix(); }); DEF_GATE(SparseMatrixGate, Fp, "Specific class of sparse matrix gate.") - .def("matrix", [](const SparseMatrixGate& gate) { return gate->get_matrix(); }) + .def("matrix", [](const SparseMatrixGate& gate) { return gate->get_matrix(); }) .def("sparse_matrix", [](const SparseMatrixGate& gate) { return gate->get_sparse_matrix(); }); DEF_GATE(DenseMatrixGate, Fp, "Specific class of dense matrix gate.") diff --git a/include/scaluq/gate/gate_pauli.hpp b/include/scaluq/gate/gate_pauli.hpp index 6a4b649c..71157c63 100644 --- a/include/scaluq/gate/gate_pauli.hpp +++ b/include/scaluq/gate/gate_pauli.hpp @@ -27,6 +27,7 @@ class PauliGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override { return this->_pauli.get_matrix(); } void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -59,6 +60,7 @@ class PauliRotationGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; diff --git a/include/scaluq/gate/gate_probablistic.hpp b/include/scaluq/gate/gate_probablistic.hpp index 72bd8f6f..5d656ccd 100644 --- a/include/scaluq/gate/gate_probablistic.hpp +++ b/include/scaluq/gate/gate_probablistic.hpp @@ -11,7 +11,7 @@ namespace internal { template class ProbablisticGateImpl : public GateBase { std::vector _distribution; - std::vector _cumlative_distribution; + std::vector _cumulative_distribution; std::vector> _gate_list; public: @@ -59,6 +59,7 @@ class ProbablisticGateImpl : public GateBase { } void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; diff --git a/include/scaluq/gate/gate_standard.hpp b/include/scaluq/gate/gate_standard.hpp index a160be0f..50dc2668 100644 --- a/include/scaluq/gate/gate_standard.hpp +++ b/include/scaluq/gate/gate_standard.hpp @@ -17,6 +17,7 @@ class IGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -40,6 +41,7 @@ class GlobalPhaseGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -74,6 +76,7 @@ class XGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -95,6 +98,7 @@ class YGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -116,6 +120,7 @@ class ZGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -137,6 +142,7 @@ class HGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -175,6 +181,7 @@ class SGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -196,6 +203,7 @@ class SdagGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -217,6 +225,7 @@ class TGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -238,6 +247,7 @@ class TdagGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -261,6 +271,7 @@ class SqrtXGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -282,6 +293,7 @@ class SqrtXdagGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -305,6 +317,7 @@ class SqrtYGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -326,6 +339,7 @@ class SqrtYdagGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -347,6 +361,7 @@ class P0GateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -368,6 +383,7 @@ class P1GateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -390,6 +406,7 @@ class RXGateImpl : public RotationGateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -413,6 +430,7 @@ class RYGateImpl : public RotationGateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -436,6 +454,7 @@ class RZGateImpl : public RotationGateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -464,6 +483,7 @@ class U1GateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -494,6 +514,7 @@ class U2GateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -525,6 +546,7 @@ class U3GateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; @@ -549,6 +571,7 @@ class SwapGateImpl : public GateBase { internal::ComplexMatrix get_matrix() const override; void update_quantum_state(StateVector& state_vector) const override; + void update_quantum_state(StateVectorBatched& states) const override; std::string to_string(const std::string& indent) const override; diff --git a/include/scaluq/gate/param_gate.hpp b/include/scaluq/gate/param_gate.hpp index cfa0b0b9..cbdacb73 100644 --- a/include/scaluq/gate/param_gate.hpp +++ b/include/scaluq/gate/param_gate.hpp @@ -1,6 +1,7 @@ #pragma once #include "../state/state_vector.hpp" +#include "../state/state_vector_batched.hpp" #include "../types.hpp" namespace scaluq { @@ -62,6 +63,7 @@ class ParamGateBase { std::uint64_t _target_mask, _control_mask; Fp _pcoef; void check_qubit_mask_within_bounds(const StateVector& state_vector) const; + void check_qubit_mask_within_bounds(const StateVectorBatched& states) const; std::string get_qubit_info_as_string(const std::string& indent) const; @@ -90,6 +92,8 @@ class ParamGateBase { [[nodiscard]] virtual internal::ComplexMatrix get_matrix(Fp param) const = 0; virtual void update_quantum_state(StateVector& state_vector, Fp param) const = 0; + virtual void update_quantum_state(StateVectorBatched& states, + std::vector params) const = 0; [[nodiscard]] virtual std::string to_string(const std::string& indent = "") const = 0; @@ -230,6 +234,13 @@ namespace internal { FLOAT param) { param_gate->update_quantum_state(state_vector, param); }, \ "Apply gate to `state_vector` with holding the parameter. `state_vector` in args is " \ "directly updated.") \ + .def( \ + "update_quantum_state", \ + [](const PARAM_GATE_TYPE& param_gate, \ + StateVectorBatched& states, \ + std::vector params) { param_gate->update_quantum_state(states, params); }, \ + "Apply gate to `states` with holding the parameter. `states` in args is directly " \ + "updated.") \ .def( \ "get_matrix", \ [](const PARAM_GATE_TYPE& gate, FLOAT param) { \ diff --git a/include/scaluq/gate/param_gate_pauli.hpp b/include/scaluq/gate/param_gate_pauli.hpp index 2dce8ce6..ee59261e 100644 --- a/include/scaluq/gate/param_gate_pauli.hpp +++ b/include/scaluq/gate/param_gate_pauli.hpp @@ -29,6 +29,8 @@ class ParamPauliRotationGateImpl : public ParamGateBase { } internal::ComplexMatrix get_matrix(Fp param) const override; void update_quantum_state(StateVector& state_vector, Fp param) const override; + void update_quantum_state(StateVectorBatched& states, + std::vector params) const override; std::string to_string(const std::string& indent) const override; void get_as_json(Json& j) const override { diff --git a/include/scaluq/gate/param_gate_probablistic.hpp b/include/scaluq/gate/param_gate_probablistic.hpp index 182812f9..035f4be1 100644 --- a/include/scaluq/gate/param_gate_probablistic.hpp +++ b/include/scaluq/gate/param_gate_probablistic.hpp @@ -13,7 +13,7 @@ template class ParamProbablisticGateImpl : public ParamGateBase { using EitherGate = std::variant, ParamGate>; std::vector _distribution; - std::vector _cumlative_distribution; + std::vector _cumulative_distribution; std::vector _gate_list; public: @@ -61,6 +61,8 @@ class ParamProbablisticGateImpl : public ParamGateBase { } void update_quantum_state(StateVector& state_vector, Fp param) const override; + void update_quantum_state(StateVectorBatched& states, + std::vector params) const override; std::string to_string(const std::string& indent) const override; diff --git a/include/scaluq/gate/param_gate_standard.hpp b/include/scaluq/gate/param_gate_standard.hpp index 477ba5e9..4ac64da2 100644 --- a/include/scaluq/gate/param_gate_standard.hpp +++ b/include/scaluq/gate/param_gate_standard.hpp @@ -19,6 +19,8 @@ class ParamRXGateImpl : public ParamGateBase { internal::ComplexMatrix get_matrix(Fp param) const override; void update_quantum_state(StateVector& state_vector, Fp param) const override; + void update_quantum_state(StateVectorBatched& states, + std::vector params) const override; std::string to_string(const std::string& indent) const override; @@ -42,6 +44,8 @@ class ParamRYGateImpl : public ParamGateBase { internal::ComplexMatrix get_matrix(Fp param) const override; void update_quantum_state(StateVector& state_vector, Fp param) const override; + void update_quantum_state(StateVectorBatched& states, + std::vector params) const override; std::string to_string(const std::string& indent) const override; @@ -65,6 +69,8 @@ class ParamRZGateImpl : public ParamGateBase { internal::ComplexMatrix get_matrix(Fp param) const override; void update_quantum_state(StateVector& state_vector, Fp param) const override; + void update_quantum_state(StateVectorBatched& states, + std::vector params) const override; std::string to_string(const std::string& indent) const override; diff --git a/include/scaluq/state/state_vector.hpp b/include/scaluq/state/state_vector.hpp index 1239051f..a3d8292e 100644 --- a/include/scaluq/state/state_vector.hpp +++ b/include/scaluq/state/state_vector.hpp @@ -27,6 +27,7 @@ class StateVector { Kokkos::View _raw; StateVector() = default; StateVector(std::uint64_t n_qubits); + StateVector(Kokkos::View view); StateVector(const StateVector& other) = default; StateVector& operator=(const StateVector& other) = default; diff --git a/pyproject.toml b/pyproject.toml index 01e731f1..a3916a86 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,8 @@ readme = "README.md" requires-python = ">=3.8" license = { file = "LICENSE" } dependencies = [ - "numpy >= 1.22.0" + "numpy >= 1.22.0", + "scipy >= 1.11.0" ] classifiers = [ "Development Status :: 3 - Alpha", diff --git a/python/binding.cpp b/python/binding.cpp index 39aa414d..d9059a10 100644 --- a/python/binding.cpp +++ b/python/binding.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include diff --git a/src/gate/gate.cpp b/src/gate/gate.cpp index b8d9b3e1..b5552977 100644 --- a/src/gate/gate.cpp +++ b/src/gate/gate.cpp @@ -13,6 +13,15 @@ void GateBase::check_qubit_mask_within_bounds(const StateVector& state_v "Target/Control qubit exceeds the number of qubits in the system."); } } +FLOAT(Fp) +void GateBase::check_qubit_mask_within_bounds(const StateVectorBatched& states) const { + std::uint64_t full_mask = (1ULL << states.n_qubits()) - 1; + if ((_target_mask | _control_mask) > full_mask) [[unlikely]] { + throw std::runtime_error( + "Error: Gate::update_quantum_state(StateVectorBatched& states): " + "Target/Control qubit exceeds the number of qubits in the system."); + } +} FLOAT(Fp) std::string GateBase::get_qubit_info_as_string(const std::string& indent) const { diff --git a/src/gate/gate_matrix.cpp b/src/gate/gate_matrix.cpp index 11d05f46..384afccb 100644 --- a/src/gate/gate_matrix.cpp +++ b/src/gate/gate_matrix.cpp @@ -4,49 +4,6 @@ #include "update_ops.hpp" namespace scaluq::internal { -FLOAT(Fp) -ComplexMatrix OneTargetMatrixGateImpl::get_matrix() const { - internal::ComplexMatrix mat(2, 2); - mat << this->_matrix[0][0], this->_matrix[0][1], this->_matrix[1][0], this->_matrix[1][1]; - return mat; -} -FLOAT(Fp) -void OneTargetMatrixGateImpl::update_quantum_state(StateVector& state_vector) const { - this->check_qubit_mask_within_bounds(state_vector); - one_target_dense_matrix_gate(this->_target_mask, this->_control_mask, _matrix, state_vector); -} -FLOAT(Fp) -std::string OneTargetMatrixGateImpl::to_string(const std::string& indent) const { - std::ostringstream ss; - ss << indent << "Gate Type: OneTargetMatrix\n"; - ss << this->get_qubit_info_as_string(indent); - return ss.str(); -} -FLOAT_DECLARE_CLASS(OneTargetMatrixGateImpl) - -FLOAT(Fp) -ComplexMatrix TwoTargetMatrixGateImpl::get_matrix() const { - internal::ComplexMatrix mat(4, 4); - mat << this->_matrix[0][0], this->_matrix[0][1], this->_matrix[0][2], this->_matrix[0][3], - this->_matrix[1][0], this->_matrix[1][1], this->_matrix[1][2], this->_matrix[1][3], - this->_matrix[2][0], this->_matrix[2][1], this->_matrix[2][2], this->_matrix[2][3], - this->_matrix[3][0], this->_matrix[3][1], this->_matrix[3][2], this->_matrix[3][3]; - return mat; -} -FLOAT(Fp) -void TwoTargetMatrixGateImpl::update_quantum_state(StateVector& state_vector) const { - this->check_qubit_mask_within_bounds(state_vector); - two_target_dense_matrix_gate(this->_target_mask, this->_control_mask, _matrix, state_vector); -} -FLOAT(Fp) -std::string TwoTargetMatrixGateImpl::to_string(const std::string& indent) const { - std::ostringstream ss; - ss << indent << "Gate Type: TwoTargetMatrix\n"; - ss << this->get_qubit_info_as_string(indent); - return ss.str(); -} -FLOAT_DECLARE_CLASS(TwoTargetMatrixGateImpl) - FLOAT(Fp) DenseMatrixGateImpl::DenseMatrixGateImpl(std::uint64_t target_mask, std::uint64_t control_mask, @@ -83,6 +40,11 @@ void DenseMatrixGateImpl::update_quantum_state(StateVector& state_vector dense_matrix_gate(this->_target_mask, this->_control_mask, _matrix, state_vector); } FLOAT(Fp) +void DenseMatrixGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + dense_matrix_gate(this->_target_mask, this->_control_mask, _matrix, states); +} +FLOAT(Fp) std::string DenseMatrixGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: DenseMatrix\n"; @@ -128,6 +90,11 @@ void SparseMatrixGateImpl::update_quantum_state(StateVector& state_vecto sparse_matrix_gate(this->_target_mask, this->_control_mask, _matrix, state_vector); } FLOAT(Fp) +void SparseMatrixGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + sparse_matrix_gate(this->_target_mask, this->_control_mask, _matrix, states); +} +FLOAT(Fp) std::string SparseMatrixGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: SparseMatrix\n"; diff --git a/src/gate/gate_pauli.cpp b/src/gate/gate_pauli.cpp index d059d611..9a7aaef4 100644 --- a/src/gate/gate_pauli.cpp +++ b/src/gate/gate_pauli.cpp @@ -10,6 +10,11 @@ void PauliGateImpl::update_quantum_state(StateVector& state_vector) cons apply_pauli(this->_control_mask, bit_flip_mask, phase_flip_mask, _pauli.coef(), state_vector); } FLOAT(Fp) +void PauliGateImpl::update_quantum_state(StateVectorBatched& states) const { + auto [bit_flip_mask, phase_flip_mask] = _pauli.get_XZ_mask_representation(); + apply_pauli(this->_control_mask, bit_flip_mask, phase_flip_mask, _pauli.coef(), states); +} +FLOAT(Fp) std::string PauliGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; auto controls = this->control_qubit_list(); @@ -40,6 +45,12 @@ void PauliRotationGateImpl::update_quantum_state(StateVector& state_vect this->_control_mask, bit_flip_mask, phase_flip_mask, _pauli.coef(), _angle, state_vector); } FLOAT(Fp) +void PauliRotationGateImpl::update_quantum_state(StateVectorBatched& states) const { + auto [bit_flip_mask, phase_flip_mask] = _pauli.get_XZ_mask_representation(); + apply_pauli_rotation( + this->_control_mask, bit_flip_mask, phase_flip_mask, _pauli.coef(), _angle, states); +} +FLOAT(Fp) std::string PauliRotationGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; auto controls = this->control_qubit_list(); diff --git a/src/gate/gate_probablistic.cpp b/src/gate/gate_probablistic.cpp index 6975bc3e..8f6b0f4d 100644 --- a/src/gate/gate_probablistic.cpp +++ b/src/gate/gate_probablistic.cpp @@ -14,9 +14,10 @@ ProbablisticGateImpl::ProbablisticGateImpl(const std::vector& distributi if (n != gate_list.size()) { throw std::runtime_error("distribution and gate_list have different size."); } - _cumlative_distribution.resize(n + 1); - std::partial_sum(distribution.begin(), distribution.end(), _cumlative_distribution.begin() + 1); - if (std::abs(_cumlative_distribution.back() - 1.) > 1e-6) { + _cumulative_distribution.resize(n + 1); + std::partial_sum( + distribution.begin(), distribution.end(), _cumulative_distribution.begin() + 1); + if (std::abs(_cumulative_distribution.back() - 1.) > 1e-6) { throw std::runtime_error("Sum of distribution must be equal to 1."); } } @@ -33,13 +34,33 @@ FLOAT(Fp) void ProbablisticGateImpl::update_quantum_state(StateVector& state_vector) const { Random random; Fp r = random.uniform(); - std::uint64_t i = std::distance(_cumlative_distribution.begin(), - std::ranges::upper_bound(_cumlative_distribution, r)) - + std::uint64_t i = std::distance(_cumulative_distribution.begin(), + std::ranges::upper_bound(_cumulative_distribution, r)) - 1; if (i >= _gate_list.size()) i = _gate_list.size() - 1; _gate_list[i]->update_quantum_state(state_vector); } FLOAT(Fp) +void ProbablisticGateImpl::update_quantum_state(StateVectorBatched& states) const { + std::vector indices(states.batch_size()); + std::vector r(states.batch_size()); + + Random random; + for (std::size_t i = 0; i < states.batch_size(); ++i) { + r[i] = random.uniform(); + indices[i] = std::distance(_cumulative_distribution.begin(), + std::ranges::upper_bound(_cumulative_distribution, r[i])) - + 1; + if (indices[i] >= _gate_list.size()) indices[i] = _gate_list.size() - 1; + auto state_vector = StateVector(Kokkos::subview(states._raw, i, Kokkos::ALL)); + _gate_list[indices[i]]->update_quantum_state(state_vector); + Kokkos::parallel_for( + "update_states", states.dim(), KOKKOS_CLASS_LAMBDA(const int j) { + states._raw(i, j) = state_vector._raw(j); + }); + } +} +FLOAT(Fp) std::string ProbablisticGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; const auto dist = distribution(); diff --git a/src/gate/gate_standard.cpp b/src/gate/gate_standard.cpp index 3a1b9aad..6e3845c3 100644 --- a/src/gate/gate_standard.cpp +++ b/src/gate/gate_standard.cpp @@ -13,6 +13,10 @@ void IGateImpl::update_quantum_state(StateVector& state_vector) const { i_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void IGateImpl::update_quantum_state(StateVectorBatched& states) const { + i_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string IGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: I\n"; @@ -31,6 +35,11 @@ void GlobalPhaseGateImpl::update_quantum_state(StateVector& state_vector global_phase_gate(this->_target_mask, this->_control_mask, _phase, state_vector); } FLOAT(Fp) +void GlobalPhaseGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + global_phase_gate(this->_target_mask, this->_control_mask, _phase, states); +} +FLOAT(Fp) std::string GlobalPhaseGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: GlobalPhase\n"; @@ -52,6 +61,11 @@ void XGateImpl::update_quantum_state(StateVector& state_vector) const { x_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void XGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + x_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string XGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: X\n"; @@ -72,6 +86,11 @@ void YGateImpl::update_quantum_state(StateVector& state_vector) const { y_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void YGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + y_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string YGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: Y\n"; @@ -92,6 +111,11 @@ void ZGateImpl::update_quantum_state(StateVector& state_vector) const { z_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void ZGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + z_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string ZGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: Z\n"; @@ -113,6 +137,11 @@ void HGateImpl::update_quantum_state(StateVector& state_vector) const { h_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void HGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + h_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string HGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: H\n"; @@ -133,6 +162,11 @@ void SGateImpl::update_quantum_state(StateVector& state_vector) const { s_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void SGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + s_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string SGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: S\n"; @@ -153,6 +187,11 @@ void SdagGateImpl::update_quantum_state(StateVector& state_vector) const sdag_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void SdagGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + sdag_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string SdagGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: Sdag\n"; @@ -173,6 +212,11 @@ void TGateImpl::update_quantum_state(StateVector& state_vector) const { t_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void TGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + t_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string TGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: T\n"; @@ -193,6 +237,11 @@ void TdagGateImpl::update_quantum_state(StateVector& state_vector) const tdag_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void TdagGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + tdag_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string TdagGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: Tdag\n"; @@ -214,6 +263,11 @@ void SqrtXGateImpl::update_quantum_state(StateVector& state_vector) cons sqrtx_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void SqrtXGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + sqrtx_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string SqrtXGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: SqrtX\n"; @@ -235,6 +289,11 @@ void SqrtXdagGateImpl::update_quantum_state(StateVector& state_vector) c sqrtxdag_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void SqrtXdagGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + sqrtxdag_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string SqrtXdagGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: SqrtXdag\n"; @@ -256,6 +315,11 @@ void SqrtYGateImpl::update_quantum_state(StateVector& state_vector) cons sqrty_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void SqrtYGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + sqrty_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string SqrtYGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: SqrtY\n"; @@ -277,6 +341,11 @@ void SqrtYdagGateImpl::update_quantum_state(StateVector& state_vector) c sqrtydag_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void SqrtYdagGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + sqrtydag_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string SqrtYdagGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: SqrtYdag\n"; @@ -297,6 +366,11 @@ void P0GateImpl::update_quantum_state(StateVector& state_vector) const { p0_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void P0GateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + p0_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string P0GateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: P0\n"; @@ -317,6 +391,11 @@ void P1GateImpl::update_quantum_state(StateVector& state_vector) const { p1_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void P1GateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + p1_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string P1GateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: P1\n"; @@ -338,6 +417,11 @@ void RXGateImpl::update_quantum_state(StateVector& state_vector) const { rx_gate(this->_target_mask, this->_control_mask, this->_angle, state_vector); } FLOAT(Fp) +void RXGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + rx_gate(this->_target_mask, this->_control_mask, this->_angle, states); +} +FLOAT(Fp) std::string RXGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: RX\n"; @@ -360,6 +444,11 @@ void RYGateImpl::update_quantum_state(StateVector& state_vector) const { ry_gate(this->_target_mask, this->_control_mask, this->_angle, state_vector); } FLOAT(Fp) +void RYGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + ry_gate(this->_target_mask, this->_control_mask, this->_angle, states); +} +FLOAT(Fp) std::string RYGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: RY\n"; @@ -382,6 +471,11 @@ void RZGateImpl::update_quantum_state(StateVector& state_vector) const { rz_gate(this->_target_mask, this->_control_mask, this->_angle, state_vector); } FLOAT(Fp) +void RZGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + rz_gate(this->_target_mask, this->_control_mask, this->_angle, states); +} +FLOAT(Fp) std::string RZGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: RZ\n"; @@ -402,6 +496,11 @@ void U1GateImpl::update_quantum_state(StateVector& state_vector) const { u1_gate(this->_target_mask, this->_control_mask, _lambda, state_vector); } FLOAT(Fp) +void U1GateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + u1_gate(this->_target_mask, this->_control_mask, _lambda, states); +} +FLOAT(Fp) std::string U1GateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: U1\n"; @@ -426,6 +525,11 @@ void U2GateImpl::update_quantum_state(StateVector& state_vector) const { u2_gate(this->_target_mask, this->_control_mask, _phi, _lambda, state_vector); } FLOAT(Fp) +void U2GateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + u2_gate(this->_target_mask, this->_control_mask, _phi, _lambda, states); +} +FLOAT(Fp) std::string U2GateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: U2\n"; @@ -449,6 +553,11 @@ void U3GateImpl::update_quantum_state(StateVector& state_vector) const { u3_gate(this->_target_mask, this->_control_mask, _theta, _phi, _lambda, state_vector); } FLOAT(Fp) +void U3GateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + u3_gate(this->_target_mask, this->_control_mask, _theta, _phi, _lambda, states); +} +FLOAT(Fp) std::string U3GateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: U3\n"; @@ -469,6 +578,11 @@ void SwapGateImpl::update_quantum_state(StateVector& state_vector) const swap_gate(this->_target_mask, this->_control_mask, state_vector); } FLOAT(Fp) +void SwapGateImpl::update_quantum_state(StateVectorBatched& states) const { + this->check_qubit_mask_within_bounds(states); + swap_gate(this->_target_mask, this->_control_mask, states); +} +FLOAT(Fp) std::string SwapGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: Swap\n"; diff --git a/src/gate/param_gate.cpp b/src/gate/param_gate.cpp index 920943cf..5174e546 100644 --- a/src/gate/param_gate.cpp +++ b/src/gate/param_gate.cpp @@ -13,6 +13,15 @@ void ParamGateBase::check_qubit_mask_within_bounds(const StateVector& st } } FLOAT(Fp) +void ParamGateBase::check_qubit_mask_within_bounds(const StateVectorBatched& states) const { + std::uint64_t full_mask = (1ULL << states.n_qubits()) - 1; + if ((_target_mask | _control_mask) > full_mask) [[unlikely]] { + throw std::runtime_error( + "Error: ParamGate::update_quantum_state(StateVector& state): " + "Target/Control qubit exceeds the number of qubits in the system."); + } +} +FLOAT(Fp) std::string ParamGateBase::get_qubit_info_as_string(const std::string& indent) const { std::ostringstream ss; auto targets = target_qubit_list(); diff --git a/src/gate/param_gate_pauli.cpp b/src/gate/param_gate_pauli.cpp index bfd51019..7938a8db 100644 --- a/src/gate/param_gate_pauli.cpp +++ b/src/gate/param_gate_pauli.cpp @@ -27,6 +27,18 @@ void ParamPauliRotationGateImpl::update_quantum_state(StateVector& state state_vector); } FLOAT(Fp) +void ParamPauliRotationGateImpl::update_quantum_state(StateVectorBatched& states, + std::vector params) const { + auto [bit_flip_mask, phase_flip_mask] = _pauli.get_XZ_mask_representation(); + apply_pauli_rotation(this->_control_mask, + bit_flip_mask, + phase_flip_mask, + _pauli.coef(), + this->_pcoef, + params, + states); +} +FLOAT(Fp) std::string ParamPauliRotationGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; auto controls = this->control_qubit_list(); diff --git a/src/gate/param_gate_probablistic.cpp b/src/gate/param_gate_probablistic.cpp index a1eb314b..62301345 100644 --- a/src/gate/param_gate_probablistic.cpp +++ b/src/gate/param_gate_probablistic.cpp @@ -15,9 +15,10 @@ ParamProbablisticGateImpl::ParamProbablisticGateImpl( if (n != gate_list.size()) { throw std::runtime_error("distribution and gate_list have different size."); } - _cumlative_distribution.resize(n + 1); - std::partial_sum(distribution.begin(), distribution.end(), _cumlative_distribution.begin() + 1); - if (std::abs(_cumlative_distribution.back() - 1.) > 1e-6) { + _cumulative_distribution.resize(n + 1); + std::partial_sum( + distribution.begin(), distribution.end(), _cumulative_distribution.begin() + 1); + if (std::abs(_cumulative_distribution.back() - 1.) > 1e-6) { throw std::runtime_error("Sum of distribution must be equal to 1."); } } @@ -36,8 +37,8 @@ void ParamProbablisticGateImpl::update_quantum_state(StateVector& state_ Fp param) const { Random random; Fp r = random.uniform(); - std::uint64_t i = std::distance(_cumlative_distribution.begin(), - std::ranges::upper_bound(_cumlative_distribution, r)) - + std::uint64_t i = std::distance(_cumulative_distribution.begin(), + std::ranges::upper_bound(_cumulative_distribution, r)) - 1; if (i >= _gate_list.size()) i = _gate_list.size() - 1; const auto& gate = _gate_list[i]; @@ -48,6 +49,32 @@ void ParamProbablisticGateImpl::update_quantum_state(StateVector& state_ } } FLOAT(Fp) +void ParamProbablisticGateImpl::update_quantum_state(StateVectorBatched& states, + std::vector params) const { + Random random; + std::vector r(states.batch_size()); + std::ranges::generate(r, [&random]() { return random.uniform(); }); + std::vector indicies(states.batch_size()); + std::ranges::transform(r, indicies.begin(), [this](double r) { + return std::distance(_cumulative_distribution.begin(), + std::ranges::upper_bound(_cumulative_distribution, r)) - + 1; + }); + std::ranges::transform(indicies, indicies.begin(), [this](std::uint64_t i) { + if (i >= _gate_list.size()) i = _gate_list.size() - 1; + return i; + }); + for (std::size_t i = 0; i < states.batch_size(); ++i) { + const auto& gate = _gate_list[indicies[i]]; + auto state_vector = StateVector(Kokkos::subview(states._raw, i, Kokkos::ALL)); + if (gate.index() == 0) { + std::get<0>(gate)->update_quantum_state(state_vector); + } else { + std::get<1>(gate)->update_quantum_state(state_vector, params[i]); + } + } +} +FLOAT(Fp) std::string ParamProbablisticGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; const auto dist = distribution(); diff --git a/src/gate/param_gate_standard.cpp b/src/gate/param_gate_standard.cpp index 032ef116..149b2171 100644 --- a/src/gate/param_gate_standard.cpp +++ b/src/gate/param_gate_standard.cpp @@ -18,6 +18,12 @@ void ParamRXGateImpl::update_quantum_state(StateVector& state_vector, Fp rx_gate(this->_target_mask, this->_control_mask, this->_pcoef * param, state_vector); } FLOAT(Fp) +void ParamRXGateImpl::update_quantum_state(StateVectorBatched& states, + std::vector params) const { + this->check_qubit_mask_within_bounds(states); + rx_gate(this->_target_mask, this->_control_mask, this->_pcoef, params, states); +} +FLOAT(Fp) std::string ParamRXGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: ParamRX\n"; @@ -39,6 +45,12 @@ void ParamRYGateImpl::update_quantum_state(StateVector& state_vector, Fp ry_gate(this->_target_mask, this->_control_mask, this->_pcoef * param, state_vector); } FLOAT(Fp) +void ParamRYGateImpl::update_quantum_state(StateVectorBatched& states, + std::vector params) const { + this->check_qubit_mask_within_bounds(states); + ry_gate(this->_target_mask, this->_control_mask, this->_pcoef, params, states); +} +FLOAT(Fp) std::string ParamRYGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: ParamRY\n"; @@ -60,6 +72,12 @@ void ParamRZGateImpl::update_quantum_state(StateVector& state_vector, Fp rz_gate(this->_target_mask, this->_control_mask, this->_pcoef * param, state_vector); } FLOAT(Fp) +void ParamRZGateImpl::update_quantum_state(StateVectorBatched& states, + std::vector params) const { + this->check_qubit_mask_within_bounds(states); + rz_gate(this->_target_mask, this->_control_mask, this->_pcoef, params, states); +} +FLOAT(Fp) std::string ParamRZGateImpl::to_string(const std::string& indent) const { std::ostringstream ss; ss << indent << "Gate Type: ParamRZ\n"; diff --git a/src/gate/update_ops.hpp b/src/gate/update_ops.hpp index c977da52..cbc24c03 100644 --- a/src/gate/update_ops.hpp +++ b/src/gate/update_ops.hpp @@ -2,6 +2,7 @@ #include #include +#include #include namespace scaluq { @@ -11,48 +12,87 @@ template void none_target_dense_matrix_gate(std::uint64_t control_mask, const Matrix& matrix, StateVector& state); +template +void none_target_dense_matrix_gate(std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states); template void one_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const Matrix2x2& matrix, StateVector& state); +template +void one_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix2x2& matrix, + StateVectorBatched& states); template void two_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const Matrix4x4& matrix, StateVector& state); +template +void two_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix4x4& matrix, + StateVectorBatched& states); template void single_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const Matrix& matrix, StateVector& state); +template +void single_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states); template void double_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const Matrix& matrix, StateVector& state); +template +void double_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states); template void multi_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const Matrix& matrix, StateVector& state); +template +void multi_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states); template void dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const Matrix& matrix, StateVector& state); +template +void dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states); template void sparse_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const SparseMatrix& mat, StateVector& state); +template +void sparse_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const SparseMatrix& mat, + StateVectorBatched& states); template inline Matrix2x2 get_IBMQ_matrix(Fp _theta, Fp _phi, Fp _lambda); @@ -62,37 +102,72 @@ void one_target_diagonal_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const DiagonalMatrix2x2& diag, StateVector& state); +template +void one_target_diagonal_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const DiagonalMatrix2x2& diag, + StateVectorBatched& states); template inline void i_gate(std::uint64_t, std::uint64_t, StateVector&) {} +template +inline void i_gate(std::uint64_t, std::uint64_t, StateVectorBatched&) {} template void global_phase_gate(std::uint64_t, std::uint64_t control_mask, Fp angle, StateVector& state); +template +void global_phase_gate(std::uint64_t, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states); template void x_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state); +template +void x_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVectorBatched& states); template void y_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state); +template +void y_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVectorBatched& states); template void z_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state); +template +void z_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVectorBatched& states); template inline void h_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { one_target_dense_matrix_gate(target_mask, control_mask, HADAMARD_MATRIX(), state); } +template +inline void h_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, control_mask, HADAMARD_MATRIX(), states); +} template void one_target_phase_gate(std::uint64_t target_mask, std::uint64_t control_mask, Complex phase, StateVector& state); +template +void one_target_phase_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Complex phase, + StateVectorBatched& states); template inline void s_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { one_target_phase_gate(target_mask, control_mask, Complex(0, 1), state); } +template +inline void s_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_phase_gate(target_mask, control_mask, Complex(0, 1), states); +} template inline void sdag_gate(std::uint64_t target_mask, @@ -100,12 +175,25 @@ inline void sdag_gate(std::uint64_t target_mask, StateVector& state) { one_target_phase_gate(target_mask, control_mask, Complex(0, -1), state); } +template +inline void sdag_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_phase_gate(target_mask, control_mask, Complex(0, -1), states); +} template inline void t_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { one_target_phase_gate( target_mask, control_mask, Complex(INVERSE_SQRT2(), INVERSE_SQRT2()), state); } +template +inline void t_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_phase_gate( + target_mask, control_mask, Complex(INVERSE_SQRT2(), INVERSE_SQRT2()), states); +} template inline void tdag_gate(std::uint64_t target_mask, @@ -114,6 +202,13 @@ inline void tdag_gate(std::uint64_t target_mask, one_target_phase_gate( target_mask, control_mask, Complex(INVERSE_SQRT2(), -INVERSE_SQRT2()), state); } +template +inline void tdag_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_phase_gate( + target_mask, control_mask, Complex(INVERSE_SQRT2(), -INVERSE_SQRT2()), states); +} template inline void sqrtx_gate(std::uint64_t target_mask, @@ -121,6 +216,12 @@ inline void sqrtx_gate(std::uint64_t target_mask, StateVector& state) { one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_GATE_MATRIX(), state); } +template +inline void sqrtx_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_GATE_MATRIX(), states); +} template inline void sqrtxdag_gate(std::uint64_t target_mask, @@ -128,6 +229,12 @@ inline void sqrtxdag_gate(std::uint64_t target_mask, StateVector& state) { one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_DAG_GATE_MATRIX(), state); } +template +inline void sqrtxdag_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, control_mask, SQRT_X_DAG_GATE_MATRIX(), states); +} template inline void sqrty_gate(std::uint64_t target_mask, @@ -135,6 +242,12 @@ inline void sqrty_gate(std::uint64_t target_mask, StateVector& state) { one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_GATE_MATRIX(), state); } +template +inline void sqrty_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_GATE_MATRIX(), states); +} template inline void sqrtydag_gate(std::uint64_t target_mask, @@ -142,40 +255,96 @@ inline void sqrtydag_gate(std::uint64_t target_mask, StateVector& state) { one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_DAG_GATE_MATRIX(), state); } +template +inline void sqrtydag_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, control_mask, SQRT_Y_DAG_GATE_MATRIX(), states); +} template inline void p0_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { one_target_dense_matrix_gate(target_mask, control_mask, PROJ_0_MATRIX(), state); } +template +inline void p0_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, control_mask, PROJ_0_MATRIX(), states); +} template inline void p1_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { one_target_dense_matrix_gate(target_mask, control_mask, PROJ_1_MATRIX(), state); } +template +inline void p1_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, control_mask, PROJ_1_MATRIX(), states); +} template void rx_gate(std::uint64_t target_mask, std::uint64_t control_mask, Fp angle, StateVector& state); +template +void rx_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states); +template +void rx_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp pcoef, + std::vector params, + StateVectorBatched& states); template void ry_gate(std::uint64_t target_mask, std::uint64_t control_mask, Fp angle, StateVector& state); +template +void ry_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states); +template +void ry_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp pcoef, + std::vector params, + StateVectorBatched& states); template void rz_gate(std::uint64_t target_mask, std::uint64_t control_mask, Fp angle, StateVector& state); +template +void rz_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states); +template +void rz_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp pcoef, + std::vector params, + StateVectorBatched& states); template void u1_gate(std::uint64_t target_mask, std::uint64_t control_mask, Fp lambda, StateVector& state); +template +void u1_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp lambda, + StateVectorBatched& states); template void u2_gate(std::uint64_t target_mask, @@ -183,6 +352,12 @@ void u2_gate(std::uint64_t target_mask, Fp phi, Fp lambda, StateVector& state); +template +void u2_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp phi, + Fp lambda, + StateVectorBatched& states); template void u3_gate(std::uint64_t target_mask, @@ -191,20 +366,42 @@ void u3_gate(std::uint64_t target_mask, Fp phi, Fp lambda, StateVector& state); +template +void u3_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp theta, + Fp phi, + Fp lambda, + StateVectorBatched& states); template void swap_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state); +template +void swap_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states); template void sparse_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const SparseMatrix& matrix, StateVector& state); +template +void sparse_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const SparseMatrix& matrix, + StateVectorBatched& states); template void dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, const Matrix& matrix, StateVector& state); +template +void dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states); + } // namespace internal } // namespace scaluq diff --git a/src/gate/update_ops_dense_matrix.cpp b/src/gate/update_ops_dense_matrix.cpp index 77606b3d..f0731ce1 100644 --- a/src/gate/update_ops_dense_matrix.cpp +++ b/src/gate/update_ops_dense_matrix.cpp @@ -18,6 +18,160 @@ void none_target_dense_matrix_gate(std::uint64_t control_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void none_target_dense_matrix_gate(std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, {states.batch_size(), states.dim() >> std::popcount(control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis = insert_zero_at_mask_positions(it, control_mask) | control_mask; + states._raw(batch_id, basis) *= matrix(0, 0); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void none_target_dense_matrix_gate( \ + std::uint64_t, const Matrix&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + +template +void one_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix2x2& matrix, + StateVector& state) { + Kokkos::parallel_for( + state.dim() >> std::popcount(target_mask | control_mask), KOKKOS_LAMBDA(std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | target_mask; + Complex val0 = state._raw[basis_0]; + Complex val1 = state._raw[basis_1]; + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1; + state._raw[basis_0] = res0; + state._raw[basis_1] = res1; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void one_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix2x2&, StateVector&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + +template +void one_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix2x2& matrix, + StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void one_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix2x2&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + +template +void two_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix4x4& matrix, + StateVector& state) { + std::uint64_t lower_target_mask = -target_mask & target_mask; + std::uint64_t upper_target_mask = target_mask ^ lower_target_mask; + Kokkos::parallel_for( + state.dim() >> std::popcount(target_mask | control_mask), + KOKKOS_LAMBDA(const std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, target_mask | control_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | lower_target_mask; + std::uint64_t basis_2 = basis_0 | upper_target_mask; + std::uint64_t basis_3 = basis_1 | target_mask; + Complex val0 = state._raw[basis_0]; + Complex val1 = state._raw[basis_1]; + Complex val2 = state._raw[basis_2]; + Complex val3 = state._raw[basis_3]; + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1 + matrix[0][2] * val2 + + matrix[0][3] * val3; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1 + matrix[1][2] * val2 + + matrix[1][3] * val3; + Complex res2 = matrix[2][0] * val0 + matrix[2][1] * val1 + matrix[2][2] * val2 + + matrix[2][3] * val3; + Complex res3 = matrix[3][0] * val0 + matrix[3][1] * val1 + matrix[3][2] * val2 + + matrix[3][3] * val3; + state._raw[basis_0] = res0; + state._raw[basis_1] = res1; + state._raw[basis_2] = res2; + state._raw[basis_3] = res3; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void two_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix4x4&, StateVector&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + +template +void two_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix4x4& matrix, + StateVectorBatched& states) { + std::uint64_t lower_target_mask = -target_mask & target_mask; + std::uint64_t upper_target_mask = target_mask ^ lower_target_mask; + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, target_mask | control_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | lower_target_mask; + std::uint64_t basis_2 = basis_0 | upper_target_mask; + std::uint64_t basis_3 = basis_0 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex val2 = states._raw(batch_id, basis_2); + Complex val3 = states._raw(batch_id, basis_3); + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1 + matrix[0][2] * val2 + + matrix[0][3] * val3; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1 + matrix[1][2] * val2 + + matrix[1][3] * val3; + Complex res2 = matrix[2][0] * val0 + matrix[2][1] * val1 + matrix[2][2] * val2 + + matrix[2][3] * val3; + Complex res3 = matrix[3][0] * val0 + matrix[3][1] * val1 + matrix[3][2] * val2 + + matrix[3][3] * val3; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + states._raw(batch_id, basis_2) = res2; + states._raw(batch_id, basis_3) = res3; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void two_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix4x4&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void single_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -43,6 +197,34 @@ void single_target_dense_matrix_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void single_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex res0 = matrix(0, 0) * val0 + matrix(0, 1) * val1; + Complex res1 = matrix(1, 0) * val0 + matrix(1, 1) * val1; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void single_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void double_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -84,6 +266,49 @@ void double_target_dense_matrix_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void double_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states) { + std::uint64_t target_bit_right = -target_mask & target_mask; + std::uint64_t target_bit_left = target_mask ^ target_bit_right; + + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, target_mask | control_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | target_bit_right; + std::uint64_t basis_2 = basis_0 | target_bit_left; + std::uint64_t basis_3 = basis_0 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex val2 = states._raw(batch_id, basis_2); + Complex val3 = states._raw(batch_id, basis_3); + Complex res0 = matrix(0, 0) * val0 + matrix(0, 1) * val1 + matrix(0, 2) * val2 + + matrix(0, 3) * val3; + Complex res1 = matrix(1, 0) * val0 + matrix(1, 1) * val1 + matrix(1, 2) * val2 + + matrix(1, 3) * val3; + Complex res2 = matrix(2, 0) * val0 + matrix(2, 1) * val1 + matrix(2, 2) * val2 + + matrix(2, 3) * val3; + Complex res3 = matrix(3, 0) * val0 + matrix(3, 1) * val1 + matrix(3, 2) * val2 + + matrix(3, 3) * val3; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + states._raw(batch_id, basis_2) = res2; + states._raw(batch_id, basis_3) = res3; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void double_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void multi_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -136,6 +361,59 @@ void multi_target_dense_matrix_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void multi_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states) { + const std::uint64_t matrix_dim = 1ULL << std::popcount(target_mask); + const std::uint64_t outer_mask = ~target_mask & ((1ULL << states.n_qubits()) - 1); + + Kokkos::View**, Kokkos::LayoutRight> update( + Kokkos::ViewAllocateWithoutInitializing("update"), states.batch_size(), states.dim()); + + Kokkos::parallel_for( + Kokkos::MDRangePolicy>({0, 0}, {states.batch_size(), states.dim()}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t i) { + if ((i | control_mask) == i) { + update(batch_id, i) = 0; + } else { + update(batch_id, i) = states._raw(batch_id, i); + } + }); + Kokkos::fence(); + + Kokkos::View**, Kokkos::LayoutRight, Kokkos::MemoryTraits> + update_atomic(update); + + const std::uint64_t outer_dim = states.dim() >> std::popcount(target_mask | control_mask); + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0, 0, 0}, {states.batch_size(), outer_dim, matrix_dim, matrix_dim}), + KOKKOS_LAMBDA(const std::uint64_t batch_id, + const std::uint64_t outer_idx, + const std::uint64_t r, + const std::uint64_t c) { + std::uint64_t basis = + insert_zero_at_mask_positions(outer_idx, target_mask | control_mask) | control_mask; + + std::uint64_t dst_index = + internal::insert_zero_at_mask_positions(r, outer_mask) | basis; + + std::uint64_t src_index = + internal::insert_zero_at_mask_positions(c, outer_mask) | basis; + update_atomic(batch_id, dst_index) += matrix(r, c) * states._raw(batch_id, src_index); + }); + Kokkos::fence(); + + states._raw = update; +} +#define FUNC_MACRO(Fp) \ + template void multi_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -157,4 +435,26 @@ void dense_matrix_gate(std::uint64_t target_mask, std::uint64_t, std::uint64_t, const Matrix&, StateVector&); CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO + +template +void dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix& matrix, + StateVectorBatched& states) { + const std::uint64_t target_qubit_index_count = std::popcount(target_mask); + if (target_qubit_index_count == 0) { + none_target_dense_matrix_gate(control_mask, matrix, states); + } else if (target_qubit_index_count == 1) { + single_target_dense_matrix_gate(target_mask, control_mask, matrix, states); + } else if (target_qubit_index_count == 2) { + double_target_dense_matrix_gate(target_mask, control_mask, matrix, states); + } else { + multi_target_dense_matrix_gate(target_mask, control_mask, matrix, states); + } +} +#define FUNC_MACRO(Fp) \ + template void dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO } // namespace scaluq::internal diff --git a/src/gate/update_ops_sparse_matrix.cpp b/src/gate/update_ops_sparse_matrix.cpp index 00537e55..88694704 100644 --- a/src/gate/update_ops_sparse_matrix.cpp +++ b/src/gate/update_ops_sparse_matrix.cpp @@ -41,9 +41,60 @@ void sparse_matrix_gate(std::uint64_t target_mask, Kokkos::fence(); state._raw = update; } + +template +void sparse_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const SparseMatrix& mat, + StateVectorBatched& states) { + auto values = mat._values; + const std::uint64_t outer_mask = ~target_mask & ((1ULL << states.n_qubits()) - 1); + + Kokkos::View**, Kokkos::LayoutRight> update( + Kokkos::ViewAllocateWithoutInitializing("update"), states.batch_size(), states.dim()); + + Kokkos::parallel_for( + Kokkos::MDRangePolicy>({0, 0}, {states.batch_size(), states.dim()}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t i) { + if ((i | control_mask) == i) { + update(batch_id, i) = 0; + } else { + update(batch_id, i) = states._raw(batch_id, i); + } + }); + Kokkos::fence(); + + Kokkos::View**, Kokkos::LayoutRight, Kokkos::MemoryTraits> + update_atomic(update); + Kokkos::parallel_for( + "COO_Update", + Kokkos::MDRangePolicy>( + {0, 0, 0}, + {static_cast(states.batch_size()), + static_cast(states.dim() >> std::popcount(target_mask | control_mask)), + static_cast(values.size())}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t outer, std::uint64_t inner) { + std::uint64_t basis = + internal::insert_zero_at_mask_positions(outer, target_mask | control_mask) | + control_mask; + auto [v, r, c] = values(inner); + uint32_t src_index = internal::insert_zero_at_mask_positions(c, outer_mask) | basis; + uint32_t dst_index = internal::insert_zero_at_mask_positions(r, outer_mask) | basis; + update_atomic(batch_id, dst_index) += v * states._raw(batch_id, src_index); + }); + Kokkos::fence(); + states._raw = update; +} + #define FUNC_MACRO(Fp) \ template void sparse_matrix_gate( \ std::uint64_t, std::uint64_t, const SparseMatrix&, StateVector&); CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO + +#define FUNC_MACRO(Fp) \ + template void sparse_matrix_gate( \ + std::uint64_t, std::uint64_t, const SparseMatrix&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO } // namespace scaluq::internal diff --git a/src/gate/update_ops_standard.cpp b/src/gate/update_ops_standard.cpp index 3ff95b04..dcc1e4b4 100644 --- a/src/gate/update_ops_standard.cpp +++ b/src/gate/update_ops_standard.cpp @@ -39,6 +39,34 @@ void one_target_dense_matrix_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void one_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix2x2& matrix, + StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void one_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix2x2&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void two_target_dense_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -79,6 +107,48 @@ void two_target_dense_matrix_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void two_target_dense_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const Matrix4x4& matrix, + StateVectorBatched& states) { + std::uint64_t lower_target_mask = -target_mask & target_mask; + std::uint64_t upper_target_mask = target_mask ^ lower_target_mask; + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, target_mask | control_mask) | control_mask; + std::uint64_t basis_1 = basis_0 | lower_target_mask; + std::uint64_t basis_2 = basis_0 | upper_target_mask; + std::uint64_t basis_3 = basis_1 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex val2 = states._raw(batch_id, basis_2); + Complex val3 = states._raw(batch_id, basis_3); + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1 + matrix[0][2] * val2 + + matrix[0][3] * val3; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1 + matrix[1][2] * val2 + + matrix[1][3] * val3; + Complex res2 = matrix[2][0] * val0 + matrix[2][1] * val1 + matrix[2][2] * val2 + + matrix[2][3] * val3; + Complex res3 = matrix[3][0] * val0 + matrix[3][1] * val1 + matrix[3][2] * val2 + + matrix[3][3] * val3; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + states._raw(batch_id, basis_2) = res2; + states._raw(batch_id, basis_3) = res3; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void two_target_dense_matrix_gate( \ + std::uint64_t, std::uint64_t, const Matrix4x4&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void one_target_diagonal_matrix_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -99,6 +169,29 @@ void one_target_diagonal_matrix_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void one_target_diagonal_matrix_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + const DiagonalMatrix2x2& diag, + StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis = + insert_zero_at_mask_positions(it, target_mask | control_mask) | control_mask; + states._raw(batch_id, basis) *= diag[0]; + states._raw(batch_id, basis | target_mask) *= diag[1]; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void one_target_diagonal_matrix_gate( \ + std::uint64_t, std::uint64_t, const DiagonalMatrix2x2&, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void global_phase_gate(std::uint64_t, std::uint64_t control_mask, @@ -116,6 +209,26 @@ void global_phase_gate(std::uint64_t, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void global_phase_gate(std::uint64_t, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states) { + Complex coef = Kokkos::polar(1., angle); + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, {states.batch_size(), states.dim() >> std::popcount(control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t i) { + states._raw(batch_id, insert_zero_at_mask_positions(i, control_mask) | control_mask) *= + coef; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void global_phase_gate(std::uint64_t, std::uint64_t, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void x_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { Kokkos::parallel_for( @@ -130,6 +243,24 @@ void x_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector +void x_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t i = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + Kokkos::Experimental::swap(states._raw(batch_id, i), + states._raw(batch_id, i | target_mask)); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) template void x_gate(std::uint64_t, std::uint64_t, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void y_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { Kokkos::parallel_for( @@ -146,6 +277,26 @@ void y_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector +void y_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t i = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + states._raw(batch_id, i) *= Complex(0, 1); + states._raw(batch_id, i | target_mask) *= Complex(0, -1); + Kokkos::Experimental::swap(states._raw(batch_id, i), + states._raw(batch_id, i | target_mask)); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) template void y_gate(std::uint64_t, std::uint64_t, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void z_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { Kokkos::parallel_for( @@ -160,6 +311,23 @@ void z_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector +void z_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t i = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + states._raw(batch_id, i | target_mask) *= Complex(-1, 0); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) template void z_gate(std::uint64_t, std::uint64_t, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void one_target_phase_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -179,6 +347,28 @@ void one_target_phase_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void one_target_phase_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Complex phase, + StateVectorBatched& states) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t i = + insert_zero_at_mask_positions(it, control_mask | target_mask) | control_mask; + states._raw(batch_id, i | target_mask) *= phase; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void one_target_phase_gate( \ + std::uint64_t, std::uint64_t, Complex, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void rx_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -193,6 +383,60 @@ void rx_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void rx_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states) { + const Fp cosval = std::cos(angle / 2.); + const Fp sinval = std::sin(angle / 2.); + Matrix2x2 matrix = {cosval, Complex(0, -sinval), Complex(0, -sinval), cosval}; + one_target_dense_matrix_gate(target_mask, control_mask, matrix, states); +} +#define FUNC_MACRO(Fp) \ + template void rx_gate(std::uint64_t, std::uint64_t, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + +template +void rx_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp pcoef, + std::vector params, + StateVectorBatched& states) { + auto team_policy = Kokkos::TeamPolicy<>(states.batch_size(), Kokkos::AUTO); + Kokkos::parallel_for( + team_policy, KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type& team_member) { + const std::uint64_t batch_id = team_member.league_rank(); + const Fp angle = params[batch_id] * pcoef; + const Fp cosval = std::cos(angle / 2.); + const Fp sinval = std::sin(angle / 2.); + Matrix2x2 matrix = { + cosval, Complex(0, -sinval), Complex(0, -sinval), cosval}; + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team_member, + states.dim() >> std::popcount(target_mask | control_mask)), + [&](std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, control_mask | target_mask) | + control_mask; + std::uint64_t basis_1 = basis_0 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + }); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void rx_gate( \ + std::uint64_t, std::uint64_t, Fp, std::vector, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void ry_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -207,6 +451,59 @@ void ry_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void ry_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states) { + const Fp cosval = std::cos(angle / 2.); + const Fp sinval = std::sin(angle / 2.); + Matrix2x2 matrix = {cosval, -sinval, sinval, cosval}; + one_target_dense_matrix_gate(target_mask, control_mask, matrix, states); +} +#define FUNC_MACRO(Fp) \ + template void ry_gate(std::uint64_t, std::uint64_t, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + +template +void ry_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp pcoef, + std::vector params, + StateVectorBatched& states) { + auto team_policy = Kokkos::TeamPolicy<>(states.batch_size(), Kokkos::AUTO); + Kokkos::parallel_for( + team_policy, KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type& team_member) { + const std::uint64_t batch_id = team_member.league_rank(); + const Fp angle = params[batch_id] * pcoef; + const Fp cosval = std::cos(angle / 2.); + const Fp sinval = std::sin(angle / 2.); + Matrix2x2 matrix = {cosval, -sinval, sinval, cosval}; + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team_member, + states.dim() >> std::popcount(target_mask | control_mask)), + [&](std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, control_mask | target_mask) | + control_mask; + std::uint64_t basis_1 = basis_0 | target_mask; + Complex val0 = states._raw(batch_id, basis_0); + Complex val1 = states._raw(batch_id, basis_1); + Complex res0 = matrix[0][0] * val0 + matrix[0][1] * val1; + Complex res1 = matrix[1][0] * val0 + matrix[1][1] * val1; + states._raw(batch_id, basis_0) = res0; + states._raw(batch_id, basis_1) = res1; + }); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void ry_gate( \ + std::uint64_t, std::uint64_t, Fp, std::vector, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void rz_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -221,6 +518,56 @@ void rz_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void rz_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp angle, + StateVectorBatched& states) { + const Fp cosval = std::cos(angle / 2.); + const Fp sinval = std::sin(angle / 2.); + DiagonalMatrix2x2 diag = {Complex(cosval, -sinval), Complex(cosval, sinval)}; + one_target_diagonal_matrix_gate(target_mask, control_mask, diag, states); +} +#define FUNC_MACRO(Fp) \ + template void rz_gate(std::uint64_t, std::uint64_t, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + +template +void rz_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp pcoef, + std::vector params, + StateVectorBatched& states) { + auto team_policy = Kokkos::TeamPolicy<>(states.batch_size(), Kokkos::AUTO); + Kokkos::parallel_for( + team_policy, KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type& team_member) { + const std::uint64_t batch_id = team_member.league_rank(); + const Fp angle = params[batch_id] * pcoef; + const Fp cosval = std::cos(angle / 2.); + const Fp sinval = std::sin(angle / 2.); + DiagonalMatrix2x2 diag = {Complex(cosval, -sinval), + Complex(cosval, sinval)}; + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team_member, + states.dim() >> std::popcount(target_mask | control_mask)), + [&](std::uint64_t it) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(it, control_mask | target_mask) | + control_mask; + std::uint64_t basis_1 = basis_0 | target_mask; + states._raw(batch_id, basis_0) *= diag[0]; + states._raw(batch_id, basis_1) *= diag[1]; + }); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void rz_gate( \ + std::uint64_t, std::uint64_t, Fp, std::vector, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void u1_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -241,6 +588,29 @@ void u1_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void u1_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp lambda, + StateVectorBatched& states) { + Complex exp_val = Kokkos::exp(Complex(0, lambda)); + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t i = + internal::insert_zero_at_mask_positions(it, target_mask | control_mask) | + control_mask; + states._raw(batch_id, i | target_mask) *= exp_val; + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void u1_gate(std::uint64_t, std::uint64_t, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void u2_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -257,6 +627,22 @@ void u2_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void u2_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp phi, + Fp lambda, + StateVectorBatched& states) { + one_target_dense_matrix_gate(target_mask, + control_mask, + get_IBMQ_matrix((Fp)Kokkos::numbers::pi / 2, phi, lambda), + states); +} +#define FUNC_MACRO(Fp) \ + template void u2_gate(std::uint64_t, std::uint64_t, Fp, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void u3_gate(std::uint64_t target_mask, std::uint64_t control_mask, @@ -272,6 +658,21 @@ void u3_gate(std::uint64_t target_mask, CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +template +void u3_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + Fp theta, + Fp phi, + Fp lambda, + StateVectorBatched& states) { + one_target_dense_matrix_gate( + target_mask, control_mask, get_IBMQ_matrix(theta, phi, lambda), states); +} +#define FUNC_MACRO(Fp) \ + template void u3_gate(std::uint64_t, std::uint64_t, Fp, Fp, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO + template void swap_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVector& state) { // '- target' is used for bit manipulation on unsigned type, not for its numerical meaning. @@ -289,4 +690,27 @@ void swap_gate(std::uint64_t target_mask, std::uint64_t control_mask, StateVecto #define FUNC_MACRO(Fp) template void swap_gate(std::uint64_t, std::uint64_t, StateVector&); CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO + +template +void swap_gate(std::uint64_t target_mask, + std::uint64_t control_mask, + StateVectorBatched& states) { + std::uint64_t lower_target_mask = target_mask & -target_mask; + std::uint64_t upper_target_mask = target_mask ^ lower_target_mask; + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, + {states.batch_size(), states.dim() >> std::popcount(target_mask | control_mask)}), + KOKKOS_LAMBDA(std::uint64_t batch_id, std::uint64_t it) { + std::uint64_t basis = + insert_zero_at_mask_positions(it, target_mask | control_mask) | control_mask; + Kokkos::Experimental::swap(states._raw(batch_id, basis | lower_target_mask), + states._raw(batch_id, basis | upper_target_mask)); + }); + Kokkos::fence(); +} +#define FUNC_MACRO(Fp) \ + template void swap_gate(std::uint64_t, std::uint64_t, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO } // namespace scaluq::internal diff --git a/src/operator/apply_pauli.cpp b/src/operator/apply_pauli.cpp index 31bd4e46..e11fd2af 100644 --- a/src/operator/apply_pauli.cpp +++ b/src/operator/apply_pauli.cpp @@ -42,11 +42,58 @@ void apply_pauli(std::uint64_t control_mask, }); Kokkos::fence(); } +template +void apply_pauli(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + StateVectorBatched& states) { + if (bit_flip_mask == 0) { + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, {states.batch_size(), states.dim() >> std::popcount(control_mask)}), + KOKKOS_LAMBDA(const std::uint64_t batch_id, const std::uint64_t i) { + std::uint64_t state_idx = + insert_zero_at_mask_positions(i, control_mask) | control_mask; + if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) { + states._raw(batch_id, state_idx) *= -coef; + } else { + states._raw(batch_id, state_idx) *= coef; + } + }); + Kokkos::fence(); + return; + } + std::uint64_t pivot = sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; + std::uint64_t global_phase_90rot_count = std::popcount(bit_flip_mask & phase_flip_mask); + Complex global_phase = PHASE_M90ROT()[global_phase_90rot_count % 4]; + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, {states.batch_size(), states.dim() >> (std::popcount(control_mask) + 1)}), + KOKKOS_LAMBDA(const std::uint64_t batch_id, const std::uint64_t i) { + std::uint64_t basis_0 = + insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) | control_mask; + std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; + Complex tmp1 = states._raw(batch_id, basis_0) * global_phase; + Complex tmp2 = states._raw(batch_id, basis_1) * global_phase; + if (Kokkos::popcount(basis_0 & phase_flip_mask) & 1) tmp2 = -tmp2; + if (Kokkos::popcount(basis_1 & phase_flip_mask) & 1) tmp1 = -tmp1; + states._raw(batch_id, basis_0) = tmp2 * coef; + states._raw(batch_id, basis_1) = tmp1 * coef; + }); + Kokkos::fence(); +} + #define FUNC_MACRO(Fp) \ template void apply_pauli( \ std::uint64_t, std::uint64_t, std::uint64_t, Complex, StateVector&); CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +#define FUNC_MACRO(Fp) \ + template void apply_pauli( \ + std::uint64_t, std::uint64_t, std::uint64_t, Complex, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO template void apply_pauli_rotation(std::uint64_t control_mask, @@ -104,9 +151,144 @@ void apply_pauli_rotation(std::uint64_t control_mask, Kokkos::fence(); } } +template +void apply_pauli_rotation(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + Fp angle, + StateVectorBatched& states) { + std::uint64_t global_phase_90_rot_count = std::popcount(bit_flip_mask & phase_flip_mask); + Complex true_angle = angle * coef; + const Complex cosval = Kokkos::cos(-true_angle / 2); + const Complex sinval = Kokkos::sin(-true_angle / 2); + if (bit_flip_mask == 0) { + const Complex cval_min = cosval - Complex(0, 1) * sinval; + const Complex cval_pls = cosval + Complex(0, 1) * sinval; + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, {states.batch_size(), states.dim() >> std::popcount(control_mask)}), + KOKKOS_LAMBDA(const std::uint64_t batch_id, const std::uint64_t i) { + std::uint64_t state_idx = + insert_zero_at_mask_positions(i, control_mask) | control_mask; + if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) { + states._raw(batch_id, state_idx) *= cval_min; + } else { + states._raw(batch_id, state_idx) *= cval_pls; + } + }); + Kokkos::fence(); + return; + } else { + std::uint64_t pivot = sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; + Kokkos::parallel_for( + Kokkos::MDRangePolicy>( + {0, 0}, {states.batch_size(), states.dim() >> (std::popcount(control_mask) + 1)}), + KOKKOS_LAMBDA(const std::uint64_t batch_id, const std::uint64_t i) { + std::uint64_t basis_0 = + internal::insert_zero_at_mask_positions(i, control_mask | 1ULL << pivot) | + control_mask; + std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; + + int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) & 1; + int bit_parity_1 = Kokkos::popcount(basis_1 & phase_flip_mask) & 1; + + // fetch values + Complex cval_0 = states._raw(batch_id, basis_0); + Complex cval_1 = states._raw(batch_id, basis_1); + + // set values + states._raw(batch_id, basis_0) = + cosval * cval_0 + + Complex(0, 1) * sinval * cval_1 * + PHASE_M90ROT()[(global_phase_90_rot_count + bit_parity_0 * 2) % 4]; + states._raw(batch_id, basis_1) = + cosval * cval_1 + + Complex(0, 1) * sinval * cval_0 * + PHASE_M90ROT()[(global_phase_90_rot_count + bit_parity_1 * 2) % 4]; + }); + Kokkos::fence(); + } +} +template +void apply_pauli_rotation(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + Fp pcoef, + std::vector params, + StateVectorBatched& states) { + std::uint64_t global_phase_90_rot_count = std::popcount(bit_flip_mask & phase_flip_mask); + auto team_policy = Kokkos::TeamPolicy(states.batch_size(), Kokkos::AUTO); + Kokkos::parallel_for( + team_policy, KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type& team) { + const std::uint64_t batch_id = team.league_rank(); + const Fp angle = pcoef * params[batch_id]; + const Complex true_angle = angle * coef; + const Complex cosval = Kokkos::cos(-true_angle / 2); + const Complex sinval = Kokkos::sin(-true_angle / 2); + if (bit_flip_mask == 0) { + const Complex cval_min = cosval - Complex(0, 1) * sinval; + const Complex cval_pls = cosval + Complex(0, 1) * sinval; + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team, states.dim() >> std::popcount(control_mask)), + [&](const std::uint64_t i) { + std::uint64_t state_idx = + insert_zero_at_mask_positions(i, control_mask) | control_mask; + if (Kokkos::popcount(state_idx & phase_flip_mask) & 1) { + states._raw(batch_id, state_idx) *= cval_min; + } else { + states._raw(batch_id, state_idx) *= cval_pls; + } + }); + } else { + std::uint64_t pivot = + sizeof(std::uint64_t) * 8 - std::countl_zero(bit_flip_mask) - 1; + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team, + states.dim() >> (std::popcount(control_mask) + 1)), + [&](const std::uint64_t i) { + std::uint64_t basis_0 = internal::insert_zero_at_mask_positions( + i, control_mask | 1ULL << pivot) | + control_mask; + std::uint64_t basis_1 = basis_0 ^ bit_flip_mask; + int bit_parity_0 = Kokkos::popcount(basis_0 & phase_flip_mask) & 1; + int bit_parity_1 = Kokkos::popcount(basis_1 & phase_flip_mask) & 1; + Complex cval_0 = states._raw(batch_id, basis_0); + Complex cval_1 = states._raw(batch_id, basis_1); + states._raw(batch_id, basis_0) = + cosval * cval_0 + + Complex(0, 1) * sinval * cval_1 * + PHASE_M90ROT()[(global_phase_90_rot_count + bit_parity_0 * 2) % + 4]; + states._raw(batch_id, basis_1) = + cosval * cval_1 + + Complex(0, 1) * sinval * cval_0 * + PHASE_M90ROT()[(global_phase_90_rot_count + bit_parity_1 * 2) % + 4]; + }); + } + }); + Kokkos::fence(); +} #define FUNC_MACRO(Fp) \ template void apply_pauli_rotation( \ std::uint64_t, std::uint64_t, std::uint64_t, Complex, Fp, StateVector&); CALL_MACRO_FOR_FLOAT(FUNC_MACRO) #undef FUNC_MACRO +#define FUNC_MACRO(Fp) \ + template void apply_pauli_rotation( \ + std::uint64_t, std::uint64_t, std::uint64_t, Complex, Fp, StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO +#define FUNC_MACRO(Fp) \ + template void apply_pauli_rotation(std::uint64_t, \ + std::uint64_t, \ + std::uint64_t, \ + Complex, \ + Fp, \ + std::vector, \ + StateVectorBatched&); +CALL_MACRO_FOR_FLOAT(FUNC_MACRO) +#undef FUNC_MACRO } // namespace scaluq::internal diff --git a/src/operator/apply_pauli.hpp b/src/operator/apply_pauli.hpp index 9f50a953..572c2b56 100644 --- a/src/operator/apply_pauli.hpp +++ b/src/operator/apply_pauli.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include namespace scaluq::internal { @@ -10,7 +11,12 @@ void apply_pauli(std::uint64_t control_mask, std::uint64_t phase_flip_mask, Complex coef, StateVector& state_vector); - +template +void apply_pauli(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + StateVectorBatched& states); template void apply_pauli_rotation(std::uint64_t control_mask, std::uint64_t bit_flip_mask, @@ -18,4 +24,19 @@ void apply_pauli_rotation(std::uint64_t control_mask, Complex coef, Fp angle, StateVector& state_vector); +template +void apply_pauli_rotation(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + Fp angle, + StateVectorBatched& states); +template +void apply_pauli_rotation(std::uint64_t control_mask, + std::uint64_t bit_flip_mask, + std::uint64_t phase_flip_mask, + Complex coef, + Fp pcoef, + std::vector params, + StateVectorBatched& states); } // namespace scaluq::internal diff --git a/src/state/state_vector.cpp b/src/state/state_vector.cpp index 088bc300..b81f6969 100644 --- a/src/state/state_vector.cpp +++ b/src/state/state_vector.cpp @@ -11,6 +11,9 @@ StateVector::StateVector(std::uint64_t n_qubits) set_zero_state(); } FLOAT(Fp) +StateVector::StateVector(Kokkos::View view) + : _n_qubits(std::bit_width(view.extent(0)) - 1), _dim(view.extent(0)), _raw(view) {} +FLOAT(Fp) void StateVector::set_amplitude_at(std::uint64_t index, ComplexType c) { Kokkos::View host_view("single_value"); host_view() = c; @@ -157,7 +160,6 @@ Fp StateVector::get_entropy() const { ent); return ent; } - FLOAT(Fp) void StateVector::add_state_vector_with_coef(ComplexType coef, const StateVector& state) { Kokkos::parallel_for( diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 8665a542..b91dd1a6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,8 @@ add_executable(scaluq_test EXCLUDE_FROM_ALL circuit/circuit_test.cpp circuit/param_circuit_test.cpp gate/gate_test.cpp + gate/batched_gate_test.cpp + # gate/merge_test.cpp gate/merge_test.cpp gate/param_gate_test.cpp operator/test_pauli_operator.cpp diff --git a/tests/gate/batched_gate_test.cpp b/tests/gate/batched_gate_test.cpp new file mode 100644 index 00000000..4392237a --- /dev/null +++ b/tests/gate/batched_gate_test.cpp @@ -0,0 +1,1053 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../test_environment.hpp" +#include "../util/util.hpp" + +using namespace scaluq; + +const std::uint64_t BATCH_SIZE = 10; + +template (*QuantumGateConstructor)()> +void run_random_batched_gate_apply(std::uint64_t n_qubits) { + const int dim = 1ULL << n_qubits; + + ComplexVector test_state = ComplexVector::Zero(dim); + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (int i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + const Gate gate = QuantumGateConstructor(); + gate->update_quantum_state(states); + auto states_cp = states.get_amplitudes(); + + test_state = test_state; + + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (int i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template (*QuantumGateConstructor)(Fp, const std::vector&)> +void run_random_batched_gate_apply(std::uint64_t n_qubits) { + const int dim = 1ULL << n_qubits; + Random random; + + ComplexVector test_state = ComplexVector::Zero(dim); + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (int i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + const Fp angle = std::numbers::pi * random.uniform(); + const Gate gate = QuantumGateConstructor(angle, {}); + gate->update_quantum_state(states); + + test_state = std::polar(1, angle) * test_state; + + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + auto states_cp = states.get_amplitudes(); + for (int i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template (*QuantumGateConstructor)(std::uint64_t, const std::vector&)> +void run_random_batched_gate_apply(std::uint64_t n_qubits, + std::function()> matrix_factory) { + const auto matrix = matrix_factory(); + const int dim = 1ULL << n_qubits; + Random random; + + ComplexVector test_state = ComplexVector::Zero(dim); + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (int i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + const std::uint64_t target = random.int64() % n_qubits; + const Gate gate = QuantumGateConstructor(target, {}); + gate->update_quantum_state(states); + + test_state = + get_expanded_eigen_matrix_with_identity(target, matrix, n_qubits) * test_state; + + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (int i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template (*QuantumGateConstructor)(std::uint64_t, Fp, const std::vector&)> +void run_random_batched_gate_apply(std::uint64_t n_qubits, + std::function(Fp)> matrix_factory) { + const int dim = 1ULL << n_qubits; + Random random; + + ComplexVector test_state = ComplexVector::Zero(dim); + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (int i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + const Fp angle = std::numbers::pi * random.uniform(); + const auto matrix = matrix_factory(angle); + const std::uint64_t target = random.int64() % n_qubits; + const Gate gate = QuantumGateConstructor(target, angle, {}); + gate->update_quantum_state(states); + + test_state = + get_expanded_eigen_matrix_with_identity(target, matrix, n_qubits) * test_state; + + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (int i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template +void run_random_batched_gate_apply_IBMQ( + std::uint64_t n_qubits, std::function(Fp, Fp, Fp)> matrix_factory) { + const int dim = 1ULL << n_qubits; + Random random; + + ComplexVector test_state = ComplexVector::Zero(dim); + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + for (int gate_type = 0; gate_type < 3; gate_type++) { + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (int i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + Fp theta = std::numbers::pi * random.uniform(); + Fp phi = std::numbers::pi * random.uniform(); + Fp lambda = std::numbers::pi * random.uniform(); + if (gate_type == 0) { + theta = 0; + phi = 0; + } else if (gate_type == 1) { + theta = std::numbers::pi / 2; + } + const auto matrix = matrix_factory(theta, phi, lambda); + const std::uint64_t target = random.int64() % n_qubits; + Gate gate; + if (gate_type == 0) { + gate = gate::U1(target, lambda, {}); + } else if (gate_type == 1) { + gate = gate::U2(target, phi, lambda, {}); + } else { + gate = gate::U3(target, theta, phi, lambda, {}); + } + gate->update_quantum_state(states); + + test_state = + get_expanded_eigen_matrix_with_identity(target, matrix, n_qubits) * test_state; + + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (int i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } + } +} + +template +void run_random_batched_gate_apply_two_target(std::uint64_t n_qubits) { + const int dim = 1ULL << n_qubits; + Random random; + + ComplexVector test_state = ComplexVector::Zero(dim); + std::function(std::uint64_t, std::uint64_t, std::uint64_t)> func_eig; + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + for (int g = 0; g < 2; g++) { + Gate gate; + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (int i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + std::uint64_t target = random.int64() % n_qubits; + std::uint64_t control = random.int64() % n_qubits; + if (target == control) target = (target + 1) % n_qubits; + if (g == 0) { + gate = gate::CX(control, target); + func_eig = get_eigen_matrix_full_qubit_CX; + } else { + gate = gate::CZ(control, target); + func_eig = get_eigen_matrix_full_qubit_CZ; + } + gate->update_quantum_state(states); + + ComplexMatrix test_mat = func_eig(control, target, n_qubits); + test_state = test_mat * test_state; + + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (int i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } + } + + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (int i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + std::uint64_t target1 = random.int64() % n_qubits; + std::uint64_t target2 = random.int64() % n_qubits; + if (target1 == target2) target1 = (target1 + 1) % n_qubits; + auto gate = gate::Swap(target1, target2); + gate->update_quantum_state(states); + + ComplexMatrix test_mat = + get_eigen_matrix_full_qubit_Swap(target1, target2, n_qubits); + test_state = test_mat * test_state; + + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (int i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template +void run_random_batched_gate_apply_pauli(std::uint64_t n_qubits) { + const std::uint64_t dim = 1ULL << n_qubits; + Random random; + ComplexVector test_state = ComplexVector::Zero(dim); + ComplexMatrix matrix; + + // Test for PauliGate + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + auto states_bef = states.copy(); + + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + std::vector target_vec, pauli_id_vec; + for (std::uint64_t target = 0; target < n_qubits; target++) { + target_vec.emplace_back(target); + pauli_id_vec.emplace_back(random.int64() % 4); + } + + if (pauli_id_vec[0] == 0) { + matrix = make_I(); + } else if (pauli_id_vec[0] == 1) { + matrix = make_X(); + } else if (pauli_id_vec[0] == 2) { + matrix = make_Y(); + } else if (pauli_id_vec[0] == 3) { + matrix = make_Z(); + } + for (int i = 1; i < (int)n_qubits; i++) { + if (pauli_id_vec[i] == 0) { + matrix = internal::kronecker_product(make_I(), matrix); + } else if (pauli_id_vec[i] == 1) { + matrix = internal::kronecker_product(make_X(), matrix); + } else if (pauli_id_vec[i] == 2) { + matrix = internal::kronecker_product(make_Y(), matrix); + } else if (pauli_id_vec[i] == 3) { + matrix = internal::kronecker_product(make_Z(), matrix); + } + } + + PauliOperator pauli(target_vec, pauli_id_vec, 1.0); + Gate pauli_gate = gate::Pauli(pauli); + pauli_gate->update_quantum_state(states); + + test_state = matrix * test_state; + + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + + auto states_bef_cp = states_bef.get_amplitudes(); + Gate pauli_gate_inv = pauli_gate->get_inverse(); + pauli_gate_inv->update_quantum_state(states); + states_cp = states.get_amplitudes(); + + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], + (StdComplex)states_bef_cp[batch_id][i]); + } + } + } + + // Test for PauliRotationGate + for (int repeat = 0; repeat < 10; repeat++) { + auto states = StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + auto states_bef = states.copy(); + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + + const Fp angle = std::numbers::pi * random.uniform(); + std::vector target_vec, pauli_id_vec; + for (std::uint64_t target = 0; target < n_qubits; target++) { + target_vec.emplace_back(target); + pauli_id_vec.emplace_back(random.int64() % 4); + } + + if (pauli_id_vec[0] == 0) { + matrix = make_I(); + } else if (pauli_id_vec[0] == 1) { + matrix = make_X(); + } else if (pauli_id_vec[0] == 2) { + matrix = make_Y(); + } else if (pauli_id_vec[0] == 3) { + matrix = make_Z(); + } + for (int i = 1; i < (int)n_qubits; i++) { + if (pauli_id_vec[i] == 0) { + matrix = internal::kronecker_product(make_I(), matrix); + } else if (pauli_id_vec[i] == 1) { + matrix = internal::kronecker_product(make_X(), matrix); + } else if (pauli_id_vec[i] == 2) { + matrix = internal::kronecker_product(make_Y(), matrix); + } else if (pauli_id_vec[i] == 3) { + matrix = internal::kronecker_product(make_Z(), matrix); + } + } + matrix = std::cos(angle / 2) * ComplexMatrix::Identity(dim, dim) - + Complex(0, 1) * std::sin(angle / 2) * matrix; + + PauliOperator pauli(target_vec, pauli_id_vec, 1.0); + Gate pauli_gate = gate::PauliRotation(pauli, angle); + pauli_gate->update_quantum_state(states); + + test_state = matrix * test_state; + + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + + Gate pauli_gate_inv = pauli_gate->get_inverse(); + pauli_gate_inv->update_quantum_state(states); + states_cp = states.get_amplitudes(); + auto states_bef_cp = states_bef.get_amplitudes(); + + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], + (StdComplex)states_bef_cp[batch_id][i]); + } + } + } +} + +template +void run_random_batched_gate_apply_none_dense(std::uint64_t n_qubits) { + const std::uint64_t dim = 1ULL << n_qubits; + const std::uint64_t max_repeat = 10; + Eigen::Matrix, 1, 1, Eigen::RowMajor> U; + Random random; + ComplexVector test_state = ComplexVector::Zero(dim); + for (std::uint64_t rep = 0; rep < max_repeat; rep++) { + StateVectorBatched states = + StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + std::vector target_list = {}; + auto re = random.uniform(); + auto im = random.uniform(); + auto val = StdComplex(re, im); + auto norm = std::sqrt(std::norm(val)); + U(0, 0) = val / norm; + ComplexMatrix mat = ComplexMatrix::Identity(dim, dim); + mat *= val / norm; + std::vector control_list = {}; + Gate dense_gate = gate::DenseMatrix(target_list, U, control_list); + dense_gate->update_quantum_state(states); + test_state = mat * test_state; + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template +void run_random_batched_gate_apply_single_dense(std::uint64_t n_qubits) { + const std::uint64_t dim = 1ULL << n_qubits; + const std::uint64_t max_repeat = 10; + + Eigen::Matrix, 2, 2, Eigen::RowMajor> U; + std::uint64_t target; + Kokkos::View**> mat_view("mat_view", 2, 2); + Random random; + ComplexVector test_state = ComplexVector::Zero(dim); + for (std::uint64_t rep = 0; rep < max_repeat; rep++) { + StateVectorBatched states = + StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + target = random.int64() % n_qubits; + std::vector target_list = {target}; + U = get_eigen_matrix_random_one_target_unitary(); + ComplexMatrix mat(U.rows(), U.cols()); + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 2; j++) { + mat(i, j) = U(i, j); + } + } + std::vector control_list = {}; + Gate dense_gate = gate::DenseMatrix(target_list, mat, control_list); + dense_gate->update_quantum_state(states); + test_state = get_expanded_eigen_matrix_with_identity(target, U, n_qubits) * test_state; + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template +void run_random_batched_gate_apply_sparse(std::uint64_t n_qubits) { + const std::uint64_t dim = 1ULL << n_qubits; + const std::uint64_t max_repeat = 10; + + ComplexVector test_state = ComplexVector::Zero(dim); + internal::SparseComplexMatrix mat; + std::vector targets(3); + std::vector index_list; + std::random_device seed_gen; + std::mt19937 engine(seed_gen()); + Eigen::Matrix, 2, 2, Eigen::RowMajor> u1, u2, u3; + Eigen::Matrix, 8, 8, Eigen::RowMajor> Umerge; + for (std::uint64_t i = 0; i < n_qubits; i++) { + index_list.push_back(i); + } + for (std::uint64_t rep = 0; rep < max_repeat; rep++) { + StateVectorBatched states = + StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + std::shuffle(index_list.begin(), index_list.end(), engine); + targets[0] = index_list[0]; + targets[1] = index_list[1]; + targets[2] = index_list[2]; + u1 = get_eigen_matrix_random_one_target_unitary(); + u2 = get_eigen_matrix_random_one_target_unitary(); + u3 = get_eigen_matrix_random_one_target_unitary(); + std::vector target_list = {targets[0], targets[1], targets[2]}; + std::vector control_list = {}; + + test_state = get_expanded_eigen_matrix_with_identity(target_list[2], u3, n_qubits) * + get_expanded_eigen_matrix_with_identity(target_list[1], u2, n_qubits) * + get_expanded_eigen_matrix_with_identity(target_list[0], u1, n_qubits) * + test_state; + + Umerge = internal::kronecker_product(u3, internal::kronecker_product(u2, u1)); + mat = Umerge.sparseView(); + Gate sparse_gate = gate::SparseMatrix(target_list, mat, control_list); + sparse_gate->update_quantum_state(states); + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } +} + +template +void run_random_batched_gate_apply_general_dense(std::uint64_t n_qubits) { + const std::uint64_t dim = 1ULL << n_qubits; + const std::uint64_t max_repeat = 10; + + ComplexVector test_state = ComplexVector::Zero(dim); + Eigen::Matrix, 2, 2, Eigen::RowMajor> U1, U2, U3; + std::vector targets(3); + std::vector index_list; + std::random_device seed_gen; + std::mt19937 engine(seed_gen()); + for (std::uint64_t i = 0; i < n_qubits; i++) { + index_list.push_back(i); + } + // general single + { + for (std::uint64_t rep = 0; rep < max_repeat; rep++) { + StateVectorBatched states = + StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + U1 = get_eigen_matrix_random_one_target_unitary(); + std::shuffle(index_list.begin(), index_list.end(), engine); + targets[0] = index_list[0]; + test_state = + get_expanded_eigen_matrix_with_identity(targets[0], U1, n_qubits) * test_state; + std::vector target_list = {targets[0]}; + std::vector control_list = {}; + Gate dense_gate = gate::DenseMatrix(target_list, U1, control_list); + dense_gate->update_quantum_state(states); + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } + } + // general double + { + Eigen::Matrix, 4, 4, Eigen::RowMajor> Umerge; + for (std::uint64_t rep = 0; rep < max_repeat; rep++) { + StateVectorBatched states = + StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + U1 = get_eigen_matrix_random_one_target_unitary(); + U2 = get_eigen_matrix_random_one_target_unitary(); + + std::shuffle(index_list.begin(), index_list.end(), engine); + targets[0] = index_list[0]; + targets[1] = index_list[1]; + Umerge = internal::kronecker_product(U2, U1); + test_state = get_expanded_eigen_matrix_with_identity(targets[1], U2, n_qubits) * + get_expanded_eigen_matrix_with_identity(targets[0], U1, n_qubits) * + test_state; + std::vector target_list = {targets[0], targets[1]}; + std::vector control_list = {}; + Gate dense_gate = gate::DenseMatrix(target_list, Umerge, control_list); + dense_gate->update_quantum_state(states); + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } + } + // general triple + { + Eigen::Matrix, 8, 8, Eigen::RowMajor> Umerge; + for (std::uint64_t rep = 0; rep < max_repeat; rep++) { + StateVectorBatched states = + StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto state_cp = states.get_state_vector_at(0).get_amplitudes(); + for (std::uint64_t i = 0; i < dim; i++) { + test_state[i] = state_cp[i]; + } + U1 = get_eigen_matrix_random_one_target_unitary(); + U2 = get_eigen_matrix_random_one_target_unitary(); + U3 = get_eigen_matrix_random_one_target_unitary(); + + std::shuffle(index_list.begin(), index_list.end(), engine); + targets[0] = index_list[0]; + targets[1] = index_list[1]; + targets[2] = index_list[2]; + Umerge = internal::kronecker_product(U3, internal::kronecker_product(U2, U1)); + + test_state = get_expanded_eigen_matrix_with_identity(targets[2], U3, n_qubits) * + get_expanded_eigen_matrix_with_identity(targets[1], U2, n_qubits) * + get_expanded_eigen_matrix_with_identity(targets[0], U1, n_qubits) * + test_state; + std::vector target_list = {targets[0], targets[1], targets[2]}; + std::vector control_list = {}; + Gate dense_gate = gate::DenseMatrix(target_list, Umerge, control_list); + dense_gate->update_quantum_state(states); + auto states_cp = states.get_amplitudes(); + for (std::uint64_t batch_id = 0; batch_id < states.batch_size(); batch_id++) { + for (std::uint64_t i = 0; i < dim; i++) { + check_near((StdComplex)states_cp[batch_id][i], test_state[i]); + } + } + } + } +} + +TEST(BatchedGateTest, ApplyI) { + run_random_batched_gate_apply>(5); + run_random_batched_gate_apply>(5); +} +TEST(BatchedGateTest, ApplyGlobalPhase) { + run_random_batched_gate_apply>(5); + run_random_batched_gate_apply>(5); +} +TEST(BatchedGateTest, ApplyX) { + run_random_batched_gate_apply>(5, make_X); + run_random_batched_gate_apply>(5, make_X); +} +TEST(BatchedGateTest, ApplyY) { + run_random_batched_gate_apply>(5, make_Y); + run_random_batched_gate_apply>(5, make_Y); +} +TEST(BatchedGateTest, ApplyZ) { + run_random_batched_gate_apply>(5, make_Z); + run_random_batched_gate_apply>(5, make_Z); +} +TEST(BatchedGateTest, ApplyH) { + run_random_batched_gate_apply>(5, make_H); + run_random_batched_gate_apply>(5, make_H); +} +TEST(BatchedGateTest, ApplyS) { + run_random_batched_gate_apply>(5, make_S); + run_random_batched_gate_apply>(5, make_S); +} +TEST(BatchedGateTest, ApplySdag) { + run_random_batched_gate_apply>(5, make_Sdag); + run_random_batched_gate_apply>(5, make_Sdag); +} +TEST(BatchedGateTest, ApplyT) { + run_random_batched_gate_apply>(5, make_T); + run_random_batched_gate_apply>(5, make_T); +} +TEST(BatchedGateTest, ApplyTdag) { + run_random_batched_gate_apply>(5, make_Tdag); + run_random_batched_gate_apply>(5, make_Tdag); +} +TEST(BatchedGateTest, ApplySqrtX) { + run_random_batched_gate_apply>(5, make_SqrtX); + run_random_batched_gate_apply>(5, make_SqrtX); +} +TEST(BatchedGateTest, ApplySqrtY) { + run_random_batched_gate_apply>(5, make_SqrtY); + run_random_batched_gate_apply>(5, make_SqrtY); +} +TEST(BatchedGateTest, ApplySqrtXdag) { + run_random_batched_gate_apply>(5, make_SqrtXdag); + run_random_batched_gate_apply>(5, make_SqrtXdag); +} +TEST(BatchedGateTest, ApplySqrtYdag) { + run_random_batched_gate_apply>(5, make_SqrtYdag); + run_random_batched_gate_apply>(5, make_SqrtYdag); +} +TEST(BatchedGateTest, ApplyP0) { + run_random_batched_gate_apply>(5, make_P0); + run_random_batched_gate_apply>(5, make_P0); +} +TEST(BatchedGateTest, ApplyP1) { + run_random_batched_gate_apply>(5, make_P1); + run_random_batched_gate_apply>(5, make_P1); +} +TEST(BatchedGateTest, ApplyRX) { + run_random_batched_gate_apply>(5, make_RX); + run_random_batched_gate_apply>(5, make_RX); +} +TEST(BatchedGateTest, ApplyRY) { + run_random_batched_gate_apply>(5, make_RY); + run_random_batched_gate_apply>(5, make_RY); +} +TEST(BatchedGateTest, ApplyRZ) { + run_random_batched_gate_apply>(5, make_RZ); + run_random_batched_gate_apply>(5, make_RZ); +} + +TEST(BatchedGateTest, ApplyIBMQ) { + run_random_batched_gate_apply_IBMQ(5, make_U); + run_random_batched_gate_apply_IBMQ(5, make_U); +} +TEST(BatchedGateTest, ApplySparseMatrixGate) { + run_random_batched_gate_apply_sparse(6); + run_random_batched_gate_apply_sparse(6); +} +TEST(BatchedGateTest, ApplyDenseMatrixGate) { + run_random_batched_gate_apply_none_dense(6); + run_random_batched_gate_apply_none_dense(6); + run_random_batched_gate_apply_single_dense(6); + run_random_batched_gate_apply_single_dense(6); + run_random_batched_gate_apply_general_dense(6); + run_random_batched_gate_apply_general_dense(6); +} +TEST(BatchedGateTest, ApplyPauliGate) { + run_random_batched_gate_apply_pauli(5); + run_random_batched_gate_apply_pauli(5); +} + +TEST(BatchedGateTest, ApplyProbablisticGate) { + { + auto probgate = + gate::Probablistic({.1, .9}, {gate::X(0), gate::I()}); + StateVectorBatched states(BATCH_SIZE, 1); + std::vector> befores, afters; + std::vector x_counts(BATCH_SIZE), i_counts(BATCH_SIZE); + for ([[maybe_unused]] auto _ : std::views::iota(0, 100)) { + befores = states.sampling(1); + probgate->update_quantum_state(states); + afters = states.sampling(1); + for (std::size_t i = 0; i < BATCH_SIZE; i++) { + if (befores[i][0] != afters[i][0]) { + x_counts[i]++; + } else { + i_counts[i]++; + } + } + } + // These test is probablistic, but pass at least 99.9% cases. + for (std::size_t i = 0; i < BATCH_SIZE; i++) { + ASSERT_GT(x_counts[i], 0); + ASSERT_GT(i_counts[i], 0); + ASSERT_LT(x_counts[i], i_counts[i]); + } + } + { + auto probgate = gate::Probablistic({.1, .9}, {gate::X(0), gate::I()}); + StateVectorBatched states(BATCH_SIZE, 1); + std::vector> befores, afters; + std::vector x_counts(BATCH_SIZE), i_counts(BATCH_SIZE); + for ([[maybe_unused]] auto _ : std::views::iota(0, 100)) { + befores = states.sampling(1); + probgate->update_quantum_state(states); + afters = states.sampling(1); + for (std::size_t i = 0; i < BATCH_SIZE; i++) { + if (befores[i][0] != afters[i][0]) { + x_counts[i]++; + } else { + i_counts[i]++; + } + } + } + // These test is probablistic, but pass at least 99.9% cases. + for (std::size_t i = 0; i < BATCH_SIZE; i++) { + ASSERT_GT(x_counts[i], 0); + ASSERT_GT(i_counts[i], 0); + ASSERT_LT(x_counts[i], i_counts[i]); + } + } +} + +template +void test_batched_gate(Gate gate_control, + Gate gate_simple, + std::uint64_t n_qubits, + std::uint64_t control_mask) { + StateVectorBatched states = + StateVectorBatched::Haar_random_state(BATCH_SIZE, n_qubits, true); + auto amplitudes = states.get_amplitudes(); + StateVectorBatched states_controlled(BATCH_SIZE, n_qubits - std::popcount(control_mask)); + std::vector>> amplitudes_controlled( + BATCH_SIZE, std::vector>(states_controlled.dim())); + for (std::size_t i = 0; i < BATCH_SIZE; i++) { + for (std::uint64_t j = 0; j < states_controlled.dim(); j++) { + amplitudes_controlled[i][j] = + amplitudes[i] + [internal::insert_zero_at_mask_positions(j, control_mask) | control_mask]; + } + } + states_controlled.load(amplitudes_controlled); + gate_control->update_quantum_state(states); + gate_simple->update_quantum_state(states_controlled); + amplitudes = states.get_amplitudes(); + amplitudes_controlled = states_controlled.get_amplitudes(); + for (std::uint64_t i = 0; i < BATCH_SIZE; i++) { + for (std::uint64_t j : std::views::iota(0ULL, states_controlled.dim())) { + check_near((StdComplex)amplitudes_controlled[i][j], + (StdComplex) + amplitudes[i][internal::insert_zero_at_mask_positions(j, control_mask) | + control_mask]); + } + } +} + +template +void test_batched_standard_gate_control(Factory factory, std::uint64_t n) { + Random random; + std::vector shuffled = random.permutation(n); + std::vector targets(num_target); + for (std::uint64_t i : std::views::iota(0ULL, num_target)) { + targets[i] = shuffled[i]; + } + std::uint64_t num_control = random.int32() % (n - num_target + 1); + std::vector controls(num_control); + for (std::uint64_t i : std::views::iota(0ULL, num_control)) { + controls[i] = shuffled[num_target + i]; + } + std::uint64_t control_mask = 0ULL; + for (std::uint64_t c : controls) control_mask |= 1ULL << c; + std::vector angles(num_rotation); + for (Fp& angle : angles) angle = random.uniform() * std::numbers::pi * 2; + if constexpr (num_target == 0 && num_rotation == 1) { + Gate g1 = factory(angles[0], controls); + Gate g2 = factory(angles[0], {}); + test_batched_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 0) { + Gate g1 = factory(targets[0], controls); + Gate g2 = + factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), {}); + test_batched_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 1) { + Gate g1 = factory(targets[0], angles[0], controls); + Gate g2 = factory( + targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), angles[0], {}); + test_batched_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 2) { + Gate g1 = factory(targets[0], angles[0], angles[1], controls); + Gate g2 = factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), + angles[0], + angles[1], + {}); + test_batched_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 1 && num_rotation == 3) { + Gate g1 = factory(targets[0], angles[0], angles[1], angles[2], controls); + Gate g2 = factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), + angles[0], + angles[1], + angles[2], + {}); + test_batched_gate(g1, g2, n, control_mask); + } else if constexpr (num_target == 2 && num_rotation == 0) { + Gate g1 = factory(targets[0], targets[1], controls); + Gate g2 = factory(targets[0] - std::popcount(control_mask & ((1ULL << targets[0]) - 1)), + targets[1] - std::popcount(control_mask & ((1ULL << targets[1]) - 1)), + {}); + test_batched_gate(g1, g2, n, control_mask); + } else { + FAIL(); + } +} + +template +void test_batched_pauli_control(std::uint64_t n) { + typename PauliOperator::Data data1, data2; + std::vector controls; + std::uint64_t control_mask = 0; + std::uint64_t num_control = 0; + Random random; + for (std::uint64_t i : std::views::iota(0ULL, n)) { + std::uint64_t dat = random.int32() % 12; + if (dat < 4) { + data1.add_single_pauli(i, dat); + data2.add_single_pauli(i - num_control, dat); + } else if (dat < 8) { + controls.push_back(i); + control_mask |= 1ULL << i; + num_control++; + } + } + if constexpr (!rotation) { + Gate g1 = gate::Pauli(PauliOperator(data1), controls); + Gate g2 = gate::Pauli(PauliOperator(data2), {}); + test_batched_gate(g1, g2, n, control_mask); + } else { + Fp angle = random.uniform() * std::numbers::pi * 2; + Gate g1 = gate::PauliRotation(PauliOperator(data1), angle, controls); + Gate g2 = gate::PauliRotation(PauliOperator(data2), angle, {}); + test_batched_gate(g1, g2, n, control_mask); + } +} + +template +void test_batched_matrix_control(std::uint64_t n_qubits) { + Random random; + std::vector shuffled(n_qubits); + std::iota(shuffled.begin(), shuffled.end(), 0ULL); + for (std::uint64_t i : std::views::iota(0ULL, n_qubits) | std::views::reverse) { + std::uint64_t j = random.int32() % (i + 1); + if (i != j) std::swap(shuffled[i], shuffled[j]); + } + std::vector targets(num_target); + for (std::uint64_t i : std::views::iota(0ULL, num_target)) { + targets[i] = shuffled[i]; + } + std::uint64_t num_control = random.int32() % (n_qubits - num_target + 1); + std::vector controls(num_control); + for (std::uint64_t i : std::views::iota(0ULL, num_control)) { + controls[i] = shuffled[num_target + i]; + } + std::uint64_t control_mask = 0ULL; + for (std::uint64_t c : controls) control_mask |= 1ULL << c; + auto adjust = [](std::vector targets, std::uint64_t control_mask) { + std::vector new_targets; + for (auto i : targets) { + new_targets.push_back(i - std::popcount(control_mask & ((1ULL << i) - 1))); + } + return new_targets; + }; + auto new_targets = adjust(targets, control_mask); + if constexpr (num_target == 0) { + auto re = random.uniform(); + auto im = random.uniform(); + auto val = StdComplex(re, im); + auto norm = std::sqrt(std::norm(val)); + ComplexMatrix mat(1, 1); + mat(0, 0) = val / norm; + Gate d1 = gate::DenseMatrix(targets, mat, controls); + Gate d2 = gate::DenseMatrix(new_targets, mat, {}); + Gate s1 = gate::SparseMatrix(targets, mat.sparseView(), controls); + Gate s2 = gate::SparseMatrix(new_targets, mat.sparseView(), {}); + test_batched_gate(d1, d2, n_qubits, control_mask); + test_batched_gate(s1, s2, n_qubits, control_mask); + } else if constexpr (num_target == 1) { + Eigen::Matrix, 2, 2, Eigen::RowMajor> U = + get_eigen_matrix_random_one_target_unitary(); + ComplexMatrix mat(U.rows(), U.cols()); + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 2; j++) { + mat(i, j) = U(i, j); + } + } + Gate d1 = gate::DenseMatrix(targets, mat, controls); + Gate d2 = gate::DenseMatrix(new_targets, mat, {}); + Gate s1 = gate::SparseMatrix(targets, mat.sparseView(), controls); + Gate s2 = gate::SparseMatrix(new_targets, mat.sparseView(), {}); + test_batched_gate(d1, d2, n_qubits, control_mask); + test_batched_gate(s1, s2, n_qubits, control_mask); + } else if constexpr (num_target == 2) { + Eigen::Matrix, 2, 2, Eigen::RowMajor> U1 = + get_eigen_matrix_random_one_target_unitary(); + Eigen::Matrix, 2, 2, Eigen::RowMajor> U2 = + get_eigen_matrix_random_one_target_unitary(); + auto U = internal::kronecker_product(U2, U1); + ComplexMatrix mat(U.rows(), U.cols()); + for (int i = 0; i < 4; i++) { + for (int j = 0; j < 4; j++) { + mat(i, j) = U(i, j); + } + } + Gate d1 = gate::DenseMatrix(targets, mat, controls); + Gate d2 = gate::DenseMatrix(new_targets, mat, {}); + Gate s1 = gate::SparseMatrix(targets, mat.sparseView(), controls); + Gate s2 = gate::SparseMatrix(new_targets, mat.sparseView(), {}); + test_batched_gate(d1, d2, n_qubits, control_mask); + test_batched_gate(s1, s2, n_qubits, control_mask); + } else { + Eigen::Matrix, 2, 2, Eigen::RowMajor> U1 = + get_eigen_matrix_random_one_target_unitary(); + Eigen::Matrix, 2, 2, Eigen::RowMajor> U2 = + get_eigen_matrix_random_one_target_unitary(); + Eigen::Matrix, 2, 2, Eigen::RowMajor> U3 = + get_eigen_matrix_random_one_target_unitary(); + auto U = internal::kronecker_product(U3, internal::kronecker_product(U2, U1)); + ComplexMatrix mat(U.rows(), U.cols()); + for (int i = 0; i < 8; i++) { + for (int j = 0; j < 8; j++) { + mat(i, j) = U(i, j); + } + } + Gate d1 = gate::DenseMatrix(targets, mat, controls); + Gate d2 = gate::DenseMatrix(new_targets, mat, {}); + Gate s1 = gate::SparseMatrix(targets, mat.sparseView(), controls); + Gate s2 = gate::SparseMatrix(new_targets, mat.sparseView(), {}); + test_batched_gate(d1, d2, n_qubits, control_mask); + test_batched_gate(s1, s2, n_qubits, control_mask); + } +} +TEST(BatchGateTest, Control) { + std::uint64_t n = 10; + for ([[maybe_unused]] std::uint64_t _ : std::views::iota(0, 10)) { + test_batched_standard_gate_control(gate::GlobalPhase, n); + test_batched_standard_gate_control(gate::X, n); + test_batched_standard_gate_control(gate::Y, n); + test_batched_standard_gate_control(gate::Z, n); + test_batched_standard_gate_control(gate::S, n); + test_batched_standard_gate_control(gate::Sdag, n); + test_batched_standard_gate_control(gate::T, n); + test_batched_standard_gate_control(gate::Tdag, n); + test_batched_standard_gate_control(gate::SqrtX, n); + test_batched_standard_gate_control(gate::SqrtXdag, n); + test_batched_standard_gate_control(gate::SqrtY, n); + test_batched_standard_gate_control(gate::SqrtYdag, n); + test_batched_standard_gate_control(gate::P0, n); + test_batched_standard_gate_control(gate::P1, n); + test_batched_standard_gate_control(gate::RX, n); + test_batched_standard_gate_control(gate::RY, n); + test_batched_standard_gate_control(gate::RZ, n); + test_batched_standard_gate_control(gate::U1, n); + test_batched_standard_gate_control(gate::U2, n); + test_batched_standard_gate_control(gate::U3, n); + test_batched_standard_gate_control(gate::Swap, n); + test_batched_pauli_control(n); + test_batched_pauli_control(n); + test_batched_matrix_control(n); + test_batched_matrix_control(n); + test_batched_matrix_control(n); + test_batched_matrix_control(n); + + test_batched_standard_gate_control(gate::GlobalPhase, n); + test_batched_standard_gate_control(gate::X, n); + test_batched_standard_gate_control(gate::Y, n); + test_batched_standard_gate_control(gate::Z, n); + test_batched_standard_gate_control(gate::S, n); + test_batched_standard_gate_control(gate::Sdag, n); + test_batched_standard_gate_control(gate::T, n); + test_batched_standard_gate_control(gate::Tdag, n); + test_batched_standard_gate_control(gate::SqrtX, n); + test_batched_standard_gate_control(gate::SqrtXdag, n); + test_batched_standard_gate_control(gate::SqrtY, n); + test_batched_standard_gate_control(gate::SqrtYdag, n); + test_batched_standard_gate_control(gate::P0, n); + test_batched_standard_gate_control(gate::P1, n); + test_batched_standard_gate_control(gate::RX, n); + test_batched_standard_gate_control(gate::RY, n); + test_batched_standard_gate_control(gate::RZ, n); + test_batched_standard_gate_control(gate::U1, n); + test_batched_standard_gate_control(gate::U2, n); + test_batched_standard_gate_control(gate::U3, n); + test_batched_standard_gate_control(gate::Swap, n); + test_batched_pauli_control(n); + test_batched_pauli_control(n); + test_batched_matrix_control(n); + test_batched_matrix_control(n); + test_batched_matrix_control(n); + test_batched_matrix_control(n); + } +} diff --git a/tests/gate/gate_test.cpp b/tests/gate/gate_test.cpp index d60755ec..ae98378f 100644 --- a/tests/gate/gate_test.cpp +++ b/tests/gate/gate_test.cpp @@ -656,11 +656,6 @@ TEST(GateTest, ApplyIBMQ) { run_random_gate_apply_IBMQ(5, make_U); } -TEST(GateTest, ApplyTwoTarget) { - run_random_gate_apply_two_target(5); - run_random_gate_apply_two_target(5); -} - TEST(GateTest, ApplySparseMatrixGate) { run_random_gate_apply_sparse(6); run_random_gate_apply_sparse(6);