add function for transforming matrix, maybe fix cuda warning

qulacs · Oct 8, 2024 · bb2219b · bb2219b
1 parent 88cf10f
commit bb2219b
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 16 deletions.
diff --git a/scaluq/gate/gate_factory.hpp b/scaluq/gate/gate_factory.hpp
@@ -174,8 +174,20 @@ inline Gate DenseMatrix(const std::vector<std::uint64_t>& targets,
             "gate::DenseMatrix(const std::vector<std::uint64_t>&, const internal::ComplexMatrix&): "
             "matrix size must be 2^{n_qubits} x 2^{n_qubits}.");
     }
+    if (std::is_sorted(targets.begin(), targets.end())) {
+        return internal::GateFactory::create_gate<internal::DenseMatrixGateImpl>(
+            internal::vector_to_mask(targets),
+            internal::vector_to_mask(controls),
+            matrix,
+            is_unitary);
+    }
+    internal::ComplexMatrix matrix_transformed =
+        internal::transform_dense_matrix_by_order(matrix, targets);
     return internal::GateFactory::create_gate<internal::DenseMatrixGateImpl>(
-        internal::vector_to_mask(targets), internal::vector_to_mask(controls), matrix, is_unitary);
+        internal::vector_to_mask(targets),
+        internal::vector_to_mask(controls),
+        matrix_transformed,
+        is_unitary);
 }
 inline Gate SparseMatrix(const std::vector<std::uint64_t>& targets,
                          const internal::SparseComplexMatrix& matrix,

diff --git a/scaluq/gate/gate_matrix.hpp b/scaluq/gate/gate_matrix.hpp
@@ -174,8 +174,7 @@ class SparseMatrixGateImpl : public GateBase {
 
     Gate get_inverse() const override {
         ComplexMatrix mat_eigen = convert_coo_to_external_matrix(_matrix);
-        ComplexMatrix inv_eigen;
-        inv_eigen = mat_eigen.inverse();
+        ComplexMatrix inv_eigen = mat_eigen.inverse().eval();
         return std::make_shared<const DenseMatrixGateImpl>(_target_mask, _control_mask, inv_eigen);
     }
 

diff --git a/scaluq/util/utility.hpp b/scaluq/util/utility.hpp
@@ -184,24 +184,37 @@ inline ComplexMatrix convert_coo_to_external_matrix(SparseMatrix mat) {
     return eigen_matrix;
 }
 
-void std::vector<std::uint64_t> compute_mapping(std::vector<std::uint64_t> targets) {
-    std::vector<std::uint64_t> ret;
+inline ComplexMatrix transform_dense_matrix_by_order(const ComplexMatrix& mat,
+                                                     const std::vector<std::uint64_t>& targets) {
+    std::vector<std::uint64_t> sorted(targets);
+    std::sort(sorted.begin(), sorted.end());
+
+    const std::size_t matrix_size = mat.rows();
+
+    std::vector<std::uint64_t> targets_order(targets.size());
     for (std::size_t i = 0; i < targets.size(); i++) {
-        ret.emplace_back(targets[i]);
+        targets_order[i] =
+            std::lower_bound(sorted.begin(), sorted.end(), targets[i]) - sorted.begin();
     }
-    return ret;
-}
 
-void transform_dense_matrix_by_order(ComplexMatrix& mat, std::vector<std::uint64_t> targets) {
-    std::vector<std::uint64_t> mapping = compute_mapping(targets);
-    std::vector<std::uint64_t> prev_position = std::iota(0, targets.size());
-    for (std::size_t src = 0; src < targets.size(); src++) {
-        std::size_t dst = mapping[src];
-        if (src != dst) {
-            mat.row(src).swap(mat.row(dst));
-            std::swap(prev_position[src], prev_position[dst]);
+    auto transform_index = [&targets_order](std::size_t index) {
+        std::size_t transformed = 0;
+        for (std::size_t j = 0; j < targets_order.size(); j++) {
+            transformed |= ((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);
+        for (std::size_t j = 0; j < matrix_size; j++) {
+            std::size_t col_dst = transform_index(j);
+            ret(row_dst, col_dst) = mat(i, j);
         }
     }
+    return ret;
 }
 
 }  // namespace internal