✨ Add complex eigendecomposition for QCO matrices

mlir/include/mlir/Dialect/QCO/Utils/Eigensolver.h

@@ -0,0 +1,116 @@
+/*
+ * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+ * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+ * All rights reserved.
+ *
+ * SPDX-License-Identifier: MIT
+ *
+ * Licensed under the MIT License
+ */
+
+#pragma once
+
+#include "mlir/Dialect/QCO/Utils/Matrix.h"
+
+#include <array>
+#include <optional>
+
+namespace mlir::qco {
+
+/**
+ * @brief Computes the eigendecomposition of a real symmetric `4x4` matrix.
+ *
+ * Uses Householder tridiagonalization (EISPACK `tred2`) followed by implicit
+ * QL iteration (`tql2`) on the tridiagonal form. Adapted from John Burkardt's
+ * MIT-licensed EISPACK C port (`tred2`, `tql2`):
+ * https://people.sc.fsu.edu/~jburkardt/c_src/eispack/eispack.c
+ * Original Fortran: https://netlib.org/eispack/tred2.f,
+ * https://netlib.org/eispack/tql2.f
+ *
+ * @param symmetric Row-major real symmetric `4x4` matrix.
+ * @return Ascending eigenvalues and matching eigenvectors (as columns).
+ */
+[[nodiscard]] SymmetricEigen4
+decomposeSymmetricEigen4(const std::array<double, 16>& symmetric);
+
+/**
+ * @brief Computes the eigendecomposition of a real symmetric `4x4` matrix.
+ *
+ * @copydoc decomposeSymmetricEigen4(const std::array<double, 16>&)
+ *
+ * @param matrix Source matrix.
+ */
+[[nodiscard]] SymmetricEigen4 decomposeSymmetricEigen4(const Matrix4x4& matrix);
+
+/**
+ * @brief Computes the eigendecomposition of a `1x1` dynamic matrix.
+ *
+ * @param matrix Source matrix.
+ * @return The single eigenpair.
+ */
+[[nodiscard]] std::optional<ComplexEigen>
+decomposeComplexEigen1x1(const DynamicMatrix& matrix);
+
+/**
+ * @brief Computes the eigendecomposition of a `2x2` complex matrix using a
+ * closed-form formula.
+ *
+ * @param matrix Source matrix.
+ * @return Eigenpairs, or `std::nullopt` if the closed-form solver produces
+ * non-finite eigenvalues.
+ */
+[[nodiscard]] std::optional<ComplexEigen2>
+decomposeComplexEigen2(const Matrix2x2& matrix);
+
+/**
+ * @brief Computes the eigendecomposition of a `4x4` complex matrix.
+ *
+ * Uses stack-specialized EISPACK `corth` / `comqr2` for fixed `n = 4`
+ * (complex Hessenberg reduction and QR eigenanalysis). `pythag` and `csroot`
+ * follow John Burkardt's MIT-licensed EISPACK C port; `cdiv`, `corth`, and
+ * `comqr2` follow NETLIB EISPACK Fortran
+ * (https://netlib.org/eispack/cdiv.f, https://netlib.org/eispack/corth.f,
+ * https://netlib.org/eispack/comqr2.f). See also
+ * https://people.sc.fsu.edu/~jburkardt/c_src/eispack/eispack.c
+ *
+ * @param matrix Source matrix.
+ * @return Eigenpairs, or `std::nullopt` if the solver does not converge.
+ */
+[[nodiscard]] std::optional<ComplexEigen4>
+decomposeComplexEigen4(const Matrix4x4& matrix);
+
+/**
+ * @brief Computes the eigendecomposition of a square dynamic matrix.
+ *
+ * Uses EISPACK `corth` followed by `comqr2` for dimensions other than `1`, `2`,
+ * and `4` (which have specialized paths in @ref DynamicMatrix::complexEigen).
+ * `pythag` and `csroot` follow John Burkardt's MIT-licensed EISPACK C port;
+ * `cdiv`, `corth`, and `comqr2` follow NETLIB EISPACK Fortran
+ * (https://netlib.org/eispack/cdiv.f, https://netlib.org/eispack/corth.f,
+ * https://netlib.org/eispack/comqr2.f). See also
+ * https://people.sc.fsu.edu/~jburkardt/c_src/eispack/eispack.c
+ *
+ * @param matrix Square source matrix.
+ * @return Eigenpairs, or `std::nullopt` if the matrix is not square, its
+ * dimension exceeds `INT_MAX`, or the solver does not converge.
+ */
+[[nodiscard]] std::optional<ComplexEigen>
+decomposeComplexEigenDynamic(const DynamicMatrix& matrix);
+
+/**
+ * @brief Lifts a fixed `2x2` eigendecomposition to @ref ComplexEigen.
+ *
+ * @param eigen2 Fixed-size eigenpairs from @ref decomposeComplexEigen2.
+ * @return Dynamic-matrix eigenvector storage with the same eigenvalues.
+ */
+[[nodiscard]] ComplexEigen fromComplexEigen(const ComplexEigen2& eigen2);
+
+/**
+ * @brief Lifts a fixed `4x4` eigendecomposition to @ref ComplexEigen.
+ *
+ * @param eigen4 Fixed-size eigenpairs from @ref decomposeComplexEigen4.
+ * @return Dynamic-matrix eigenvector storage with the same eigenvalues.
+ */
+[[nodiscard]] ComplexEigen fromComplexEigen(const ComplexEigen4& eigen4);
+
+} // namespace mlir::qco

mlir/include/mlir/Dialect/QCO/Utils/Matrix.h

@@ -15,9 +15,11 @@
 
 #include <array>
 #include <complex>
+#include <concepts>
 #include <cstddef>
 #include <cstdint>
 #include <memory>
+#include <optional>
 #include <type_traits>
 
 namespace mlir::qco {
@@ -31,6 +33,9 @@ inline constexpr double MATRIX_TOLERANCE = 1e-14;
 class DynamicMatrix;
 struct Matrix4x4;
 struct SymmetricEigen4;
+struct ComplexEigen;
+struct ComplexEigen2;
+struct ComplexEigen4;
 
 /**
  * @brief 1x1 matrix for global-phase gates.
@@ -238,6 +243,14 @@ struct Matrix2x2 {
    */
   [[nodiscard]] bool assignFrom(const DynamicMatrix& src);
 
+  /**
+   * @brief Computes the eigendecomposition of this complex matrix.
+   *
+   * @return Eigenpairs, or `std::nullopt` if the closed-form solver produces
+   * non-finite eigenvalues.
+   */
+  [[nodiscard]] std::optional<ComplexEigen2> complexEigen() const;
+
   /**
    * @brief Embed this single-qubit matrix into an @p numQubits-qubit Hilbert
    * space.
@@ -536,6 +549,13 @@ struct Matrix4x4 {
    */
   [[nodiscard]] static SymmetricEigen4
   symmetricEigen4(const std::array<double, 16>& symmetric);
+
+  /**
+   * @brief Computes the eigendecomposition of this complex matrix.
+   *
+   * @return Eigenpairs, or `std::nullopt` if the solver does not converge.
+   */
+  [[nodiscard]] std::optional<ComplexEigen4> complexEigen() const;
 };
 
 /**
@@ -756,11 +776,36 @@ class DynamicMatrix {
    */
   [[nodiscard]] bool isIdentity(double tol = MATRIX_TOLERANCE) const;
 
+  /**
+   * @brief Computes the eigendecomposition of this square matrix.
+   *
+   * Uses EISPACK `corth` followed by `comqr2` (complex Hessenberg reduction
+   * and QR eigenanalysis). `pythag` and `csroot` follow John Burkardt's
+   * MIT-licensed EISPACK C port; `cdiv`, `corth`, and `comqr2` follow NETLIB
+   * EISPACK Fortran (https://netlib.org/eispack/cdiv.f,
+   * https://netlib.org/eispack/corth.f, https://netlib.org/eispack/comqr2.f).
+   * Dimensions `1` and `2` use closed-form formulas;
+   * dimension `4` delegates to @ref Matrix4x4::complexEigen; all other sizes
+   * use the general EISPACK path.
+   *
+   * @return Eigenpairs, or `std::nullopt` if this matrix is empty or the
+   * solver does not converge.
+   */
+  [[nodiscard]] std::optional<ComplexEigen> complexEigen() const;
+
 private:
   struct Impl;
   std::unique_ptr<Impl> impl_;
 };
 
+/**
+ * @brief Concept for the four supported matrix types.
+ */
+template <typename T>
+concept SupportedMatrix =
+    std::same_as<T, Matrix1x1> || std::same_as<T, Matrix2x2> ||
+    std::same_as<T, Matrix4x4> || std::same_as<T, DynamicMatrix>;
+
 /**
  * @brief Type trait for the four supported matrix types.
  *
@@ -772,9 +817,7 @@ class DynamicMatrix {
 template <typename T>
 inline constexpr bool
     is_supported_matrix_v = // NOLINT(readability-identifier-naming)
-    std::disjunction_v<std::is_same<T, Matrix1x1>, std::is_same<T, Matrix2x2>,
-                       std::is_same<T, Matrix4x4>,
-                       std::is_same<T, DynamicMatrix>>;
+    SupportedMatrix<T>;
 
 /// Scalar-on-the-left multiply `scalar * matrix` (commutes with the member
 /// `matrix * scalar`). Provided so generic code can scale a matrix from
@@ -802,4 +845,43 @@ struct SymmetricEigen4 {
   Matrix4x4 eigenvectors{};
 };
 
+/**
+ * @brief Eigenvalues and eigenvectors of a complex `4x4` matrix.
+ *
+ * `eigenvalues[i]` matches column `i` of `eigenvectors` (column `j` is the
+ * eigenvector for `eigenvalues[j]`).
+ */
+struct ComplexEigen4 {
+  /// Eigenvalues in no particular order.
+  std::array<Complex, 4> eigenvalues{};
+  /// Eigenvectors as columns (column `j` matches `eigenvalues[j]`).
+  Matrix4x4 eigenvectors{};
+};
+
+/**
+ * @brief Eigenvalues and eigenvectors of a complex `2x2` matrix.
+ *
+ * `eigenvalues[i]` matches column `i` of `eigenvectors` (column `j` is the
+ * eigenvector for `eigenvalues[j]`).
+ */
+struct ComplexEigen2 {
+  /// Eigenvalues in no particular order.
+  std::array<Complex, 2> eigenvalues{};
+  /// Eigenvectors as columns (column `j` matches `eigenvalues[j]`).
+  Matrix2x2 eigenvectors{};
+};
+
+/**
+ * @brief Eigenvalues and eigenvectors of a square complex matrix.
+ *
+ * `eigenvalues[i]` matches column `i` of `eigenvectors` (column `j` is the
+ * eigenvector for `eigenvalues[j]`).
+ */
+struct ComplexEigen {
+  /// Eigenvalues in no particular order.
+  SmallVector<Complex, 8> eigenvalues;
+  /// Eigenvectors as columns (column `j` matches `eigenvalues[j]`).
+  DynamicMatrix eigenvectors;
+};
+
 } // namespace mlir::qco

mlir/lib/Dialect/QCO/Utils/CMakeLists.txt

@@ -6,15 +6,29 @@
 #
 # Licensed under the MIT License
 
-add_mlir_library(MLIRQCOMatrix PARTIAL_SOURCES_INTENDED Matrix.cpp LINK_LIBS PUBLIC MLIRSupport)
+add_mlir_library(
+  MLIRQCOMatrix
+  PARTIAL_SOURCES_INTENDED
+  Matrix.cpp
+  Eigensolver.cpp
+  LINK_LIBS
+  PUBLIC
+  MLIRSupport)
 
 mqt_mlir_target_use_project_options(MLIRQCOMatrix)
 
-target_sources(MLIRQCOMatrix PUBLIC FILE_SET HEADERS BASE_DIRS ${MQT_MLIR_SOURCE_INCLUDE_DIR} FILES
-                                    ${MQT_MLIR_SOURCE_INCLUDE_DIR}/mlir/Dialect/QCO/Utils/Matrix.h)
+target_sources(
+  MLIRQCOMatrix
+  PUBLIC FILE_SET
+         HEADERS
+         BASE_DIRS
+         ${MQT_MLIR_SOURCE_INCLUDE_DIR}
+         FILES
+         ${MQT_MLIR_SOURCE_INCLUDE_DIR}/mlir/Dialect/QCO/Utils/Matrix.h
+         ${MQT_MLIR_SOURCE_INCLUDE_DIR}/mlir/Dialect/QCO/Utils/Eigensolver.h)
 
 file(GLOB_RECURSE UTILS_CPP "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp")
-list(FILTER UTILS_CPP EXCLUDE REGEX "Matrix\\.cpp$")
+list(FILTER UTILS_CPP EXCLUDE REGEX "(Matrix|Eigensolver)\\.cpp$")
 
 add_mlir_dialect_library(
   MLIRQCOUtils
@@ -34,7 +48,7 @@ mqt_mlir_target_use_project_options(MLIRQCOUtils)
 
 # collect header files
 file(GLOB_RECURSE UTILS_HEADERS_SOURCE "${MQT_MLIR_SOURCE_INCLUDE_DIR}/mlir/Dialect/QCO/Utils/*.h")
-list(FILTER UTILS_HEADERS_SOURCE EXCLUDE REGEX "Matrix\\.h$")
+list(FILTER UTILS_HEADERS_SOURCE EXCLUDE REGEX "(Matrix|Eigensolver)\\.h$")
 file(GLOB_RECURSE UTILS_HEADERS_BUILD "${MQT_MLIR_BUILD_INCLUDE_DIR}/mlir/Dialect/QCO/Utils/*.inc")
 
 # add public headers using file sets