test using Armadillo

Author: taminob

CMakeLists.txt

@@ -24,6 +24,24 @@ include(PackageAddTest)
 include(Cache)
 include(AddMQTCoreLibrary)
 
+# include(FetchContent)
+# FetchContent_Declare(
+#   armadillo
+#   GIT_REPOSITORY https://gitlab.com/conradsnicta/armadillo-code.git
+#   GIT_TAG 15.0.1
+# )
+# FetchContent_MakeAvailable(armadillo)
+find_package(Armadillo 15 REQUIRED)
+if(Armadillo_FOUND AND NOT TARGET Armadillo::Armadillo)
+  add_library(Armadillo::Armadillo INTERFACE IMPORTED)
+  set_target_properties(
+    Armadillo::Armadillo
+    PROPERTIES
+      INTERFACE_LINK_LIBRARIES "${ARMADILLO_LIBRARIES}"
+      INTERFACE_INCLUDE_DIRECTORIES "${ARMADILLO_INCLUDE_DIRS}"
+  )
+endif()
+
 option(BUILD_MQT_CORE_BINDINGS "Build the MQT Core Python bindings" OFF)
 if(BUILD_MQT_CORE_BINDINGS)
   # ensure that the BINDINGS option is set

mlir/lib/Dialect/MQTOpt/Transforms/CMakeLists.txt

@@ -7,7 +7,7 @@
 # Licensed under the MIT License
 
 get_property(dialect_libs GLOBAL PROPERTY MLIR_DIALECT_LIBS)
-set(LIBRARIES ${dialect_libs} MQT::CoreIR MQT::CoreDD)
+set(LIBRARIES ${dialect_libs} MQT::CoreIR MQT::CoreDD Armadillo::Armadillo)
 add_compile_options(-fexceptions)
 
 file(GLOB TRANSFORMS_SOURCES *.cpp)

mlir/lib/Dialect/MQTOpt/Transforms/GateDecompositionPattern.cpp

@@ -25,7 +25,7 @@
 #include <mlir/Support/LogicalResult.h>
 #include <string>
 
-#include "b.h"
+#include "d.h"
 
 namespace mqt::ir::opt {
 
@@ -324,7 +324,7 @@ struct GateDecompositionPattern final
 
     auto [U, S] = self_adjoint_evd(A);
 
-    return {U, S};
+    return std::make_pair(U, S);
   }
 
   static std::tuple<matrix2x2, matrix2x2, fp>
@@ -564,6 +564,15 @@ struct GateDecompositionPattern final
       llvm::errs() << "===== M2 =====\n";
       helpers::print(m2);
 
+      arma::Mat<qfp> U(4, 4);
+  for (int i = 0; i < 4; ++i) {
+    for (int j = 0; j < 4; ++j) {
+      U.at(j, i) = u_p[j * 4 + i];
+    }
+  }
+      auto x = U.st() * U;
+      std::cerr << "ARMA\n" << U.t() << "\n\n" << U << "\n\n" << x << std::endl;
+
       // M2 is a symmetric complex matrix. We need to decompose it as M2 = P D
       // P^T where P ∈ SO(4), D is diagonal with unit-magnitude elements.
       //
@@ -603,7 +612,7 @@ struct GateDecompositionPattern final
         llvm::transform(m2, m2_real.begin(), [&](const qfp& val) {
           return rand_a * val.real() + rand_b * val.imag();
         });
-        auto p_inner_real = self_adjoint_eigen_lower(m2_real).first;
+        rmatrix4x4 p_inner_real = self_adjoint_eigen_lower(m2_real).first;
         matrix4x4 p_inner;
         llvm::transform(p_inner_real, p_inner.begin(),
                         [](auto&& x) { return qfp(x, 0.0); });
@@ -638,6 +647,7 @@ struct GateDecompositionPattern final
       std::array<fp, d.size()> d_real;
       llvm::transform(d, d_real.begin(),
                       [](auto&& x) { return -std::arg(x) / 2.0; });
+      helpers::print(d_real, "D_REAL");
       d_real[3] = -d_real[0] - d_real[1] - d_real[2];
       std::array<fp, 3> cs;
       for (std::size_t i = 0; i < cs.size(); ++i) {
@@ -680,6 +690,7 @@ struct GateDecompositionPattern final
       temp[1 * 4 + 1] = std::exp(IM * d_real[1]);
       temp[2 * 4 + 2] = std::exp(IM * d_real[2]);
       temp[3 * 4 + 3] = std::exp(IM * d_real[3]);
+      helpers::print(temp, "TEMP");
       auto k1 = magic_basis_transform(dot(dot(u_p, p), temp),
                                       MagicBasisTransform::Into);
       auto k2 = magic_basis_transform(transpose(p), MagicBasisTransform::Into);

mlir/lib/Dialect/MQTOpt/Transforms/Helpers.h

@@ -40,10 +40,27 @@ constexpr qfp M_IM{0., -1.};
 namespace mqt::ir::opt::helpers {
 
 // TODO: remove
-template <std::size_t N> void print(std::array<std::complex<fp>, N> matrix) {
+template <std::size_t N> void print(std::array<std::complex<fp>, N> matrix, std::string s = "") {
   int i{};
+  if (!s.empty()) {
+    llvm::errs() << "=== " << s << " ===\n";
+  }
+  for (auto&& a : matrix) {
+    std::cerr << std::setprecision(17) << a.real() << 'i' << a.imag() << ' ';
+    if (++i % 4 == 0) {
+      llvm::errs() << '\n';
+    }
+  }
+  llvm::errs() << '\n';
+}
+
+template <std::size_t N> void print(std::array<fp, N> matrix, std::string s = "") {
+  int i{};
+  if (!s.empty()) {
+    llvm::errs() << "=== " << s << " ===\n";
+  }
   for (auto&& a : matrix) {
-    std::cerr << std::setprecision(50) << a.real() << 'i' << a.imag() << ' ';
+    std::cerr << std::setprecision(17) << a << ' ';
     if (++i % 4 == 0) {
       llvm::errs() << '\n';
     }

mlir/lib/Dialect/MQTOpt/Transforms/b.h

@@ -425,7 +425,7 @@ void householderSequenceEval(rmatrix4x4& m_vectors,
 }
 
 // https://docs.rs/faer/latest/faer/mat/generic/struct.Mat.html#method.self_adjoint_eigen
-static std::pair<rmatrix4x4, rdiagonal4x4> self_adjoint_evd(rmatrix4x4 matrix) {
+inline std::pair<rmatrix4x4, rdiagonal4x4> self_adjoint_evd(rmatrix4x4 matrix) {
 
   auto lowerTriangularView = [](auto matrix) {
     const int n = std::sqrt(matrix.size());

mlir/lib/Dialect/MQTOpt/Transforms/d.h

@@ -0,0 +1,31 @@
+#pragma once
+
+#include "Helpers.h"
+
+#include <armadillo>
+
+namespace mqt::ir::opt {
+
+auto self_adjoint_evd(rmatrix4x4 A) {
+  arma::Mat<fp> a(4, 4, arma::fill::scalar_holder<fp>(0.0));
+  for (int i = 0; i < 4; ++i) {
+    for (int j = 0; j < 4; ++j) {
+      a.at(j, i) = A[j * 4 + i];
+    }
+  }
+  arma::Mat<fp> vecs;
+  arma::Col<fp> vals;
+  arma::eig_sym(vals, vecs, a, "std");
+  rmatrix4x4 rvecs;
+  rdiagonal4x4 rvals;
+  for (int i = 0; i < 4; ++i) {
+    for (int j = 0; j < 4; ++j) {
+      rvecs[j * 4 + i] = vecs.at(j, i);
+    }
+    rvals[i] = vals.at(i);
+  }
+  std::cerr << "========\n" << vecs << "========\n" << std::endl;
+  std::cerr << "========\n" << vals << "========\n" << std::endl;
+  return std::make_pair(rvecs, rvals);
+}
+} // namespace mqt::ir::opt