optimize transform_dense_matrix_by_order, fix macos warning, add tran…

…sform_sparse_matrix_by_order,

exe/main.cpp

@@ -35,6 +35,15 @@ void run() {
 
     auto pprob_gate = gate::ParamProbablistic({0.7, 0.3}, {prx_gate, pry_gate});
     std::cout << pprob_gate << std::endl;
+
+    Eigen::Matrix<StdComplex, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> mat(4, 4);
+    mat << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15;
+
+    auto dense_gate = gate::DenseMatrix({1, 3}, mat);
+    std::cout << dense_gate << std::endl;
+
+    auto sparse_gate = gate::SparseMatrix({2, 0}, mat.sparseView());
+    std::cout << sparse_gate << std::endl;
 }
 
 int main() {

scaluq/gate/gate_matrix.hpp

@@ -134,7 +134,7 @@ class DenseMatrixGateImpl : public GateBase {
         if (_is_unitary) {
             inv_eigen = mat_eigen.adjoint();
         } else {
-            inv_eigen = mat_eigen.inverse();
+            inv_eigen = mat_eigen.inverse().eval();
         }
         return std::make_shared<const DenseMatrixGateImpl>(
             _target_mask, _control_mask, inv_eigen, _is_unitary);
@@ -165,17 +165,27 @@ class DenseMatrixGateImpl : public GateBase {
 
 class SparseMatrixGateImpl : public GateBase {
     SparseMatrix _matrix;
+    std::uint64_t num_nnz;
 
 public:
     SparseMatrixGateImpl(std::uint64_t target_mask,
                          std::uint64_t control_mask,
                          const SparseComplexMatrix& mat)
-        : GateBase(target_mask, control_mask), _matrix(SparseMatrix(mat)) {}
+        : GateBase(target_mask, control_mask),
+          _matrix(SparseMatrix(mat)),
+          num_nnz(mat.nonZeros()) {}
 
     Gate get_inverse() const override {
-        ComplexMatrix mat_eigen = convert_coo_to_external_matrix(_matrix);
-        ComplexMatrix inv_eigen = mat_eigen.inverse().eval();
-        return std::make_shared<const DenseMatrixGateImpl>(_target_mask, _control_mask, inv_eigen);
+        Kokkos::View<SparseValue*, Kokkos::HostSpace> vec_h("h_view", num_nnz);
+        Kokkos::deep_copy(vec_h, _matrix._values);
+        // conversion to Eigen matrix (COO format)
+        ComplexMatrix eigen_matrix(_matrix._row, _matrix._col);
+        for (std::size_t i = 0; i < vec_h.extent(0); i++) {
+            eigen_matrix(vec_h(i).r, vec_h(i).c) = vec_h(i).val;
+        }
+
+        return std::make_shared<const DenseMatrixGateImpl>(
+            _target_mask, _control_mask, eigen_matrix.inverse().eval());
     }
 
     Matrix get_matrix_internal() const {

scaluq/gate/update_ops_sparse_matrix.cpp

@@ -30,7 +30,9 @@ void sparse_matrix_gate(std::uint64_t target_mask,
     Kokkos::parallel_for(
         "COO_Update",
         Kokkos::MDRangePolicy<Kokkos::Rank<2>>(
-            {0, 0}, {state.dim() >> std::popcount(target_mask | control_mask), values.size()}),
+            {0, 0},
+            {static_cast<std::int64_t>(state.dim() >> std::popcount(target_mask | control_mask)),
+             static_cast<std::int64_t>(values.size())}),
         KOKKOS_LAMBDA(std::uint64_t outer, std::uint64_t inner) {
             std::uint64_t basis =
                 internal::insert_zero_at_mask_positions(outer, target_mask | control_mask) |

scaluq/util/utility.hpp

@@ -197,26 +197,37 @@ inline ComplexMatrix transform_dense_matrix_by_order(const ComplexMatrix& mat,
             std::lower_bound(sorted.begin(), sorted.end(), targets[i]) - sorted.begin();
     }
 
-    auto transform_index = [&targets_order](std::size_t index) {
-        std::size_t transformed = 0;
+    // transform_indices
+    std::vector<std::uint64_t> transformed(targets_order.size());
+    for (std::size_t index = 0; index < targets_order.size(); index++) {
         for (std::size_t j = 0; j < targets_order.size(); j++) {
-            transformed |= ((index & (1ULL << targets_order[j])) >> targets_order[j]) << j;
+            transformed[j] |= ((index & (1ULL << targets_order[j])) >> targets_order[j]) << j;
         }
-        return transformed;
-    };
+    }
 
     ComplexMatrix ret(matrix_size, matrix_size);
 
     for (std::size_t i = 0; i < matrix_size; i++) {
-        std::size_t row_dst = transform_index(i);
+        std::size_t row_dst = transformed[i];
         for (std::size_t j = 0; j < matrix_size; j++) {
-            std::size_t col_dst = transform_index(j);
+            std::size_t col_dst = transformed[j];
             ret(row_dst, col_dst) = mat(i, j);
         }
     }
     return ret;
 }
 
+inline SparseComplexMatrix transform_sparse_matrix_by_order(
+    // This is temporary implementation.
+    // SparseComplexMatrix will be replaced with std::vector<std::vector<std::Complex<double>>>
+    // hence this function will be refactored.
+    const SparseComplexMatrix& mat,
+    const std::vector<std::uint64_t>& targets) {
+    ComplexMatrix dense_mat = mat;
+    ComplexMatrix transformed = transform_dense_matrix_by_order(dense_mat, targets);
+    return transformed.sparseView();
+}
+
 }  // namespace internal
 
 }  // namespace scaluq

tests/gate/gate_test.cpp

@@ -498,9 +498,9 @@ void run_random_gate_apply_general_dense(std::uint64_t n_qubits) {
             std::shuffle(index_list.begin(), index_list.end(), engine);
             targets[0] = index_list[0];
             targets[1] = index_list[1];
-            if (targets[0] > targets[1]) {
-                std::swap(targets[0], targets[1]);
-            }
+            // if (targets[0] > targets[1]) {
+            //     std::swap(targets[0], targets[1]);
+            // }
             Umerge = internal::kronecker_product(U2, U1);
             internal::ComplexMatrix mat(Umerge.rows(), Umerge.cols());
             test_state = get_expanded_eigen_matrix_with_identity(targets[1], U2, n_qubits) *