Merge pull request #34 from Qulacs-Osaka/29-add-swap

29 add swap

qulacs/CMakeLists.txt

@@ -1,13 +1,17 @@
 cmake_minimum_required(VERSION 3.21)
 
 target_sources(qulacs PRIVATE
+    gate/gate_npair_qubit.cpp
     gate/gate_one_control_one_target.cpp
     gate/gate_one_qubit.cpp
     gate/gate_quantum_matrix.cpp
-    gate/update_ops_matrix_dense.cpp
+    gate/gate_two_qubit.cpp
+    gate/update_ops_dense_matrix.cpp
+    gate/update_ops_npair_qubit.cpp
+    gate/update_ops_one_control_one_target.cpp
     gate/update_ops_one_qubit.cpp
-    gate/update_ops_one_target_one_control.cpp
     gate/update_ops_quantum_matrix.cpp
+    gate/update_ops_two_qubit.cpp
     state/state_vector.cpp
     util/random.cpp
 )

qulacs/all.hpp

@@ -0,0 +1,11 @@
+#include "gate/constant.hpp"
+#include "gate/gate.hpp"
+#include "gate/gate_npair_qubit.hpp"
+#include "gate/gate_one_control_one_target.hpp"
+#include "gate/gate_one_qubit.hpp"
+#include "gate/gate_quantum_matrix.hpp"
+#include "gate/gate_two_qubit.hpp"
+#include "gate/update_ops.hpp"
+#include "state/state_vector.hpp"
+#include "types.hpp"
+#include "util/random.hpp"

qulacs/gate/constant.hpp

@@ -27,47 +27,72 @@ KOKKOS_INLINE_FUNCTION
 double SINPI8() { return 0.382683432365090; }
 
 //! identity matrix
-matrix_2_2 I_GATE = {1, 0, 0, 1};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 I_GATE() { return {1, 0, 0, 1}; }
 //! Pauli matrix X
-matrix_2_2 X_GATE = {0, 1, 1, 0};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 X_GATE() { return {0, 1, 1, 0}; }
 //! Pauli matrix Y
-matrix_2_2 Y_GATE = {0, Complex(0, -1), Complex(0, 1), 0};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 Y_GATE() { return {0, Complex(0, -1), Complex(0, 1), 0}; }
 //! Pauli matrix Z
-matrix_2_2 Z_GATE = {1, 0, 0, -1};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 Z_GATE() { return {1, 0, 0, -1}; }
 
 //! list of Pauli matrix I,X,Y,Z
 // std::array<matrix_2_2, 4> PAULI_MATRIX = {I_GATE, X_GATE, Y_GATE, Z_GATE};
 
 //! S-gate
-matrix_2_2 S_GATE_MATRIX = {1, 0, 0, Complex(0, 1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 S_GATE_MATRIX() { return {1, 0, 0, Complex(0, 1)}; }
 //! Sdag-gate
-matrix_2_2 S_DAG_GATE_MATRIX = {1, 0, 0, Complex(0, -1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 S_DAG_GATE_MATRIX() { return {1, 0, 0, Complex(0, -1)}; }
 //! T-gate
-matrix_2_2 T_GATE_MATRIX = {
-    COSPI8() - Complex(0, SINPI8()), 0., 0., COSPI8() + Complex(0, SINPI8()) * SINPI8()};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 T_GATE_MATRIX() {
+    return {COSPI8() - Complex(0, SINPI8()), 0., 0., COSPI8() + Complex(0, SINPI8()) * SINPI8()};
+}
 //! Tdag-gate
-matrix_2_2 T_DAG_GATE_MATRIX = {
-    COSPI8() + Complex(0, SINPI8()), 0., 0., COSPI8() - Complex(0, SINPI8())};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 T_DAG_GATE_MATRIX() {
+    return {COSPI8() + Complex(0, SINPI8()), 0., 0., COSPI8() - Complex(0, SINPI8())};
+}
 //! Hadamard gate
-matrix_2_2 HADAMARD_MATRIX = {INVERSE_SQRT2(), INVERSE_SQRT2(), INVERSE_SQRT2(), -INVERSE_SQRT2()};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 HADAMARD_MATRIX() {
+    return {INVERSE_SQRT2(), INVERSE_SQRT2(), INVERSE_SQRT2(), -INVERSE_SQRT2()};
+}
 //! square root of X gate
-matrix_2_2 SQRT_X_GATE_MATRIX = {
-    Complex(0.5, 0.5), Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(0.5, 0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 SQRT_X_GATE_MATRIX() {
+    return {Complex(0.5, 0.5), Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(0.5, 0.5)};
+}
 //! square root of Y gate
-matrix_2_2 SQRT_Y_GATE_MATRIX = {
-    Complex(0.5, 0.5), Complex(-0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 SQRT_Y_GATE_MATRIX() {
+    return {Complex(0.5, 0.5), Complex(-0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5)};
+}
 //! square root dagger of X gate
-matrix_2_2 SQRT_X_DAG_GATE_MATRIX = {
-    Complex(0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5), Complex(0.5, -0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 SQRT_X_DAG_GATE_MATRIX() {
+    return {Complex(0.5, -0.5), Complex(0.5, 0.5), Complex(0.5, 0.5), Complex(0.5, -0.5)};
+}
 //! square root dagger of Y gate
-matrix_2_2 SQRT_Y_DAG_GATE_MATRIX = {
-    Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(-0.5, 0.5), Complex(0.5, -0.5)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 SQRT_Y_DAG_GATE_MATRIX() {
+    return {Complex(0.5, -0.5), Complex(0.5, -0.5), Complex(-0.5, 0.5), Complex(0.5, -0.5)};
+}
 //! Projection to 0
-matrix_2_2 PROJ_0_MATRIX = {1, 0, 0, 0};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PROJ_0_MATRIX() { return {1, 0, 0, 0}; }
 //! Projection to 1
-matrix_2_2 PROJ_1_MATRIX = {0, 0, 0, 1};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PROJ_1_MATRIX() { return {0, 0, 0, 1}; }
 //! complex values for exp(j * i*pi/4 )
-matrix_2_2 PHASE_90ROT = {1., Complex(0, 1), -1, Complex(0, -1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PHASE_90ROT() { return {1., Complex(0, 1), -1, Complex(0, -1)}; }
 //! complex values for exp(-j * i*pi/4 )
-matrix_2_2 PHASE_M90ROT = {1., Complex(0, -1), -1, Complex(0, 1)};
+KOKKOS_INLINE_FUNCTION
+matrix_2_2 PHASE_M90ROT() { return {1., Complex(0, -1), -1, Complex(0, 1)}; }
 }  // namespace qulacs

qulacs/gate/gate.hpp

@@ -7,5 +7,6 @@ namespace qulacs {
 class QuantumGate {
 public:
     virtual void update_quantum_state(StateVector& state_vector) const = 0;
+    virtual ~QuantumGate() {}
 };
 }  // namespace qulacs

qulacs/gate/gate_npair_qubit.cpp

@@ -0,0 +1,9 @@
+#include "gate_npair_qubit.hpp"
+
+#include "update_ops.hpp"
+
+namespace qulacs {
+void FusedSWAP::update_quantum_state(StateVector& state_vector) const {
+    fusedswap_gate(this->qubit_index1, this->qubit_index2, this->block_size, state_vector);
+}
+}  // namespace qulacs

qulacs/gate/gate_npair_qubit.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+#include <cassert>
+
+#include "gate.hpp"
+
+namespace qulacs {
+class FusedSWAP : public QuantumGate {
+    UINT qubit_index1, qubit_index2, block_size;
+
+public:
+    FusedSWAP(UINT _qubit_index1, UINT _qubit_index2, UINT _block_size)
+        : qubit_index1(_qubit_index1), qubit_index2(_qubit_index2), block_size(_block_size) {
+        UINT upper_index = std::max(_qubit_index1, _qubit_index2);
+        UINT lower_index = std::min(_qubit_index1, _qubit_index2);
+        if (upper_index <= (lower_index + block_size - 1)) {
+            throw std::runtime_error(
+                "FusedSwap: upper index must be bigger than lower_index + block_size - 1");
+        }
+    };
+
+    void update_quantum_state(StateVector& state_vector) const override;
+};
+}  // namespace qulacs

qulacs/gate/gate_two_qubit.cpp

@@ -0,0 +1,9 @@
+#include "gate_two_qubit.hpp"
+
+#include "update_ops.hpp"
+
+namespace qulacs {
+void SWAP::update_quantum_state(StateVector& state_vector) const {
+    swap_gate(this->_target1, this->_target2, state_vector);
+}
+}  // namespace qulacs

qulacs/gate/gate_two_qubit.hpp

@@ -0,0 +1,15 @@
+#pragma once
+
+#include <vector>
+
+#include "gate.hpp"
+
+namespace qulacs {
+class SWAP : public QuantumGate {
+    UINT _target1, _target2;
+
+public:
+    SWAP(UINT target1, UINT target2) : _target1(target1), _target2(target2){};
+    void update_quantum_state(StateVector& state_vector) const override;
+};
+}  // namespace qulacs

qulacs/gate/update_ops.hpp

@@ -52,4 +52,20 @@ void single_qubit_dense_matrix_gate(UINT target_qubit_index,
                                     StateVector& state);
 
 void u_gate(UINT target_qubit_index, const matrix_2_2& matrix, StateVector& state);
+
+void swap_gate(UINT target1, UINT target2, StateVector& state);
+
+void fusedswap_gate(UINT target_qubit_index_0,
+                    UINT target_qubit_index_1,
+                    UINT block_size,
+                    StateVector& state);
+
+void pauli_gate(std::vector<UINT> target_qubit_index_list,
+                std::vector<UINT> pauli_id_index,
+                StateVector& state);
+
+void pauli_rotation_gate(std::vector<UINT> target_qubit_index_list,
+                         std::vector<UINT> pauli_id_index,
+                         double angle,
+                         StateVector& state);
 }  // namespace qulacs

qulacs/gate/update_ops_npair_qubit.cpp

@@ -0,0 +1,37 @@
+#include <Kokkos_Core.hpp>
+#include <Kokkos_StdAlgorithms.hpp>
+#include <cassert>
+
+#include "../types.hpp"
+#include "update_ops.hpp"
+
+namespace qulacs {
+void fusedswap_gate(UINT target_qubit_index_0,
+                    UINT target_qubit_index_1,
+                    UINT block_size,
+                    StateVector& state) {
+    UINT n_qubits = state.n_qubits();
+    UINT upper_index, lower_index;
+    upper_index = std::max(target_qubit_index_0, target_qubit_index_1);
+    lower_index = std::min(target_qubit_index_0, target_qubit_index_1);
+    if (n_qubits <= (upper_index + block_size - 1)) {
+        throw std::runtime_error(
+            "FusedSwap: num of qubits must be bigger than upper_index + block_size - 1");
+    }
+    const UINT mask_block = (1 << block_size) - 1;
+    auto amplitudes = state.amplitudes_raw();
+    const UINT kblk_mask = mask_block << upper_index;
+    const UINT jblk_mask = mask_block << lower_index;
+    const UINT else_mask = (1 << n_qubits) - 1 - kblk_mask - jblk_mask;
+
+    Kokkos::parallel_for(
+        1 << n_qubits, KOKKOS_LAMBDA(const UINT& i) {
+            const UINT kblk = (i & kblk_mask) >> upper_index;
+            const UINT jblk = (i & jblk_mask) >> lower_index;
+            if (jblk > kblk) {
+                const UINT index = (i & else_mask) | jblk << upper_index | kblk << lower_index;
+                Kokkos::Experimental::swap(amplitudes[i], amplitudes[index]);
+            }
+        });
+}
+}  // namespace qulacs

qulacs/gate/update_ops_one_target_one_control.cpp → qulacs/gate/update_ops_one_control_one_target.cpp

@@ -12,9 +12,11 @@ void cnot_gate(UINT control_qubit_index, UINT target_qubit_index, StateVector& s
     const UINT control_mask = 1ULL << control_qubit_index;
     auto [min_qubit_index, max_qubit_index] =
         Kokkos::minmax(control_qubit_index, target_qubit_index);
-    const UINT low_mask = (1ULL << min_qubit_index) - 1;
-    const UINT mid_mask = ((1ULL << (max_qubit_index - 1)) - 1) ^ low_mask;
-    const UINT high_mask = ~(low_mask | mid_mask);
+    const UINT min_qubit_mask = 1ULL << min_qubit_index;
+    const UINT max_qubit_mask = 1ULL << (max_qubit_index - 1);
+    const UINT low_mask = min_qubit_mask - 1;
+    const UINT mid_mask = (max_qubit_mask - 1) ^ low_mask;
+    const UINT high_mask = ~(max_qubit_mask - 1);
 
     auto amplitudes = state.amplitudes_raw();
     Kokkos::parallel_for(
@@ -39,7 +41,7 @@ void cz_gate(UINT control_qubit_index, UINT target_qubit_index, StateVector& sta
     auto amplitudes = state.amplitudes_raw();
     Kokkos::parallel_for(
         1ULL << (n_qubits - 2), KOKKOS_LAMBDA(const UINT& it) {
-            UINT i = ((it & high_mask) << 2) | ((it & mid_mask) << 1) | (it & low_mask) |
+            UINT i = (it & low_mask) | ((it & mid_mask) << 1) | ((it & high_mask) << 2) |
                      control_mask | target_mask;
             amplitudes[i] *= -1;
         });

qulacs/gate/update_ops_one_qubit.cpp

@@ -90,27 +90,27 @@ void tdag_gate(UINT target_qubit_index, StateVector& state) {
 }
 
 void sqrtx_gate(UINT target_qubit_index, StateVector& state) {
-    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_X_GATE_MATRIX, state);
+    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_X_GATE_MATRIX(), state);
 }
 
 void sqrtxdag_gate(UINT target_qubit_index, StateVector& state) {
-    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_X_DAG_GATE_MATRIX, state);
+    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_X_DAG_GATE_MATRIX(), state);
 }
 
 void sqrty_gate(UINT target_qubit_index, StateVector& state) {
-    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_Y_GATE_MATRIX, state);
+    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_Y_GATE_MATRIX(), state);
 }
 
 void sqrtydag_gate(UINT target_qubit_index, StateVector& state) {
-    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_Y_DAG_GATE_MATRIX, state);
+    single_qubit_dense_matrix_gate(target_qubit_index, SQRT_Y_DAG_GATE_MATRIX(), state);
 }
 
 void p0_gate(UINT target_qubit_index, StateVector& state) {
-    single_qubit_dense_matrix_gate(target_qubit_index, PROJ_0_MATRIX, state);
+    single_qubit_dense_matrix_gate(target_qubit_index, PROJ_0_MATRIX(), state);
 }
 
 void p1_gate(UINT target_qubit_index, StateVector& state) {
-    single_qubit_dense_matrix_gate(target_qubit_index, PROJ_1_MATRIX, state);
+    single_qubit_dense_matrix_gate(target_qubit_index, PROJ_1_MATRIX(), state);
 }
 
 void rx_gate(UINT target_qubit_index, double angle, StateVector& state) {

qulacs/gate/update_ops_two_qubit.cpp

@@ -0,0 +1,30 @@
+#include <Kokkos_Core.hpp>
+#include <Kokkos_StdAlgorithms.hpp>
+
+#include "../types.hpp"
+#include "constant.hpp"
+#include "update_ops.hpp"
+
+namespace qulacs {
+void swap_gate(UINT target0, UINT target1, StateVector& state) {
+    UINT n_qubits = state.n_qubits();
+    UINT mask_0 = 1ULL << target0;
+    UINT mask_1 = 1ULL << target1;
+    UINT mask = mask_0 + mask_1;
+    UINT min_qubit_index = std::min(target0, target1);
+    UINT max_qubit_index = std::max(target0, target1);
+    UINT min_qubit_mask = 1ULL << min_qubit_index;
+    UINT max_qubit_mask = 1ULL << (max_qubit_index - 1);
+    UINT low_mask = min_qubit_mask - 1;
+    UINT mid_mask = (max_qubit_mask - 1) ^ low_mask;
+    UINT high_mask = ~(max_qubit_mask - 1);
+    auto amplitudes = state.amplitudes_raw();
+    Kokkos::parallel_for(
+        1ULL << (n_qubits - 2), KOKKOS_LAMBDA(const UINT& it) {
+            UINT basis_index_0 =
+                (it & low_mask) + ((it & mid_mask) << 1) + ((it & high_mask) << 2) + mask_0;
+            UINT basis_index_1 = basis_index_0 ^ mask;
+            Kokkos::Experimental::swap(amplitudes[basis_index_0], amplitudes[basis_index_1]);
+        });
+}
+}  // namespace qulacs