[on hold] Switch convert_2q_block_matrix.rs to use matmul from faer and an explicit version of kron

crates/accelerate/src/convert_2q_block_matrix.rs

@@ -15,11 +15,16 @@ use pyo3::wrap_pyfunction;
 use pyo3::Python;
 
 use num_complex::Complex64;
-use numpy::ndarray::linalg::kron;
 use numpy::ndarray::{aview2, Array2, ArrayView2};
 use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2};
 use smallvec::SmallVec;
 
+use faer::modules::core::mul::matmul;
+use faer::perm::swap_rows;
+use faer::prelude::*;
+use faer::{Mat, Parallelism};
+use faer_ext::{IntoFaerComplex, IntoNdarrayComplex};
+
 static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [
     [Complex64::new(1., 0.), Complex64::new(0., 0.)],
     [Complex64::new(0., 0.), Complex64::new(1., 0.)],
@@ -32,32 +37,61 @@ pub fn blocks_to_matrix(
     py: Python,
     op_list: Vec<(PyReadonlyArray2<Complex64>, SmallVec<[u8; 2]>)>,
 ) -> PyResult<Py<PyArray2<Complex64>>> {
-    let identity = aview2(&ONE_QUBIT_IDENTITY);
-    let input_matrix = op_list[0].0.as_array();
-    let mut matrix: Array2<Complex64> = match op_list[0].1.as_slice() {
-        [0] => kron(&identity, &input_matrix),
-        [1] => kron(&input_matrix, &identity),
+    let identity = aview2(&ONE_QUBIT_IDENTITY).into_faer_complex();
+    let input_matrix = op_list[0].0.as_array().into_faer_complex();
+
+    let mut matrix = match op_list[0].1.as_slice() {
+        [0] => identity.kron(input_matrix.as_ref()),
+        [1] => input_matrix.kron(identity.as_ref()),
         [0, 1] => input_matrix.to_owned(),
-        [1, 0] => change_basis(input_matrix),
-        [] => Array2::eye(4),
+        [1, 0] => change_basis_faer(input_matrix),
+        [] => Mat::<c64>::identity(4, 4),
         _ => unreachable!(),
     };
     for (op_matrix, q_list) in op_list.into_iter().skip(1) {
-        let op_matrix = op_matrix.as_array();
+        let op_matrix = op_matrix.as_array().into_faer_complex();
 
         let result = match q_list.as_slice() {
-            [0] => Some(kron(&identity, &op_matrix)),
-            [1] => Some(kron(&op_matrix, &identity)),
-            [1, 0] => Some(change_basis(op_matrix)),
-            [] => Some(Array2::eye(4)),
-            _ => None,
-        };
-        matrix = match result {
-            Some(result) => result.dot(&matrix),
-            None => op_matrix.dot(&matrix),
+            [0] => identity.kron(op_matrix.as_ref()),
+            [1] => op_matrix.kron(identity.as_ref()),
+            [1, 0] => change_basis_faer(op_matrix),
+            [] => Mat::<c64>::identity(4, 4),
+            _ => op_matrix.to_owned(),
         };
+
+        let aux = matrix.clone();
+        matmul(
+            matrix.as_mut(),
+            result.as_ref(),
+            aux.as_ref(),
+            None,
+            c64::new(1., 0.),
+            Parallelism::None,
+        );
     }
-    Ok(matrix.into_pyarray_bound(py).unbind())
+
+    Ok(matrix
+        .as_ref()
+        .into_ndarray_complex()
+        .to_owned()
+        .into_pyarray_bound(py)
+        .unbind())
+}
+
+/// Switches the order of qubits in a two qubit operation.
+/// This function will substitue `change_basis` once the
+/// `two_qubit_decompose.rs` uses Mat<c64> instead of ArrayView2
+#[inline]
+pub fn change_basis_faer(matrix: MatRef<c64>) -> Mat<c64> {
+    let mut trans_matrix: Mat<c64> = matrix.transpose().to_owned();
+    let (row1, row2) = trans_matrix.as_mut().two_rows_mut(1, 2);
+    swap_rows(row1, row2);
+
+    trans_matrix = trans_matrix.transpose().to_owned();
+    let (row1, row2) = trans_matrix.as_mut().two_rows_mut(1, 2);
+    swap_rows(row1, row2);
+
+    trans_matrix
 }
 
 /// Switches the order of qubits in a two qubit operation.