From 239ed3452bc85a0102611bc58c416809513768c5 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 10:50:40 +0200 Subject: [PATCH 01/10] =?UTF-8?q?=E2=9C=A8=20Enhance=20Matrix=20functional?= =?UTF-8?q?ity=20with=20new=20methods=20and=20optimizations?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/include/mlir/Dialect/QCO/Utils/Matrix.h | 196 ++++++ mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 605 +++++++++++++++++- .../Dialect/QCO/Utils/test_unitary_matrix.cpp | 277 +++++++- 3 files changed, 1046 insertions(+), 32 deletions(-) diff --git a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h index 8f49d372d7..402e86ac4a 100644 --- a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h +++ b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h @@ -183,6 +183,12 @@ struct Matrix2x2 { */ [[nodiscard]] Matrix2x2 adjoint() const; + /** + * @brief Returns the (non-conjugate) transpose of this matrix. + * @return Transposed matrix `A^T`. + */ + [[nodiscard]] Matrix2x2 transpose() const; + /** * @brief Returns the trace of this matrix. * @return Sum of diagonal entries. @@ -195,6 +201,13 @@ struct Matrix2x2 { */ [[nodiscard]] Complex determinant() const; + /** + * @brief Checks whether this matrix is approximately the identity. + * @param tol Maximum allowed complex modulus of each entry difference. + * @return True if every entry is within @p tol of the identity. + */ + [[nodiscard]] bool isIdentity(double tol = MATRIX_TOLERANCE) const; + /** * @brief Checks approximate equality using an absolute entry-wise tolerance. * @@ -322,6 +335,12 @@ struct Matrix4x4 { */ [[nodiscard]] Matrix4x4 adjoint() const; + /** + * @brief Returns the (non-conjugate) transpose of this matrix. + * @return Transposed matrix `A^T`. + */ + [[nodiscard]] Matrix4x4 transpose() const; + /** * @brief Returns the trace of this matrix. * @return Sum of diagonal entries. @@ -334,6 +353,53 @@ struct Matrix4x4 { */ [[nodiscard]] Complex determinant() const; + /** + * @brief Checks whether this matrix is approximately the identity. + * @param tol Maximum allowed complex modulus of each entry difference. + * @return True if every entry is within @p tol of the identity. + */ + [[nodiscard]] bool isIdentity(double tol = MATRIX_TOLERANCE) const; + + /** + * @brief Returns the four diagonal entries `(m00, m11, m22, m33)`. + * @return Array of diagonal entries. + */ + [[nodiscard]] std::array diagonal() const; + + /** + * @brief Builds a diagonal matrix from four diagonal entries. + * @param diagonalEntries Diagonal entries `(m00, m11, m22, m33)`. + * @return Diagonal matrix with the given entries. + */ + [[nodiscard]] static Matrix4x4 + fromDiagonal(const std::array& diagonalEntries); + + /** + * @brief Returns the entries of column @p col, top to bottom. + * @param col Column index in `[0, K_COLS)`. + * @return Array of the four column entries. + */ + [[nodiscard]] std::array column(std::size_t col) const; + + /** + * @brief Overwrites column @p col with @p values. + * @param col Column index in `[0, K_COLS)`. + * @param values New column entries, top to bottom. + */ + void setColumn(std::size_t col, const std::array& values); + + /** + * @brief Returns the element-wise real parts in row-major order. + * @return Real parts of all entries. + */ + [[nodiscard]] std::array realPart() const; + + /** + * @brief Returns the element-wise imaginary parts in row-major order. + * @return Imaginary parts of all entries. + */ + [[nodiscard]] std::array imagPart() const; + /** * @brief Checks approximate equality using an absolute entry-wise tolerance. * @@ -539,6 +605,41 @@ class DynamicMatrix { [[nodiscard]] bool isApprox(const DynamicMatrix& other, double tol = MATRIX_TOLERANCE) const; + /** + * @brief Returns the trace of this matrix. + * @return Sum of diagonal entries. + */ + [[nodiscard]] Complex trace() const; + + /** + * @brief Matrix product `*this * rhs`. + * @param rhs Right-hand factor. + * @return Product of the two matrices. + */ + [[nodiscard]] DynamicMatrix operator*(const DynamicMatrix& rhs) const; + + /** + * @brief Element-wise scaling by a complex scalar. + * @param scalar Factor applied to every matrix entry. + * @return Scaled copy of this matrix. + */ + [[nodiscard]] DynamicMatrix operator*(const Complex& scalar) const; + + /** + * @brief Element-wise in-place scaling by a complex scalar. + * @param scalar Factor applied to every matrix entry. + * @return Reference to this matrix. + */ + DynamicMatrix& operator*=(const Complex& scalar); + + /** + * @brief Checks whether this matrix is approximately the identity. + * @param tol Maximum allowed complex modulus of each off-diagonal entry and + * each diagonal deviation from one. + * @return True when the matrix is close to the identity. + */ + [[nodiscard]] bool isIdentity(double tol = MATRIX_TOLERANCE) const; + private: struct Impl; std::unique_ptr impl_; @@ -558,4 +659,99 @@ inline constexpr bool std::disjunction_v, std::is_same, std::is_same, std::is_same>; + +/** + * @brief Kronecker product `lhs (x) rhs` of two single-qubit matrices. + * + * Uses the computational-basis bit order where the first operand labels the + * high bit, matching `UnitaryOpInterface::getUnitaryMatrix4x4`. + * + * @param lhs Left factor (acts on the high bit / qubit 0). + * @param rhs Right factor (acts on the low bit / qubit 1). + * @return The `4x4` Kronecker product. + */ +[[nodiscard]] Matrix4x4 kron(const Matrix2x2& lhs, const Matrix2x2& rhs); + +/// Scalar-on-the-left multiply `scalar * matrix` (commutes with the member +/// `matrix * scalar`). Provided so generic code can scale a matrix from +/// either side. +[[nodiscard]] Matrix2x2 operator*(const Complex& scalar, + const Matrix2x2& matrix); +/// @copydoc operator*(const Complex&, const Matrix2x2&) +[[nodiscard]] Matrix4x4 operator*(const Complex& scalar, + const Matrix4x4& matrix); +/// @copydoc operator*(const Complex&, const Matrix2x2&) +[[nodiscard]] DynamicMatrix operator*(const Complex& scalar, + const DynamicMatrix& matrix); + +/** + * @brief Eigenvalues and eigenvectors of a real symmetric `4x4` matrix. + * + * `eigenvalues` are sorted ascending and `eigenvectors` holds the + * corresponding orthonormal eigenvectors as columns (column `j` is the + * eigenvector for `eigenvalues[j]`). + */ +struct SymmetricEigen4 { + std::array eigenvalues{}; + Matrix4x4 eigenvectors{}; +}; + +/** + * @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. + * + * @pre @p symmetric is real and symmetric: `symmetric[i,j] == symmetric[j,i]` + * for all `i, j`. Only the lower triangle (including the diagonal) is read, + * but supplying a non-symmetric matrix yields undefined numerical results. + * + * @param symmetric Row-major real symmetric `4x4` matrix. + * @return Ascending eigenvalues and matching eigenvectors (as columns). + */ +[[nodiscard]] SymmetricEigen4 +symmetricEigen4(const std::array& symmetric); + +/** + * @brief Computes the eigendecomposition of a real symmetric `4x4` matrix. + * + * @copydoc symmetricEigen4(const std::array&) + * + * @pre Entries of @p symmetric are real (imaginary parts must be negligible). + * The real parts must form a symmetric matrix; imaginary parts are ignored. + */ +[[nodiscard]] SymmetricEigen4 symmetricEigen4(const Matrix4x4& symmetric); + +/** + * @brief Reorder a two-qubit matrix to act on qubits `{0, 1}`. + * + * @param q0Index Wire index of operand 0; @p q1Index wire index of operand 1. + */ +[[nodiscard]] Matrix4x4 reorderTwoQubitMatrix(const Matrix4x4& matrix, + std::size_t q0Index, + std::size_t q1Index); + +/** + * @brief Embed a single-qubit matrix into an @p numQubits-qubit Hilbert space. + * + * Qubit @p qubitIndex uses the same MSB-first convention as @ref kron + * (high bit first operand, low bit second). For each basis pair whose untouched + * wires match, copies @p matrix at the target qubit's row/column bits. + */ +[[nodiscard]] DynamicMatrix embedSingleQubitInNqubit(const Matrix2x2& matrix, + std::size_t numQubits, + std::size_t qubitIndex); + +/** + * @brief Embed a two-qubit matrix into an @p numQubits-qubit Hilbert space. + * + * Operand 0 labels the high bit of the pair and acts on @p q0Index; operand 1 + * labels the low bit and acts on @p q1Index. For each basis pair whose other + * wires match, copies @p matrix at the packed two-qubit row/column indices. + */ +[[nodiscard]] DynamicMatrix embedTwoQubitInNqubit(const Matrix4x4& matrix, + std::size_t numQubits, + std::size_t q0Index, + std::size_t q1Index); + } // namespace mlir::qco diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 8df45840a3..7a53347c04 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -11,6 +11,7 @@ #include "mlir/Dialect/QCO/Utils/Matrix.h" #include +#include #include #include @@ -81,21 +82,32 @@ assignFromDynamicImpl(const DynamicMatrix& src, } /// Writes the row-major product `lhs * rhs` into @p out (2x2, fully unrolled). -static void -multiply2x2(const std::array& lhs, - const std::array& rhs, - std::array& out) { +static void multiply2x2(const ArrayRef lhs, + const ArrayRef rhs, + const MutableArrayRef out) { + assert(lhs.size() == Matrix2x2::K_SIZE_AT_COMPILE_TIME && + rhs.size() == Matrix2x2::K_SIZE_AT_COMPILE_TIME && + out.size() == Matrix2x2::K_SIZE_AT_COMPILE_TIME); out[0] = lhs[0] * rhs[0] + lhs[1] * rhs[2]; out[1] = lhs[0] * rhs[1] + lhs[1] * rhs[3]; out[2] = lhs[2] * rhs[0] + lhs[3] * rhs[2]; out[3] = lhs[2] * rhs[1] + lhs[3] * rhs[3]; } -/// Writes the row-major product `lhs * rhs` into @p out (4x4, unrolled rows). static void -multiply4x4(const std::array& lhs, - const std::array& rhs, - std::array& out) { +multiply2x2(const std::array& lhs, + const std::array& rhs, + std::array& out) { + multiply2x2(ArrayRef(lhs), ArrayRef(rhs), MutableArrayRef(out)); +} + +/// Writes the row-major product `lhs * rhs` into @p out (4x4, unrolled rows). +static void multiply4x4(const ArrayRef lhs, + const ArrayRef rhs, + const MutableArrayRef out) { + assert(lhs.size() == Matrix4x4::K_SIZE_AT_COMPILE_TIME && + rhs.size() == Matrix4x4::K_SIZE_AT_COMPILE_TIME && + out.size() == Matrix4x4::K_SIZE_AT_COMPILE_TIME); for (std::size_t row = 0; row < Matrix4x4::K_ROWS; ++row) { const std::size_t rowBase = row * Matrix4x4::K_COLS; const Complex& a0 = lhs[rowBase + 0]; @@ -109,6 +121,34 @@ multiply4x4(const std::array& lhs, } } +static void +multiply4x4(const std::array& lhs, + const std::array& rhs, + std::array& out) { + multiply4x4(ArrayRef(lhs), ArrayRef(rhs), MutableArrayRef(out)); +} + +/// Returns true if @p data is approximately the @p dim x @p dim identity +/// matrix. +template +[[nodiscard]] static bool +isIdentityEntries(const std::array& data, + const double tol) { + for (std::size_t row = 0; row < Dim; ++row) { + for (std::size_t col = 0; col < Dim; ++col) { + const Complex& entry = data[(row * Dim) + col]; + if (row == col) { + if (std::abs(entry - Complex{1.0, 0.0}) > tol) { + return false; + } + } else if (std::abs(entry) > tol) { + return false; + } + } + } + return true; +} + /// Returns @p dim as `size_t` after asserting it is non-negative and squarable. [[nodiscard]] static std::size_t checkedDim(const std::int64_t dim) { assert(dim >= 0 && "DynamicMatrix dimension must be non-negative"); @@ -231,12 +271,20 @@ Matrix2x2 Matrix2x2::adjoint() const { std::conj(data[1]), std::conj(data[3])); } +Matrix2x2 Matrix2x2::transpose() const { + return fromElements(data[0], data[2], data[1], data[3]); +} + Complex Matrix2x2::trace() const { return data[0] + data[3]; } Complex Matrix2x2::determinant() const { return data[0] * data[3] - data[1] * data[2]; } +bool Matrix2x2::isIdentity(const double tol) const { + return isIdentityEntries(data, tol); +} + bool Matrix2x2::isApprox(const Matrix2x2& other, const double tol) const { return entriesAreApprox(data, other.data, tol); } @@ -296,6 +344,12 @@ Matrix4x4 Matrix4x4::adjoint() const { return out; } +Matrix4x4 Matrix4x4::transpose() const { + return fromElements(data[0], data[4], data[8], data[12], data[1], data[5], + data[9], data[13], data[2], data[6], data[10], data[14], + data[3], data[7], data[11], data[15]); +} + Complex Matrix4x4::trace() const { return data[0] + data[5] + data[10] + data[15]; } @@ -317,6 +371,52 @@ Complex Matrix4x4::determinant() const { data[12], data[13], data[14]); } +bool Matrix4x4::isIdentity(const double tol) const { + return isIdentityEntries(data, tol); +} + +std::array Matrix4x4::diagonal() const { + return {data[0], data[5], data[10], data[15]}; +} + +Matrix4x4 +Matrix4x4::fromDiagonal(const std::array& diagonalEntries) { + Matrix4x4 out{}; + for (std::size_t i = 0; i < K_ROWS; ++i) { + out.data[(i * K_COLS) + i] = diagonalEntries[i]; + } + return out; +} + +std::array +Matrix4x4::column(const std::size_t col) const { + return {data[col], data[K_COLS + col], data[(2 * K_COLS) + col], + data[(3 * K_COLS) + col]}; +} + +void Matrix4x4::setColumn(const std::size_t col, + const std::array& values) { + for (std::size_t row = 0; row < K_ROWS; ++row) { + data[(row * K_COLS) + col] = values[row]; + } +} + +std::array +Matrix4x4::realPart() const { + return {data[0].real(), data[1].real(), data[2].real(), data[3].real(), + data[4].real(), data[5].real(), data[6].real(), data[7].real(), + data[8].real(), data[9].real(), data[10].real(), data[11].real(), + data[12].real(), data[13].real(), data[14].real(), data[15].real()}; +} + +std::array +Matrix4x4::imagPart() const { + return {data[0].imag(), data[1].imag(), data[2].imag(), data[3].imag(), + data[4].imag(), data[5].imag(), data[6].imag(), data[7].imag(), + data[8].imag(), data[9].imag(), data[10].imag(), data[11].imag(), + data[12].imag(), data[13].imag(), data[14].imag(), data[15].imag()}; +} + bool Matrix4x4::isApprox(const Matrix4x4& other, const double tol) const { return entriesAreApprox(data, other.data, tol); } @@ -453,4 +553,493 @@ bool DynamicMatrix::isApprox(const DynamicMatrix& other, return entriesAreApprox(impl_->data, other.impl_->data, tol); } +Complex DynamicMatrix::trace() const { + Complex sum{0.0, 0.0}; + const auto udim = checkedDim(impl_->dim); + for (std::size_t i = 0; i < udim; ++i) { + sum += impl_->data[(i * udim) + i]; + } + return sum; +} + +DynamicMatrix DynamicMatrix::operator*(const DynamicMatrix& rhs) const { + if (impl_->dim != rhs.impl_->dim) { + llvm::reportFatalInternalError( + "DynamicMatrix multiply requires matching dimensions"); + } + DynamicMatrix out(impl_->dim); + if (impl_->dim == static_cast(Matrix2x2::K_ROWS)) { + multiply2x2(impl_->data, rhs.impl_->data, out.impl_->data); + return out; + } + if (impl_->dim == static_cast(Matrix4x4::K_ROWS)) { + multiply4x4(impl_->data, rhs.impl_->data, out.impl_->data); + return out; + } + + const auto udim = checkedDim(impl_->dim); + for (std::size_t row = 0; row < udim; ++row) { + for (std::size_t col = 0; col < udim; ++col) { + Complex sum{0.0, 0.0}; + for (std::size_t k = 0; k < udim; ++k) { + sum += + impl_->data[(row * udim) + k] * rhs.impl_->data[(k * udim) + col]; + } + out.impl_->data[(row * udim) + col] = sum; + } + } + return out; +} + +DynamicMatrix DynamicMatrix::operator*(const Complex& scalar) const { + DynamicMatrix out(impl_->dim); + for (std::size_t i = 0; i < impl_->data.size(); ++i) { + out.impl_->data[i] = impl_->data[i] * scalar; + } + return out; +} + +DynamicMatrix& DynamicMatrix::operator*=(const Complex& scalar) { + for (Complex& entry : impl_->data) { + entry *= scalar; + } + return *this; +} + +bool DynamicMatrix::isIdentity(const double tol) const { + const auto udim = checkedDim(impl_->dim); + for (std::size_t row = 0; row < udim; ++row) { + for (std::size_t col = 0; col < udim; ++col) { + const Complex& entry = impl_->data[(row * udim) + col]; + if (row == col) { + if (std::abs(entry - Complex{1.0, 0.0}) > tol) { + return false; + } + } else if (std::abs(entry) > tol) { + return false; + } + } + } + return true; +} + +/** + * @brief Returns the @p qubitIndex bit of a computational-basis label. + * + * Qubit 0 is the MSB of @p stateIndex, matching @ref kron and + * @ref embedSingleQubitInNqubit. + */ +[[nodiscard]] static std::size_t qubitBitAt(const std::size_t stateIndex, + const std::size_t numQubits, + const std::size_t qubitIndex) { + return (stateIndex >> (numQubits - 1 - qubitIndex)) & 1U; +} + +/** + * @brief True when row and col agree on every wire except @p skipA and @p + * skipB. + * + * Used when embedding a gate: untouched qubits must match or the matrix entry + * is zero. For a single-qubit embed, pass @p skipB = @p numQubits so only @p + * skipA is skipped. + */ +[[nodiscard]] static bool otherQubitBitsMatch(const std::size_t row, + const std::size_t col, + const std::size_t numQubits, + const std::size_t skipA, + const std::size_t skipB) { + for (std::size_t q = 0; q < numQubits; ++q) { + if (q == skipA || q == skipB) { + continue; + } + if (qubitBitAt(row, numQubits, q) != qubitBitAt(col, numQubits, q)) { + return false; + } + } + return true; +} + +Matrix4x4 kron(const Matrix2x2& lhs, const Matrix2x2& rhs) { + const auto& a = lhs.data; + const auto& b = rhs.data; + return Matrix4x4::fromElements( + a[0] * b[0], a[0] * b[1], a[1] * b[0], a[1] * b[1], a[0] * b[2], + a[0] * b[3], a[1] * b[2], a[1] * b[3], a[2] * b[0], a[2] * b[1], + a[3] * b[0], a[3] * b[1], a[2] * b[2], a[2] * b[3], a[3] * b[2], + a[3] * b[3]); +} + +Matrix2x2 operator*(const Complex& scalar, const Matrix2x2& matrix) { + return matrix * scalar; +} + +Matrix4x4 operator*(const Complex& scalar, const Matrix4x4& matrix) { + return matrix * scalar; +} + +DynamicMatrix operator*(const Complex& scalar, const DynamicMatrix& matrix) { + return matrix * scalar; +} + +[[nodiscard]] static double pythag2(const double a, const double b) { + const double absA = std::abs(a); + const double absB = std::abs(b); + if (absA > absB) { + const double ratio = absB / absA; + return absA * std::sqrt(1.0 + (ratio * ratio)); + } + if (absB == 0.0) { + return 0.0; + } + const double ratio = absA / absB; + return absB * std::sqrt(1.0 + (ratio * ratio)); +} + +[[nodiscard]] static double copySign(const double magnitude, + const double signSource) { + return (signSource >= 0.0) ? std::abs(magnitude) : -std::abs(magnitude); +} + +// Adapted from John Burkardt's MIT-licensed EISPACK C port (`tred2` / `tql2`): +// https://people.sc.fsu.edu/~jburkardt/c_src/eispack/eispack.c +// Specialized to `n = 4`; input is row-major, accumulator `z` is column-major. + +/// EISPACK `tred2` for `n = 4` (column-major `z[row + col*n]`). +static void symmetricTred24(const std::array& input, + std::array& z, + std::array& diag, + std::array& subdiag) { + constexpr std::size_t n = 4; + const auto zAt = [&z](const std::size_t row, + const std::size_t col) -> double& { + return z[row + (col * n)]; + }; + double h = 0.0; + + for (std::size_t col = 0; col < n; ++col) { + for (std::size_t row = col; row < n; ++row) { + zAt(row, col) = input[(row * n) + col]; + } + diag[col] = input[((n - 1) * n) + col]; + } + + for (int i = static_cast(n) - 1; i >= 1; --i) { + const auto ui = static_cast(i); + const std::size_t l = ui - 1; + h = 0.0; + double scale = 0.0; + for (std::size_t k = 0; k <= l; ++k) { + scale += std::abs(diag[k]); + } + if (scale == 0.0) { + subdiag[ui] = diag[l]; + for (std::size_t j = 0; j <= l; ++j) { + diag[j] = zAt(l, j); + zAt(ui, j) = 0.0; + zAt(j, ui) = 0.0; + } + diag[ui] = 0.0; + continue; + } + for (std::size_t k = 0; k <= l; ++k) { + diag[k] /= scale; + } + for (std::size_t k = 0; k <= l; ++k) { + h += diag[k] * diag[k]; + } + const double f = diag[l]; + const double g = -std::sqrt(h) * copySign(1.0, f); + subdiag[ui] = scale * g; + h -= f * g; + diag[l] = f - g; + + for (std::size_t k = 0; k <= l; ++k) { + subdiag[k] = 0.0; + } + for (std::size_t j = 0; j <= l; ++j) { + const double fj = diag[j]; + zAt(j, ui) = fj; + double gj = subdiag[j] + (zAt(j, j) * fj); + for (std::size_t k = j + 1; k <= l; ++k) { + gj += zAt(k, j) * diag[k]; + subdiag[k] += zAt(k, j) * fj; + } + subdiag[j] = gj; + } + double ff = 0.0; + for (std::size_t k = 0; k <= l; ++k) { + subdiag[k] /= h; + ff += subdiag[k] * diag[k]; + } + const double hh = 0.5 * ff / h; + for (std::size_t k = 0; k <= l; ++k) { + subdiag[k] -= hh * diag[k]; + } + for (std::size_t j = 0; j <= l; ++j) { + const double fj = diag[j]; + const double gj = subdiag[j]; + for (std::size_t k = j; k <= l; ++k) { + zAt(k, j) -= (fj * subdiag[k]) + (gj * diag[k]); + } + diag[j] = zAt(l, j); + zAt(ui, j) = 0.0; + } + diag[ui] = h; + } + + for (std::size_t i = 1; i < n; ++i) { + const std::size_t l = i - 1; + zAt(n - 1, l) = zAt(l, l); + zAt(l, l) = 1.0; + h = diag[i]; + if (h != 0.0) { + for (std::size_t k = 0; k <= l; ++k) { + diag[k] = zAt(k, i) / h; + } + for (std::size_t j = 0; j <= l; ++j) { + double g = 0.0; + for (std::size_t k = 0; k <= l; ++k) { + g += zAt(k, i) * zAt(k, j); + } + for (std::size_t k = 0; k <= l; ++k) { + zAt(k, j) -= g * diag[k]; + } + } + } + for (std::size_t k = 0; k <= l; ++k) { + zAt(k, i) = 0.0; + } + } + + for (std::size_t j = 0; j < n; ++j) { + diag[j] = zAt(n - 1, j); + } + for (std::size_t j = 0; j < n - 1; ++j) { + zAt(n - 1, j) = 0.0; + } + zAt(n - 1, n - 1) = 1.0; + subdiag[0] = 0.0; +} + +/// EISPACK `tql2` for `n = 4` (column-major `z[row + col*n]`). +static void symmetricTql24(std::array& diag, + std::array& subdiag, + std::array& z) { + constexpr std::size_t n = 4; + const auto zAt = [&z](const std::size_t row, + const std::size_t col) -> double& { + return z[row + (col * n)]; + }; + + for (std::size_t i = 1; i < n; ++i) { + subdiag[i - 1] = subdiag[i]; + } + double f = 0.0; + double tst1 = 0.0; + subdiag[n - 1] = 0.0; + + for (std::size_t l = 0; l < n; ++l) { + int j = 0; + const double h = std::abs(diag[l]) + std::abs(subdiag[l]); + tst1 = std::max(tst1, h); + + std::size_t m = l; + for (; m < n; ++m) { + const double tst2 = tst1 + std::abs(subdiag[m]); + if (tst2 == tst1) { + break; + } + } + + if (m != l) { + while (true) { + if (j == 30) { + llvm::reportFatalInternalError("symmetricTql2_4: failed to converge"); + } + ++j; + + const std::size_t l1 = l + 1; + const std::size_t l2 = l1 + 1; + const double g = diag[l]; + const double p = (diag[l1] - g) / (2.0 * subdiag[l]); + const double r = pythag2(p, 1.0); + diag[l] = subdiag[l] / (p + copySign(std::abs(r), p)); + diag[l1] = subdiag[l] * (p + copySign(std::abs(r), p)); + const double dl1 = diag[l1]; + const double hh = g - diag[l]; + for (std::size_t i = l2; i < n; ++i) { + diag[i] -= hh; + } + f += hh; + + double pv = diag[m]; + double c = 1.0; + double c2 = c; + const double el1 = subdiag[l1]; + double s = 0.0; + double c3 = 1.0; + double s2 = 0.0; + const std::size_t mml = m - l; + for (std::size_t ii = 1; ii <= mml; ++ii) { + c3 = c2; + c2 = c; + s2 = s; + const std::size_t i = m - ii; + const double gi = c * subdiag[i]; + const double hi = c * pv; + const double ri = pythag2(pv, subdiag[i]); + subdiag[i + 1] = s * ri; + s = subdiag[i] / ri; + c = pv / ri; + pv = (c * diag[i]) - (s * gi); + diag[i + 1] = hi + (s * ((c * gi) + (s * diag[i]))); + for (std::size_t k = 0; k < n; ++k) { + const double zkI1 = zAt(k, i + 1); + zAt(k, i + 1) = (s * zAt(k, i)) + (c * zkI1); + zAt(k, i) = (c * zAt(k, i)) - (s * zkI1); + } + } + pv = -s * s2 * c3 * el1 * subdiag[l] / dl1; + subdiag[l] = s * pv; + diag[l] = c * pv; + const double tst2 = tst1 + std::abs(subdiag[l]); + if (tst2 > tst1) { + continue; + } + break; + } + } + diag[l] += f; + } + + for (std::size_t ii = 1; ii < n; ++ii) { + const std::size_t i = ii - 1; + std::size_t k = i; + double p = diag[i]; + for (std::size_t j = ii; j < n; ++j) { + if (diag[j] < p) { + k = j; + p = diag[j]; + } + } + if (k == i) { + continue; + } + diag[k] = diag[i]; + diag[i] = p; + for (std::size_t j = 0; j < n; ++j) { + const double tmp = zAt(j, i); + zAt(j, i) = zAt(j, k); + zAt(j, k) = tmp; + } + } +} + +SymmetricEigen4 symmetricEigen4(const std::array& symmetric) { + constexpr std::size_t n = 4; + + std::array z{}; + std::array diag{}; + std::array subdiag{}; + symmetricTred24(symmetric, z, diag, subdiag); + symmetricTql24(diag, subdiag, z); + + SymmetricEigen4 result; + for (std::size_t col = 0; col < n; ++col) { + result.eigenvalues[col] = diag[col]; + for (std::size_t row = 0; row < n; ++row) { + result.eigenvectors(row, col) = Complex{z[row + (col * n)], 0.0}; + } + } + return result; +} + +SymmetricEigen4 symmetricEigen4(const Matrix4x4& symmetric) { + return symmetricEigen4(symmetric.realPart()); +} + +static Matrix4x4 embedSingleQubitInTwoQubit(const Matrix2x2& matrix, + const std::size_t qubitIndex) { + if (qubitIndex == 0) { + return kron(matrix, Matrix2x2::identity()); + } + if (qubitIndex == 1) { + return kron(Matrix2x2::identity(), matrix); + } + llvm::reportFatalInternalError("Invalid qubit index for single-qubit embed"); +} + +Matrix4x4 reorderTwoQubitMatrix(const Matrix4x4& matrix, + const std::size_t q0Index, + const std::size_t q1Index) { + if (q0Index == 0 && q1Index == 1) { + return matrix; + } + if (q0Index == 1 && q1Index == 0) { + // Conjugate by SWAP: out[i, j] = matrix[pi(i), pi(j)] with pi swapping |01> + // and |10> (basis indices 1 and 2). + const auto& m = matrix.data; + return Matrix4x4::fromElements(m[0], m[2], m[1], m[3], m[8], m[10], m[9], + m[11], m[4], m[6], m[5], m[7], m[12], m[14], + m[13], m[15]); + } + llvm::reportFatalInternalError("Invalid qubit indices for two-qubit reorder"); +} + +DynamicMatrix embedSingleQubitInNqubit(const Matrix2x2& matrix, + const std::size_t numQubits, + const std::size_t qubitIndex) { + if (qubitIndex >= numQubits) { + llvm::reportFatalInternalError( + "Invalid qubit index for single-qubit embed"); + } + if (numQubits == 2) { + return DynamicMatrix(embedSingleQubitInTwoQubit(matrix, qubitIndex)); + } + const auto dim = static_cast(1ULL << numQubits); + DynamicMatrix out(dim); + const auto udim = static_cast(dim); + for (std::size_t row = 0; row < udim; ++row) { + for (std::size_t col = 0; col < udim; ++col) { + if (!otherQubitBitsMatch(row, col, numQubits, qubitIndex, numQubits)) { + continue; + } + const std::size_t rowBit = qubitBitAt(row, numQubits, qubitIndex); + const std::size_t colBit = qubitBitAt(col, numQubits, qubitIndex); + out(static_cast(row), static_cast(col)) = + matrix(rowBit, colBit); + } + } + return out; +} + +DynamicMatrix embedTwoQubitInNqubit(const Matrix4x4& matrix, + const std::size_t numQubits, + const std::size_t q0Index, + const std::size_t q1Index) { + if (q0Index >= numQubits || q1Index >= numQubits || q0Index == q1Index) { + llvm::reportFatalInternalError("Invalid qubit indices for two-qubit embed"); + } + if (numQubits == 2) { + return DynamicMatrix(reorderTwoQubitMatrix(matrix, q0Index, q1Index)); + } + const auto dim = static_cast(1ULL << numQubits); + DynamicMatrix out(dim); + const auto udim = static_cast(dim); + for (std::size_t row = 0; row < udim; ++row) { + for (std::size_t col = 0; col < udim; ++col) { + if (!otherQubitBitsMatch(row, col, numQubits, q0Index, q1Index)) { + continue; + } + const std::size_t rowPair = (qubitBitAt(row, numQubits, q0Index) << 1) | + qubitBitAt(row, numQubits, q1Index); + const std::size_t colPair = (qubitBitAt(col, numQubits, q0Index) << 1) | + qubitBitAt(col, numQubits, q1Index); + out(static_cast(row), static_cast(col)) = + matrix(rowPair, colPair); + } + } + return out; +} + } // namespace mlir::qco diff --git a/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp b/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp index afa0792415..22ce0315f7 100644 --- a/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp +++ b/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp @@ -12,8 +12,11 @@ #include +#include #include #include +#include +#include #include using namespace mlir::qco; @@ -130,6 +133,7 @@ TEST(UnitaryMatrix4x4, MultiplyAdjointTraceDeterminant) { EXPECT_TRUE(scaled.isApprox(swap * 2.0)); EXPECT_EQ(identity.trace(), Complex(4.0, 0.0)); EXPECT_EQ(identity.determinant(), Complex(1.0, 0.0)); + EXPECT_EQ(swap.determinant(), Complex(-1.0, 0.0)); } TEST(UnitaryMatrix4x4, PremultiplyBy) { @@ -186,6 +190,7 @@ TEST(DynamicMatrix, IdentityAndElementAccess) { EXPECT_EQ(identity(0, 0), 1.0); EXPECT_EQ(identity(1, 2), 0.0); EXPECT_EQ(identity(2, 2), 1.0); + EXPECT_TRUE(identity.isIdentity()); DynamicMatrix mutableMatrix = identity; mutableMatrix(1, 1) = 0.5; @@ -202,6 +207,13 @@ TEST(DynamicMatrix, FromAdjoint) { EXPECT_TRUE(DynamicMatrix(swapMatrix()).isApprox(swapMatrix())); } +TEST(DynamicMatrix, ScalarMultiplyAssign) { + DynamicMatrix matrix = DynamicMatrix::identity(2); + matrix *= std::exp(Complex{0.0, 0.5}); + EXPECT_TRUE(matrix.isApprox(DynamicMatrix::identity(2) * + std::exp(Complex{0.0, 0.5}))); +} + TEST(DynamicMatrix, AssignFrom) { DynamicMatrix dynamic; @@ -267,6 +279,38 @@ TEST(DynamicMatrix, IsApproxRejectsMismatchedExtents) { EXPECT_FALSE(DynamicMatrix::identity(1).isApprox(DynamicMatrix::identity(2))); } +TEST(DynamicMatrix, IsApproxOverloads) { + const Matrix1x1 phase = Matrix1x1::fromElements(Complex{0.25, 0.5}); + const Matrix2x2 x = pauliX(); + const Matrix4x4 swap = swapMatrix(); + + DynamicMatrix as1x1; + as1x1.assignFrom(phase); + EXPECT_TRUE(as1x1.isApprox(phase)); + EXPECT_FALSE(as1x1.isApprox(Matrix1x1::fromElements(1.0))); + + DynamicMatrix as2x2; + as2x2.assignFrom(x); + EXPECT_TRUE(as2x2.isApprox(x)); + EXPECT_FALSE(as2x2.isApprox(Matrix2x2::identity())); + + DynamicMatrix as4x4; + as4x4.assignFrom(swap); + EXPECT_TRUE(as4x4.isApprox(swap)); + EXPECT_FALSE(as4x4.isApprox(Matrix4x4::identity())); + + DynamicMatrix wrongDim = DynamicMatrix::identity(3); + EXPECT_FALSE(wrongDim.isApprox(phase)); + EXPECT_FALSE(wrongDim.isApprox(x)); + EXPECT_FALSE(wrongDim.isApprox(swap)); + + const DynamicMatrix a = DynamicMatrix::identity(2); + DynamicMatrix b = a; + b(1, 0) += 1e-15; + EXPECT_TRUE(a.isApprox(b)); + EXPECT_FALSE(a.isApprox(DynamicMatrix::identity(3))); +} + TEST(Matrix1x1, AssignFromDynamicMatrix) { const Matrix1x1 phase = Matrix1x1::fromElements(Complex{0.25, 0.5}); @@ -303,34 +347,219 @@ TEST(Matrix4x4, AssignFromDynamicMatrix) { EXPECT_FALSE(out.assignFrom(DynamicMatrix::identity(2))); } -TEST(DynamicMatrix, IsApproxOverloads) { - const Matrix1x1 phase = Matrix1x1::fromElements(Complex{0.25, 0.5}); +TEST(UnitaryMatrix2x2, TransposeAndIsIdentity) { + const Matrix2x2 m = Matrix2x2::fromElements(1, 2i, 3, 4); + EXPECT_TRUE(m.transpose().isApprox(Matrix2x2::fromElements(1, 3, 2i, 4))); + EXPECT_TRUE(Matrix2x2::identity().isIdentity()); + EXPECT_FALSE(pauliX().isIdentity()); +} + +TEST(UnitaryMatrix4x4, TransposeAndIsIdentity) { + Matrix4x4 m = Matrix4x4::identity(); + m(0, 3) = 2i; + m(3, 0) = 5.0; + const Matrix4x4 t = m.transpose(); + EXPECT_EQ(t(3, 0), 2i); + EXPECT_EQ(t(0, 3), 5.0); + EXPECT_TRUE(Matrix4x4::identity().isIdentity()); + EXPECT_FALSE(swapMatrix().isIdentity()); +} + +TEST(UnitaryMatrix4x4, DiagonalColumnsAndParts) { + Matrix4x4 m = + Matrix4x4::fromElements(Complex{1, 1}, 0, 0, 0, 0, Complex{2, 2}, 0, 0, 0, + 0, Complex{3, 3}, 0, 0, 0, 0, Complex{4, 4}); + const auto diag = m.diagonal(); + EXPECT_EQ(diag[0], (Complex{1, 1})); + EXPECT_EQ(diag[3], (Complex{4, 4})); + EXPECT_TRUE(Matrix4x4::fromDiagonal(diag).isApprox(m)); + + const auto col1 = m.column(1); + EXPECT_EQ(col1[1], (Complex{2, 2})); + Matrix4x4 n = Matrix4x4::identity(); + n.setColumn(2, {1i, 2i, 3i, 4i}); + EXPECT_EQ(n(0, 2), 1i); + EXPECT_EQ(n(3, 2), 4i); + + const auto re = m.realPart(); + const auto im = m.imagPart(); + EXPECT_EQ(re[0], 1.0); + EXPECT_EQ(im[0], 1.0); + EXPECT_EQ(re[15], 4.0); + EXPECT_EQ(im[15], 4.0); +} + +TEST(UnitaryMatrix4x4, KroneckerProduct) { const Matrix2x2 x = pauliX(); - const Matrix4x4 swap = swapMatrix(); + // X (x) I should swap the high bit. + const Matrix4x4 xi = kron(x, Matrix2x2::identity()); + EXPECT_TRUE(xi.isApprox(Matrix4x4::fromElements(0, 0, 1, 0, // row 0 + 0, 0, 0, 1, // row 1 + 1, 0, 0, 0, // row 2 + 0, 1, 0, 0))); + // I (x) X swaps the low bit. + const Matrix4x4 ix = kron(Matrix2x2::identity(), x); + EXPECT_TRUE(ix.isApprox(Matrix4x4::fromElements(0, 1, 0, 0, // row 0 + 1, 0, 0, 0, // row 1 + 0, 0, 0, 1, // row 2 + 0, 0, 1, 0))); +} - DynamicMatrix as1x1; - as1x1.assignFrom(phase); - EXPECT_TRUE(as1x1.isApprox(phase)); - EXPECT_FALSE(as1x1.isApprox(Matrix1x1::fromElements(1.0))); +TEST(UnitaryMatrix4x4, ReorderTwoQubitMatrix) { + const Matrix2x2 x = pauliX(); + const Matrix4x4 onHigh = kron(x, Matrix2x2::identity()); + const Matrix4x4 onLow = kron(Matrix2x2::identity(), x); - DynamicMatrix as2x2; - as2x2.assignFrom(x); - EXPECT_TRUE(as2x2.isApprox(x)); - EXPECT_FALSE(as2x2.isApprox(Matrix2x2::identity())); + EXPECT_TRUE(reorderTwoQubitMatrix(onHigh, 0, 1).isApprox(onHigh)); + EXPECT_TRUE(reorderTwoQubitMatrix(onHigh, 1, 0).isApprox(onLow)); + EXPECT_TRUE(reorderTwoQubitMatrix(onLow, 1, 0).isApprox(onHigh)); +} - DynamicMatrix as4x4; - as4x4.assignFrom(swap); - EXPECT_TRUE(as4x4.isApprox(swap)); - EXPECT_FALSE(as4x4.isApprox(Matrix4x4::identity())); +TEST(UnitaryDynamicMatrix, NQubitEmbedMatchesTwoQubitSpecialization) { + const Matrix2x2 x = pauliX(); + const Matrix4x4 cx = Matrix4x4::fromElements(1, 0, 0, 0, // + 0, 1, 0, 0, // + 0, 0, 0, 1, // + 0, 0, 1, 0); + EXPECT_TRUE(embedSingleQubitInNqubit(x, 2, 0).isApprox( + kron(x, Matrix2x2::identity()))); + EXPECT_TRUE(embedSingleQubitInNqubit(x, 2, 1).isApprox( + kron(Matrix2x2::identity(), x))); + EXPECT_TRUE(embedTwoQubitInNqubit(cx, 2, 0, 1) + .isApprox(reorderTwoQubitMatrix(cx, 0, 1))); + const DynamicMatrix cxOn01 = embedTwoQubitInNqubit(cx, 3, 0, 1); + const DynamicMatrix cxOn12 = embedTwoQubitInNqubit(cx, 3, 1, 2); + EXPECT_EQ(cxOn01.rows(), 8); + EXPECT_EQ(cxOn12.rows(), 8); + EXPECT_FALSE(cxOn01.isApprox(cxOn12)); +} - DynamicMatrix wrongDim = DynamicMatrix::identity(3); - EXPECT_FALSE(wrongDim.isApprox(phase)); - EXPECT_FALSE(wrongDim.isApprox(x)); - EXPECT_FALSE(wrongDim.isApprox(swap)); +TEST(UnitaryDynamicMatrix, EmbedSingleQubitOnMiddleWire) { + const Matrix2x2 x = pauliX(); + const DynamicMatrix embedded = embedSingleQubitInNqubit(x, 3, 1); + EXPECT_EQ(embedded.rows(), 8); + EXPECT_FALSE(embedded.isIdentity()); - const DynamicMatrix a = DynamicMatrix::identity(2); - DynamicMatrix b = a; - b(1, 0) += 1e-15; - EXPECT_TRUE(a.isApprox(b)); - EXPECT_FALSE(a.isApprox(DynamicMatrix::identity(3))); + const DynamicMatrix product = embedded * embedded; + EXPECT_TRUE(product.isIdentity()); + EXPECT_NEAR(product.trace().real(), 8.0, MATRIX_TOLERANCE); +} + +TEST(UnitaryDynamicMatrix, MultiplyTraceAndScalar) { + const Matrix2x2 x = pauliX(); + const DynamicMatrix embedded = embedSingleQubitInNqubit(x, 2, 0); + EXPECT_FALSE(embedded.isIdentity()); + const Complex scalar = std::exp(1i * 0.3); + EXPECT_TRUE((scalar * embedded).isApprox(embedded * scalar)); + const DynamicMatrix product = embedded * embedded; + EXPECT_TRUE(product.isIdentity()); + EXPECT_NEAR(product.trace().real(), 4.0, MATRIX_TOLERANCE); +} + +TEST(DynamicMatrix, MultiplyAdjointTraceAt4) { + const auto swap = DynamicMatrix(swapMatrix()); + EXPECT_EQ(swap.rows(), 4); + + const DynamicMatrix product = swap * swap; + EXPECT_TRUE(product.isIdentity()); + EXPECT_NEAR(product.trace().real(), 4.0, MATRIX_TOLERANCE); + + const DynamicMatrix adjoint = swap.adjoint(); + EXPECT_TRUE(adjoint.isApprox(swapMatrix())); +} + +TEST(UnitaryMatrix2x2, ScalarLeftMultiply) { + const Matrix2x2 x = pauliX(); + const Complex scalar = std::exp(1i * 0.5); + EXPECT_TRUE((scalar * x).isApprox(x * scalar)); +} + +TEST(UnitaryMatrix4x4, ScalarLeftMultiply) { + const Matrix4x4 swap = swapMatrix(); + const Complex scalar = std::exp(1i * 0.25); + EXPECT_TRUE((scalar * swap).isApprox(swap * scalar)); +} + +TEST(SymmetricEigensolver, DiagonalMatrix) { + std::array a{}; + a[0] = 3.0; + a[5] = 1.0; + a[10] = 4.0; + a[15] = 2.0; + const SymmetricEigen4 result = symmetricEigen4(a); + EXPECT_NEAR(result.eigenvalues[0], 1.0, MATRIX_TOLERANCE); + EXPECT_NEAR(result.eigenvalues[1], 2.0, MATRIX_TOLERANCE); + EXPECT_NEAR(result.eigenvalues[2], 3.0, MATRIX_TOLERANCE); + EXPECT_NEAR(result.eigenvalues[3], 4.0, MATRIX_TOLERANCE); +} + +TEST(SymmetricEigensolver, Matrix4x4Overload) { + std::array a{}; + a[0] = 3.0; + a[5] = 1.0; + a[10] = 4.0; + a[15] = 2.0; + Matrix4x4 matrix{}; + for (std::size_t k = 0; k < 16; ++k) { + matrix(k / 4, k % 4) = a[k]; + } + const SymmetricEigen4 fromArray = symmetricEigen4(a); + const SymmetricEigen4 fromMatrix = symmetricEigen4(matrix); + for (std::size_t i = 0; i < 4; ++i) { + EXPECT_NEAR(fromMatrix.eigenvalues[i], fromArray.eigenvalues[i], + MATRIX_TOLERANCE); + } + EXPECT_TRUE(fromMatrix.eigenvectors.isApprox(fromArray.eigenvectors)); +} + +TEST(SymmetricEigensolver, ReconstructsRandomSymmetric) { + std::mt19937 rng(0xC0FFEE); + std::uniform_real_distribution dist(-2.0, 2.0); + for (int trial = 0; trial < 50; ++trial) { + std::array a{}; + for (std::size_t i = 0; i < 4; ++i) { + for (std::size_t j = i; j < 4; ++j) { + const double value = dist(rng); + a[(i * 4) + j] = value; + a[(j * 4) + i] = value; + } + } + const SymmetricEigen4 result = symmetricEigen4(a); + + // Eigenvalues are ascending. + for (std::size_t i = 0; i + 1 < 4; ++i) { + EXPECT_LE(result.eigenvalues[i], + result.eigenvalues[i + 1] + MATRIX_TOLERANCE); + } + + // Eigenvectors are orthonormal: V^T V == I. + const Matrix4x4& v = result.eigenvectors; + EXPECT_TRUE((v.transpose() * v).isIdentity()); + + // Reconstruction: V D V^T == A. + const Matrix4x4 d = + Matrix4x4::fromDiagonal({result.eigenvalues[0], result.eigenvalues[1], + result.eigenvalues[2], result.eigenvalues[3]}); + const Matrix4x4 reconstructed = v * d * v.transpose(); + Matrix4x4 original{}; + for (std::size_t k = 0; k < 16; ++k) { + original(k / 4, k % 4) = a[k]; + } + EXPECT_TRUE(reconstructed.isApprox(original)); + } +} + +TEST(SymmetricEigensolver, HandlesDegenerateSpectrum) { + // A scalar multiple of the identity: every vector is an eigenvector, but the + // returned basis must still be orthonormal. + std::array a{}; + for (std::size_t i = 0; i < 4; ++i) { + a[(i * 4) + i] = 2.5; + } + const SymmetricEigen4 result = symmetricEigen4(a); + for (const double value : result.eigenvalues) { + EXPECT_NEAR(value, 2.5, MATRIX_TOLERANCE); + } + const Matrix4x4& v = result.eigenvectors; + EXPECT_TRUE((v.transpose() * v).isIdentity()); } From f892bafeb8176361b3d685a82870afe2afb0c6b8 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 12:20:11 +0200 Subject: [PATCH 02/10] =?UTF-8?q?=F0=9F=9A=A8=20Fix=20linter=20warnings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 7a53347c04..32932d3d6a 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -568,11 +568,11 @@ DynamicMatrix DynamicMatrix::operator*(const DynamicMatrix& rhs) const { "DynamicMatrix multiply requires matching dimensions"); } DynamicMatrix out(impl_->dim); - if (impl_->dim == static_cast(Matrix2x2::K_ROWS)) { + if (std::cmp_equal(impl_->dim, Matrix2x2::K_ROWS)) { multiply2x2(impl_->data, rhs.impl_->data, out.impl_->data); return out; } - if (impl_->dim == static_cast(Matrix4x4::K_ROWS)) { + if (std::cmp_equal(impl_->dim, Matrix4x4::K_ROWS)) { multiply4x4(impl_->data, rhs.impl_->data, out.impl_->data); return out; } From 7e3fe1aa5905cafebcba2da06f5c28202f8d23c5 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 12:51:02 +0200 Subject: [PATCH 03/10] =?UTF-8?q?=E2=98=82=EF=B8=8F=20Increase=20coverage?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 103 ++++++------------ .../Dialect/QCO/Utils/test_unitary_matrix.cpp | 20 ++++ 2 files changed, 56 insertions(+), 67 deletions(-) diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 32932d3d6a..ed55c4a212 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -27,14 +27,20 @@ namespace mlir::qco { +/// Returns true if @p lhs and @p rhs differ by at most @p tol (complex +/// modulus). +[[nodiscard]] static bool entryIsApprox(const Complex& lhs, const Complex& rhs, + const double tol) { + return std::abs(lhs - rhs) <= tol; +} + /// Returns true if every entry pair differs by at most @p tol (complex /// modulus). [[nodiscard]] static bool entriesAreApprox(ArrayRef lhs, ArrayRef rhs, double tol) { - return std::ranges::equal(lhs, rhs, - [tol](const Complex& a, const Complex& b) { - return std::abs(a - b) <= tol; - }); + return llvm::equal(lhs, rhs, [tol](const Complex& a, const Complex& b) { + return entryIsApprox(a, b, tol); + }); } /// Writes the conjugate transpose of @p in into @p out (square, row-major). @@ -130,15 +136,15 @@ multiply4x4(const std::array& lhs, /// Returns true if @p data is approximately the @p dim x @p dim identity /// matrix. -template -[[nodiscard]] static bool -isIdentityEntries(const std::array& data, - const double tol) { - for (std::size_t row = 0; row < Dim; ++row) { - for (std::size_t col = 0; col < Dim; ++col) { - const Complex& entry = data[(row * Dim) + col]; +[[nodiscard]] static bool isIdentityEntries(ArrayRef data, + const std::size_t dim, + const double tol) { + assert(data.size() >= dim * dim); + for (std::size_t row = 0; row < dim; ++row) { + for (std::size_t col = 0; col < dim; ++col) { + const Complex& entry = data[(row * dim) + col]; if (row == col) { - if (std::abs(entry - Complex{1.0, 0.0}) > tol) { + if (!entryIsApprox(entry, Complex{1.0, 0.0}, tol)) { return false; } } else if (std::abs(entry) > tol) { @@ -206,7 +212,7 @@ Complex Matrix1x1::operator()(const std::size_t row, } bool Matrix1x1::isApprox(const Matrix1x1& other, const double tol) const { - return std::abs(value - other.value) <= tol; + return entryIsApprox(value, other.value, tol); } bool Matrix1x1::assignFrom(const DynamicMatrix& src) { @@ -282,7 +288,7 @@ Complex Matrix2x2::determinant() const { } bool Matrix2x2::isIdentity(const double tol) const { - return isIdentityEntries(data, tol); + return isIdentityEntries(data, K_ROWS, tol); } bool Matrix2x2::isApprox(const Matrix2x2& other, const double tol) const { @@ -372,7 +378,7 @@ Complex Matrix4x4::determinant() const { } bool Matrix4x4::isIdentity(const double tol) const { - return isIdentityEntries(data, tol); + return isIdentityEntries(data, K_ROWS, tol); } std::array Matrix4x4::diagonal() const { @@ -533,7 +539,7 @@ bool DynamicMatrix::isApprox(const Matrix1x1& other, const double tol) const { if (impl_->dim != 1) { return false; } - return std::abs(impl_->data[0] - other.value) <= tol; + return entryIsApprox(impl_->data[0], other.value, tol); } bool DynamicMatrix::isApprox(const Matrix2x2& other, const double tol) const { @@ -563,10 +569,8 @@ Complex DynamicMatrix::trace() const { } DynamicMatrix DynamicMatrix::operator*(const DynamicMatrix& rhs) const { - if (impl_->dim != rhs.impl_->dim) { - llvm::reportFatalInternalError( - "DynamicMatrix multiply requires matching dimensions"); - } + assert(impl_->dim == rhs.impl_->dim && + "DynamicMatrix multiply requires matching dimensions"); DynamicMatrix out(impl_->dim); if (std::cmp_equal(impl_->dim, Matrix2x2::K_ROWS)) { multiply2x2(impl_->data, rhs.impl_->data, out.impl_->data); @@ -607,20 +611,7 @@ DynamicMatrix& DynamicMatrix::operator*=(const Complex& scalar) { } bool DynamicMatrix::isIdentity(const double tol) const { - const auto udim = checkedDim(impl_->dim); - for (std::size_t row = 0; row < udim; ++row) { - for (std::size_t col = 0; col < udim; ++col) { - const Complex& entry = impl_->data[(row * udim) + col]; - if (row == col) { - if (std::abs(entry - Complex{1.0, 0.0}) > tol) { - return false; - } - } else if (std::abs(entry) > tol) { - return false; - } - } - } - return true; + return isIdentityEntries(impl_->data, checkedDim(impl_->dim), tol); } /** @@ -681,25 +672,6 @@ DynamicMatrix operator*(const Complex& scalar, const DynamicMatrix& matrix) { return matrix * scalar; } -[[nodiscard]] static double pythag2(const double a, const double b) { - const double absA = std::abs(a); - const double absB = std::abs(b); - if (absA > absB) { - const double ratio = absB / absA; - return absA * std::sqrt(1.0 + (ratio * ratio)); - } - if (absB == 0.0) { - return 0.0; - } - const double ratio = absA / absB; - return absB * std::sqrt(1.0 + (ratio * ratio)); -} - -[[nodiscard]] static double copySign(const double magnitude, - const double signSource) { - return (signSource >= 0.0) ? std::abs(magnitude) : -std::abs(magnitude); -} - // Adapted from John Burkardt's MIT-licensed EISPACK C port (`tred2` / `tql2`): // https://people.sc.fsu.edu/~jburkardt/c_src/eispack/eispack.c // Specialized to `n = 4`; input is row-major, accumulator `z` is column-major. @@ -748,7 +720,7 @@ static void symmetricTred24(const std::array& input, h += diag[k] * diag[k]; } const double f = diag[l]; - const double g = -std::sqrt(h) * copySign(1.0, f); + const double g = -std::sqrt(h) * std::copysign(1.0, f); subdiag[ui] = scale * g; h -= f * g; diag[l] = f - g; @@ -862,9 +834,9 @@ static void symmetricTql24(std::array& diag, const std::size_t l2 = l1 + 1; const double g = diag[l]; const double p = (diag[l1] - g) / (2.0 * subdiag[l]); - const double r = pythag2(p, 1.0); - diag[l] = subdiag[l] / (p + copySign(std::abs(r), p)); - diag[l1] = subdiag[l] * (p + copySign(std::abs(r), p)); + const double r = std::hypot(p, 1.0); + diag[l] = subdiag[l] / (p + std::copysign(std::abs(r), p)); + diag[l1] = subdiag[l] * (p + std::copysign(std::abs(r), p)); const double dl1 = diag[l1]; const double hh = g - diag[l]; for (std::size_t i = l2; i < n; ++i) { @@ -887,7 +859,7 @@ static void symmetricTql24(std::array& diag, const std::size_t i = m - ii; const double gi = c * subdiag[i]; const double hi = c * pv; - const double ri = pythag2(pv, subdiag[i]); + const double ri = std::hypot(pv, subdiag[i]); subdiag[i + 1] = s * ri; s = subdiag[i] / ri; c = pv / ri; @@ -966,7 +938,7 @@ static Matrix4x4 embedSingleQubitInTwoQubit(const Matrix2x2& matrix, if (qubitIndex == 1) { return kron(Matrix2x2::identity(), matrix); } - llvm::reportFatalInternalError("Invalid qubit index for single-qubit embed"); + assert(false && "Invalid qubit index for single-qubit embed"); } Matrix4x4 reorderTwoQubitMatrix(const Matrix4x4& matrix, @@ -983,16 +955,14 @@ Matrix4x4 reorderTwoQubitMatrix(const Matrix4x4& matrix, m[11], m[4], m[6], m[5], m[7], m[12], m[14], m[13], m[15]); } - llvm::reportFatalInternalError("Invalid qubit indices for two-qubit reorder"); + assert(false && "Invalid qubit indices for two-qubit reorder"); } DynamicMatrix embedSingleQubitInNqubit(const Matrix2x2& matrix, const std::size_t numQubits, const std::size_t qubitIndex) { - if (qubitIndex >= numQubits) { - llvm::reportFatalInternalError( - "Invalid qubit index for single-qubit embed"); - } + assert(qubitIndex < numQubits && + "Invalid qubit index for single-qubit embed"); if (numQubits == 2) { return DynamicMatrix(embedSingleQubitInTwoQubit(matrix, qubitIndex)); } @@ -1017,9 +987,8 @@ DynamicMatrix embedTwoQubitInNqubit(const Matrix4x4& matrix, const std::size_t numQubits, const std::size_t q0Index, const std::size_t q1Index) { - if (q0Index >= numQubits || q1Index >= numQubits || q0Index == q1Index) { - llvm::reportFatalInternalError("Invalid qubit indices for two-qubit embed"); - } + assert(q0Index < numQubits && q1Index < numQubits && q0Index != q1Index && + "Invalid qubit indices for two-qubit embed"); if (numQubits == 2) { return DynamicMatrix(reorderTwoQubitMatrix(matrix, q0Index, q1Index)); } diff --git a/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp b/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp index 22ce0315f7..946dc1956e 100644 --- a/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp +++ b/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp @@ -352,6 +352,11 @@ TEST(UnitaryMatrix2x2, TransposeAndIsIdentity) { EXPECT_TRUE(m.transpose().isApprox(Matrix2x2::fromElements(1, 3, 2i, 4))); EXPECT_TRUE(Matrix2x2::identity().isIdentity()); EXPECT_FALSE(pauliX().isIdentity()); + Matrix2x2 nearIdentity = Matrix2x2::identity(); + nearIdentity(0, 1) = 1e-15; + EXPECT_TRUE(nearIdentity.isIdentity()); + nearIdentity(0, 1) = 1.0; + EXPECT_FALSE(nearIdentity.isIdentity()); } TEST(UnitaryMatrix4x4, TransposeAndIsIdentity) { @@ -468,6 +473,21 @@ TEST(DynamicMatrix, MultiplyAdjointTraceAt4) { EXPECT_TRUE(adjoint.isApprox(swapMatrix())); } +TEST(DynamicMatrix, MultiplyAt2) { + const DynamicMatrix x(pauliX()); + EXPECT_EQ(x.rows(), 2); + const DynamicMatrix product = x * x; + EXPECT_TRUE(product.isIdentity()); + EXPECT_NEAR(product.trace().real(), 2.0, MATRIX_TOLERANCE); + EXPECT_TRUE(product.isApprox(pauliX() * pauliX())); +} + +TEST(DynamicMatrix, IsIdentityOffDiagonal) { + DynamicMatrix matrix = DynamicMatrix::identity(2); + matrix(0, 1) = 1.0; + EXPECT_FALSE(matrix.isIdentity()); +} + TEST(UnitaryMatrix2x2, ScalarLeftMultiply) { const Matrix2x2 x = pauliX(); const Complex scalar = std::exp(1i * 0.5); From a2484a18f23ba74f855496216738996e68d7d3d9 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 13:43:15 +0200 Subject: [PATCH 04/10] =?UTF-8?q?=E2=9A=A1=EF=B8=8F=20Suggestions=20from?= =?UTF-8?q?=20code=20review=20to=20improve=20performance?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/include/mlir/Dialect/QCO/Utils/Matrix.h | 187 ++-- mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 865 +++++++++--------- .../Dialect/QCO/Utils/test_unitary_matrix.cpp | 58 +- 3 files changed, 578 insertions(+), 532 deletions(-) diff --git a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h index 402e86ac4a..369cd2209b 100644 --- a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h +++ b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h @@ -10,6 +10,9 @@ #pragma once +#include +#include + #include #include #include @@ -26,6 +29,8 @@ using Complex = std::complex; inline constexpr double MATRIX_TOLERANCE = 1e-14; class DynamicMatrix; +struct Matrix4x4; +struct SymmetricEigen4; /** * @brief 1x1 matrix for global-phase gates. @@ -228,6 +233,30 @@ struct Matrix2x2 { * @return `true` when @p src is 2x2. */ [[nodiscard]] bool assignFrom(const DynamicMatrix& src); + + /** + * @brief Embed this single-qubit matrix into an @p numQubits-qubit Hilbert + * space. + * + * Wire @p qubitIndex uses the same MSB-first convention as @ref + * Matrix4x4::kron (high bit first operand, low bit second). For each basis + * pair whose untouched wires match, copies this matrix at the target qubit's + * row/column bits. + * + * @param numQubits Number of qubits in the target Hilbert space. + * @param qubitIndex Wire index to act on. + * @return Embedded unitary as a dynamic matrix. + */ + [[nodiscard]] DynamicMatrix embedInNqubit(std::size_t numQubits, + std::size_t qubitIndex) const; + + /** + * @brief Embed this single-qubit matrix into a two-qubit Hilbert space. + * + * @param qubitIndex Wire index (`0` = high bit / MSB, `1` = low bit). + * @return The `4x4` embedded unitary. + */ + [[nodiscard]] Matrix4x4 embedInTwoQubit(std::size_t qubitIndex) const; }; /** @@ -368,11 +397,25 @@ struct Matrix4x4 { /** * @brief Builds a diagonal matrix from four diagonal entries. - * @param diagonalEntries Diagonal entries `(m00, m11, m22, m33)`. + * @param diagonalEntries Diagonal entries `(m00, m11, m22, m33)`; must have + * length `K_ROWS`. * @return Diagonal matrix with the given entries. */ [[nodiscard]] static Matrix4x4 - fromDiagonal(const std::array& diagonalEntries); + fromDiagonal(ArrayRef diagonalEntries); + + /** + * @brief Kronecker product `lhs (x) rhs` of two single-qubit matrices. + * + * Uses the computational-basis bit order where the first operand labels the + * high bit, matching `UnitaryOpInterface::getUnitaryMatrix4x4`. + * + * @param lhs Left factor (acts on the high bit / qubit 0). + * @param rhs Right factor (acts on the low bit / qubit 1). + * @return The `4x4` Kronecker product. + */ + [[nodiscard]] static Matrix4x4 kron(const Matrix2x2& lhs, + const Matrix2x2& rhs); /** * @brief Returns the entries of column @p col, top to bottom. @@ -384,9 +427,23 @@ struct Matrix4x4 { /** * @brief Overwrites column @p col with @p values. * @param col Column index in `[0, K_COLS)`. - * @param values New column entries, top to bottom. + * @param values New column entries, top to bottom; must have length `K_ROWS`. + */ + void setColumn(std::size_t col, ArrayRef values); + + /** + * @brief Returns the entries of row @p row, left to right. + * @param row Row index in `[0, K_ROWS)`. + * @return View over the four row entries. + */ + [[nodiscard]] ArrayRef row(std::size_t row) const; + + /** + * @brief Overwrites row @p row with @p values. + * @param row Row index in `[0, K_ROWS)`. + * @param values New row entries, left to right; must have length `K_COLS`. */ - void setColumn(std::size_t col, const std::array& values); + void setRow(std::size_t row, ArrayRef values); /** * @brief Returns the element-wise real parts in row-major order. @@ -420,6 +477,58 @@ struct Matrix4x4 { * @return `true` when @p src is 4x4. */ [[nodiscard]] bool assignFrom(const DynamicMatrix& src); + + /** + * @brief Embed this two-qubit matrix into an @p numQubits-qubit Hilbert + * space. + * + * Operand 0 labels the high bit of the pair and acts on @p q0Index; operand 1 + * labels the low bit and acts on @p q1Index. For each basis pair whose other + * wires match, copies this matrix at the packed two-qubit row/column indices. + * + * @param numQubits Number of qubits in the target Hilbert space. + * @param q0Index Wire index of operand 0. + * @param q1Index Wire index of operand 1. + * @return Embedded unitary as a dynamic matrix. + */ + [[nodiscard]] DynamicMatrix embedInNqubit(std::size_t numQubits, + std::size_t q0Index, + std::size_t q1Index) const; + + /** + * @brief Reorder this matrix to act on qubits `{0, 1}`. + * + * @param q0Index Wire index of operand 0; @p q1Index wire index of operand 1. + * @return Reordered copy of this matrix. + */ + [[nodiscard]] Matrix4x4 reorderForQubits(std::size_t q0Index, + std::size_t q1Index) const; + + /** + * @brief Computes the eigendecomposition of this real symmetric matrix. + * + * @copydoc Matrix4x4::symmetricEigen4(const std::array&) + * + * @pre Entries are real (imaginary parts must be negligible). The real parts + * must form a symmetric matrix; imaginary parts are ignored. + */ + [[nodiscard]] SymmetricEigen4 symmetricEigen4() const; + + /** + * @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. + * + * @pre @p symmetric is real and symmetric: `symmetric[i,j] == symmetric[j,i]` + * for all `i, j`. Only the lower triangle (including the diagonal) is read, + * but supplying a non-symmetric matrix yields undefined numerical results. + * + * @param symmetric Row-major real symmetric `4x4` matrix. + * @return Ascending eigenvalues and matching eigenvectors (as columns). + */ + [[nodiscard]] static SymmetricEigen4 + symmetricEigen4(const std::array& symmetric); }; /** @@ -660,18 +769,6 @@ inline constexpr bool std::is_same, std::is_same>; -/** - * @brief Kronecker product `lhs (x) rhs` of two single-qubit matrices. - * - * Uses the computational-basis bit order where the first operand labels the - * high bit, matching `UnitaryOpInterface::getUnitaryMatrix4x4`. - * - * @param lhs Left factor (acts on the high bit / qubit 0). - * @param rhs Right factor (acts on the low bit / qubit 1). - * @return The `4x4` Kronecker product. - */ -[[nodiscard]] Matrix4x4 kron(const Matrix2x2& lhs, const Matrix2x2& rhs); - /// Scalar-on-the-left multiply `scalar * matrix` (commutes with the member /// `matrix * scalar`). Provided so generic code can scale a matrix from /// either side. @@ -696,62 +793,4 @@ struct SymmetricEigen4 { Matrix4x4 eigenvectors{}; }; -/** - * @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. - * - * @pre @p symmetric is real and symmetric: `symmetric[i,j] == symmetric[j,i]` - * for all `i, j`. Only the lower triangle (including the diagonal) is read, - * but supplying a non-symmetric matrix yields undefined numerical results. - * - * @param symmetric Row-major real symmetric `4x4` matrix. - * @return Ascending eigenvalues and matching eigenvectors (as columns). - */ -[[nodiscard]] SymmetricEigen4 -symmetricEigen4(const std::array& symmetric); - -/** - * @brief Computes the eigendecomposition of a real symmetric `4x4` matrix. - * - * @copydoc symmetricEigen4(const std::array&) - * - * @pre Entries of @p symmetric are real (imaginary parts must be negligible). - * The real parts must form a symmetric matrix; imaginary parts are ignored. - */ -[[nodiscard]] SymmetricEigen4 symmetricEigen4(const Matrix4x4& symmetric); - -/** - * @brief Reorder a two-qubit matrix to act on qubits `{0, 1}`. - * - * @param q0Index Wire index of operand 0; @p q1Index wire index of operand 1. - */ -[[nodiscard]] Matrix4x4 reorderTwoQubitMatrix(const Matrix4x4& matrix, - std::size_t q0Index, - std::size_t q1Index); - -/** - * @brief Embed a single-qubit matrix into an @p numQubits-qubit Hilbert space. - * - * Qubit @p qubitIndex uses the same MSB-first convention as @ref kron - * (high bit first operand, low bit second). For each basis pair whose untouched - * wires match, copies @p matrix at the target qubit's row/column bits. - */ -[[nodiscard]] DynamicMatrix embedSingleQubitInNqubit(const Matrix2x2& matrix, - std::size_t numQubits, - std::size_t qubitIndex); - -/** - * @brief Embed a two-qubit matrix into an @p numQubits-qubit Hilbert space. - * - * Operand 0 labels the high bit of the pair and acts on @p q0Index; operand 1 - * labels the low bit and acts on @p q1Index. For each basis pair whose other - * wires match, copies @p matrix at the packed two-qubit row/column indices. - */ -[[nodiscard]] DynamicMatrix embedTwoQubitInNqubit(const Matrix4x4& matrix, - std::size_t numQubits, - std::size_t q0Index, - std::size_t q1Index); - } // namespace mlir::qco diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index ed55c4a212..806aaa51ae 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -100,13 +100,6 @@ static void multiply2x2(const ArrayRef lhs, out[3] = lhs[2] * rhs[1] + lhs[3] * rhs[3]; } -static void -multiply2x2(const std::array& lhs, - const std::array& rhs, - std::array& out) { - multiply2x2(ArrayRef(lhs), ArrayRef(rhs), MutableArrayRef(out)); -} - /// Writes the row-major product `lhs * rhs` into @p out (4x4, unrolled rows). static void multiply4x4(const ArrayRef lhs, const ArrayRef rhs, @@ -127,13 +120,6 @@ static void multiply4x4(const ArrayRef lhs, } } -static void -multiply4x4(const std::array& lhs, - const std::array& rhs, - std::array& out) { - multiply4x4(ArrayRef(lhs), ArrayRef(rhs), MutableArrayRef(out)); -} - /// Returns true if @p data is approximately the @p dim x @p dim identity /// matrix. [[nodiscard]] static bool isIdentityEntries(ArrayRef data, @@ -193,10 +179,41 @@ static void copyBottomRightCorner(const std::int64_t matrixDim, } } -struct DynamicMatrix::Impl { - std::int64_t dim = 0; - SmallVector data; -}; +/** + * @brief Returns the @p qubitIndex bit of a computational-basis label. + * + * Qubit 0 is the MSB of @p stateIndex, matching @ref Matrix4x4::kron and + * @ref Matrix2x2::embedInNqubit. + */ +[[nodiscard]] static std::size_t qubitBitAt(const std::size_t stateIndex, + const std::size_t numQubits, + const std::size_t qubitIndex) { + return (stateIndex >> (numQubits - 1 - qubitIndex)) & 1U; +} + +/** + * @brief True when row and col agree on every wire except @p skipA and @p + * skipB. + * + * Used when embedding a gate: untouched qubits must match or the matrix entry + * is zero. For a single-qubit embed, pass @p skipB = @p numQubits so only @p + * skipA is skipped. + */ +[[nodiscard]] static bool otherQubitBitsMatch(const std::size_t row, + const std::size_t col, + const std::size_t numQubits, + const std::size_t skipA, + const std::size_t skipB) { + for (std::size_t q = 0; q < numQubits; ++q) { + if (q == skipA || q == skipB) { + continue; + } + if (qubitBitAt(row, numQubits, q) != qubitBitAt(col, numQubits, q)) { + return false; + } + } + return true; +} Matrix1x1 Matrix1x1::fromElements(const Complex m00) { return {m00}; } @@ -299,6 +316,40 @@ bool Matrix2x2::assignFrom(const DynamicMatrix& src) { return assignFromDynamicImpl(src, data); } +DynamicMatrix Matrix2x2::embedInNqubit(const std::size_t numQubits, + const std::size_t qubitIndex) const { + assert(qubitIndex < numQubits && + "Invalid qubit index for single-qubit embed"); + if (numQubits == 2) { + return DynamicMatrix(embedInTwoQubit(qubitIndex)); + } + const auto dim = static_cast(1ULL << numQubits); + DynamicMatrix out(dim); + const auto udim = static_cast(dim); + for (std::size_t row = 0; row < udim; ++row) { + for (std::size_t col = 0; col < udim; ++col) { + if (!otherQubitBitsMatch(row, col, numQubits, qubitIndex, numQubits)) { + continue; + } + const std::size_t rowBit = qubitBitAt(row, numQubits, qubitIndex); + const std::size_t colBit = qubitBitAt(col, numQubits, qubitIndex); + out(static_cast(row), static_cast(col)) = + (*this)(rowBit, colBit); + } + } + return out; +} + +Matrix4x4 Matrix2x2::embedInTwoQubit(const std::size_t qubitIndex) const { + if (qubitIndex == 0) { + return Matrix4x4::kron(*this, Matrix2x2::identity()); + } + if (qubitIndex == 1) { + return Matrix4x4::kron(Matrix2x2::identity(), *this); + } + assert(false && "Invalid qubit index for single-qubit embed"); +} + Matrix4x4 Matrix4x4::fromElements(const Complex& m00, const Complex& m01, const Complex& m02, const Complex& m03, const Complex& m10, const Complex& m11, @@ -385,8 +436,9 @@ std::array Matrix4x4::diagonal() const { return {data[0], data[5], data[10], data[15]}; } -Matrix4x4 -Matrix4x4::fromDiagonal(const std::array& diagonalEntries) { +Matrix4x4 Matrix4x4::fromDiagonal(const ArrayRef diagonalEntries) { + assert(diagonalEntries.size() == K_ROWS && + "fromDiagonal requires exactly K_ROWS entries"); Matrix4x4 out{}; for (std::size_t i = 0; i < K_ROWS; ++i) { out.data[(i * K_COLS) + i] = diagonalEntries[i]; @@ -394,6 +446,15 @@ Matrix4x4::fromDiagonal(const std::array& diagonalEntries) { return out; } +Matrix4x4 Matrix4x4::kron(const Matrix2x2& lhs, const Matrix2x2& rhs) { + const auto& a = lhs.data; + const auto& b = rhs.data; + return fromElements(a[0] * b[0], a[0] * b[1], a[1] * b[0], a[1] * b[1], + a[0] * b[2], a[0] * b[3], a[1] * b[2], a[1] * b[3], + a[2] * b[0], a[2] * b[1], a[3] * b[0], a[3] * b[1], + a[2] * b[2], a[2] * b[3], a[3] * b[2], a[3] * b[3]); +} + std::array Matrix4x4::column(const std::size_t col) const { return {data[col], data[K_COLS + col], data[(2 * K_COLS) + col], @@ -401,12 +462,27 @@ Matrix4x4::column(const std::size_t col) const { } void Matrix4x4::setColumn(const std::size_t col, - const std::array& values) { + const ArrayRef values) { + assert(values.size() == K_ROWS && + "setColumn requires exactly K_ROWS entries"); for (std::size_t row = 0; row < K_ROWS; ++row) { data[(row * K_COLS) + col] = values[row]; } } +ArrayRef Matrix4x4::row(const std::size_t row) const { + assert(row < K_ROWS); + return ArrayRef(data).slice(row * K_COLS, K_COLS); +} + +void Matrix4x4::setRow(const std::size_t row, const ArrayRef values) { + assert(row < K_ROWS); + assert(values.size() == K_COLS && "setRow requires exactly K_COLS entries"); + for (std::size_t col = 0; col < K_COLS; ++col) { + data[(row * K_COLS) + col] = values[col]; + } +} + std::array Matrix4x4::realPart() const { return {data[0].real(), data[1].real(), data[2].real(), data[3].real(), @@ -431,366 +507,171 @@ bool Matrix4x4::assignFrom(const DynamicMatrix& src) { return assignFromDynamicImpl(src, data); } -DynamicMatrix::DynamicMatrix() : impl_(std::make_unique()) {} - -DynamicMatrix::DynamicMatrix(const std::int64_t dim) - : impl_(std::make_unique()) { - impl_->dim = dim; - impl_->data.assign(checkedStorageSize(dim), Complex{}); +DynamicMatrix Matrix4x4::embedInNqubit(const std::size_t numQubits, + const std::size_t q0Index, + const std::size_t q1Index) const { + assert(q0Index < numQubits && q1Index < numQubits && q0Index != q1Index && + "Invalid qubit indices for two-qubit embed"); + if (numQubits == 2) { + return DynamicMatrix(reorderForQubits(q0Index, q1Index)); + } + const auto dim = static_cast(1ULL << numQubits); + DynamicMatrix out(dim); + const auto udim = static_cast(dim); + for (std::size_t row = 0; row < udim; ++row) { + for (std::size_t col = 0; col < udim; ++col) { + if (!otherQubitBitsMatch(row, col, numQubits, q0Index, q1Index)) { + continue; + } + const std::size_t rowPair = (qubitBitAt(row, numQubits, q0Index) << 1) | + qubitBitAt(row, numQubits, q1Index); + const std::size_t colPair = (qubitBitAt(col, numQubits, q0Index) << 1) | + qubitBitAt(col, numQubits, q1Index); + out(static_cast(row), static_cast(col)) = + (*this)(rowPair, colPair); + } + } + return out; } -DynamicMatrix::DynamicMatrix(const Matrix2x2& src) - : impl_(std::make_unique()) { - assignFrom(src); +Matrix4x4 Matrix4x4::reorderForQubits(const std::size_t q0Index, + const std::size_t q1Index) const { + if (q0Index == 0 && q1Index == 1) { + return *this; + } + if (q0Index == 1 && q1Index == 0) { + // Conjugate by SWAP: out[i, j] = matrix[pi(i), pi(j)] with pi swapping |01> + // and |10> (basis indices 1 and 2). + const auto& m = data; + return fromElements(m[0], m[2], m[1], m[3], m[8], m[10], m[9], m[11], m[4], + m[6], m[5], m[7], m[12], m[14], m[13], m[15]); + } + assert(false && "Invalid qubit indices for two-qubit reorder"); } -DynamicMatrix::DynamicMatrix(const Matrix4x4& src) - : impl_(std::make_unique()) { - assignFrom(src); +SymmetricEigen4 Matrix4x4::symmetricEigen4() const { + return symmetricEigen4(realPart()); } -DynamicMatrix::DynamicMatrix(const DynamicMatrix& other) - : impl_(std::make_unique(*other.impl_)) {} +// Adapted from John Burkardt's MIT-licensed EISPACK C port (`tred2` / `tql2`): +// https://people.sc.fsu.edu/~jburkardt/c_src/eispack/eispack.c +// Specialized to `n = 4`; input is row-major, accumulator `z` is column-major. -DynamicMatrix::DynamicMatrix(DynamicMatrix&& other) noexcept = default; +/// EISPACK `tred2` for `n = 4` (column-major `z[row + col*n]`). +static void symmetricTred24(const std::array& input, + std::array& z, + std::array& diag, + std::array& subdiag) { + constexpr std::size_t n = 4; + const auto zAt = [&z](const std::size_t row, + const std::size_t col) -> double& { + return z[row + (col * n)]; + }; + double h = 0.0; -DynamicMatrix& DynamicMatrix::operator=(const DynamicMatrix& other) { - if (this != &other) { - *impl_ = *other.impl_; + for (std::size_t col = 0; col < n; ++col) { + for (std::size_t row = col; row < n; ++row) { + zAt(row, col) = input[(row * n) + col]; + } + diag[col] = input[((n - 1) * n) + col]; } - return *this; -} - -DynamicMatrix& -DynamicMatrix::operator=(DynamicMatrix&& other) noexcept = default; - -DynamicMatrix::~DynamicMatrix() = default; -std::int64_t DynamicMatrix::rows() const { return impl_->dim; } - -std::int64_t DynamicMatrix::cols() const { return impl_->dim; } + for (int i = static_cast(n) - 1; i >= 1; --i) { + const auto ui = static_cast(i); + const std::size_t l = ui - 1; + h = 0.0; + double scale = 0.0; + for (std::size_t k = 0; k <= l; ++k) { + scale += std::abs(diag[k]); + } + if (scale == 0.0) { + subdiag[ui] = diag[l]; + for (std::size_t j = 0; j <= l; ++j) { + diag[j] = zAt(l, j); + zAt(ui, j) = 0.0; + zAt(j, ui) = 0.0; + } + diag[ui] = 0.0; + continue; + } + for (std::size_t k = 0; k <= l; ++k) { + diag[k] /= scale; + } + for (std::size_t k = 0; k <= l; ++k) { + h += diag[k] * diag[k]; + } + const double f = diag[l]; + const double g = -std::sqrt(h) * std::copysign(1.0, f); + subdiag[ui] = scale * g; + h -= f * g; + diag[l] = f - g; -DynamicMatrix DynamicMatrix::identity(const std::int64_t dim) { - DynamicMatrix matrix(dim); - const auto udim = checkedDim(dim); - for (std::size_t i = 0; i < udim; ++i) { - matrix.impl_->data[(i * udim) + i] = 1.0; + for (std::size_t k = 0; k <= l; ++k) { + subdiag[k] = 0.0; + } + for (std::size_t j = 0; j <= l; ++j) { + const double fj = diag[j]; + zAt(j, ui) = fj; + double gj = subdiag[j] + (zAt(j, j) * fj); + for (std::size_t k = j + 1; k <= l; ++k) { + gj += zAt(k, j) * diag[k]; + subdiag[k] += zAt(k, j) * fj; + } + subdiag[j] = gj; + } + double ff = 0.0; + for (std::size_t k = 0; k <= l; ++k) { + subdiag[k] /= h; + ff += subdiag[k] * diag[k]; + } + const double hh = 0.5 * ff / h; + for (std::size_t k = 0; k <= l; ++k) { + subdiag[k] -= hh * diag[k]; + } + for (std::size_t j = 0; j <= l; ++j) { + const double fj = diag[j]; + const double gj = subdiag[j]; + for (std::size_t k = j; k <= l; ++k) { + zAt(k, j) -= (fj * subdiag[k]) + (gj * diag[k]); + } + diag[j] = zAt(l, j); + zAt(ui, j) = 0.0; + } + diag[ui] = h; } - return matrix; -} -DynamicMatrix DynamicMatrix::fromAdjoint(const Matrix2x2& src) { - return DynamicMatrix(src.adjoint()); -} + for (std::size_t i = 1; i < n; ++i) { + const std::size_t l = i - 1; + zAt(n - 1, l) = zAt(l, l); + zAt(l, l) = 1.0; + h = diag[i]; + if (h != 0.0) { + for (std::size_t k = 0; k <= l; ++k) { + diag[k] = zAt(k, i) / h; + } + for (std::size_t j = 0; j <= l; ++j) { + double g = 0.0; + for (std::size_t k = 0; k <= l; ++k) { + g += zAt(k, i) * zAt(k, j); + } + for (std::size_t k = 0; k <= l; ++k) { + zAt(k, j) -= g * diag[k]; + } + } + } + for (std::size_t k = 0; k <= l; ++k) { + zAt(k, i) = 0.0; + } + } -Complex& DynamicMatrix::operator()(const std::int64_t row, - const std::int64_t col) { - return impl_->data[static_cast((row * impl_->dim) + col)]; -} - -Complex DynamicMatrix::operator()(const std::int64_t row, - const std::int64_t col) const { - return impl_->data[static_cast((row * impl_->dim) + col)]; -} - -void DynamicMatrix::setBottomRightCorner(const Matrix2x2& block) { - copyBottomRightCorner(impl_->dim, impl_->data, - static_cast(Matrix2x2::K_ROWS), - block.data); -} - -void DynamicMatrix::setBottomRightCorner(const Matrix4x4& block) { - copyBottomRightCorner(impl_->dim, impl_->data, - static_cast(Matrix4x4::K_ROWS), - block.data); -} - -void DynamicMatrix::setBottomRightCorner(const DynamicMatrix& block) { - copyBottomRightCorner(impl_->dim, impl_->data, block.impl_->dim, - block.impl_->data); -} - -DynamicMatrix DynamicMatrix::adjoint() const { - DynamicMatrix out(impl_->dim); - adjointInto(impl_->data, out.impl_->data, checkedDim(impl_->dim)); - return out; -} - -void DynamicMatrix::assignFrom(const Matrix1x1& src) { - impl_->dim = 1; - impl_->data.assign({src.value}); -} - -void DynamicMatrix::assignFrom(const Matrix2x2& src) { - assignFixedImpl( - impl_->dim, impl_->data, src.data); -} - -void DynamicMatrix::assignFrom(const Matrix4x4& src) { - assignFixedImpl( - impl_->dim, impl_->data, src.data); -} - -void DynamicMatrix::assignFrom(const DynamicMatrix& src) { - *impl_ = *src.impl_; -} - -bool DynamicMatrix::isApprox(const Matrix1x1& other, const double tol) const { - if (impl_->dim != 1) { - return false; - } - return entryIsApprox(impl_->data[0], other.value, tol); -} - -bool DynamicMatrix::isApprox(const Matrix2x2& other, const double tol) const { - return isApproxFixedImpl( - impl_->dim, impl_->data, other.data, tol); -} - -bool DynamicMatrix::isApprox(const Matrix4x4& other, const double tol) const { - return isApproxFixedImpl( - impl_->dim, impl_->data, other.data, tol); -} - -bool DynamicMatrix::isApprox(const DynamicMatrix& other, - const double tol) const { - return entriesAreApprox(impl_->data, other.impl_->data, tol); -} - -Complex DynamicMatrix::trace() const { - Complex sum{0.0, 0.0}; - const auto udim = checkedDim(impl_->dim); - for (std::size_t i = 0; i < udim; ++i) { - sum += impl_->data[(i * udim) + i]; - } - return sum; -} - -DynamicMatrix DynamicMatrix::operator*(const DynamicMatrix& rhs) const { - assert(impl_->dim == rhs.impl_->dim && - "DynamicMatrix multiply requires matching dimensions"); - DynamicMatrix out(impl_->dim); - if (std::cmp_equal(impl_->dim, Matrix2x2::K_ROWS)) { - multiply2x2(impl_->data, rhs.impl_->data, out.impl_->data); - return out; - } - if (std::cmp_equal(impl_->dim, Matrix4x4::K_ROWS)) { - multiply4x4(impl_->data, rhs.impl_->data, out.impl_->data); - return out; - } - - const auto udim = checkedDim(impl_->dim); - for (std::size_t row = 0; row < udim; ++row) { - for (std::size_t col = 0; col < udim; ++col) { - Complex sum{0.0, 0.0}; - for (std::size_t k = 0; k < udim; ++k) { - sum += - impl_->data[(row * udim) + k] * rhs.impl_->data[(k * udim) + col]; - } - out.impl_->data[(row * udim) + col] = sum; - } - } - return out; -} - -DynamicMatrix DynamicMatrix::operator*(const Complex& scalar) const { - DynamicMatrix out(impl_->dim); - for (std::size_t i = 0; i < impl_->data.size(); ++i) { - out.impl_->data[i] = impl_->data[i] * scalar; - } - return out; -} - -DynamicMatrix& DynamicMatrix::operator*=(const Complex& scalar) { - for (Complex& entry : impl_->data) { - entry *= scalar; - } - return *this; -} - -bool DynamicMatrix::isIdentity(const double tol) const { - return isIdentityEntries(impl_->data, checkedDim(impl_->dim), tol); -} - -/** - * @brief Returns the @p qubitIndex bit of a computational-basis label. - * - * Qubit 0 is the MSB of @p stateIndex, matching @ref kron and - * @ref embedSingleQubitInNqubit. - */ -[[nodiscard]] static std::size_t qubitBitAt(const std::size_t stateIndex, - const std::size_t numQubits, - const std::size_t qubitIndex) { - return (stateIndex >> (numQubits - 1 - qubitIndex)) & 1U; -} - -/** - * @brief True when row and col agree on every wire except @p skipA and @p - * skipB. - * - * Used when embedding a gate: untouched qubits must match or the matrix entry - * is zero. For a single-qubit embed, pass @p skipB = @p numQubits so only @p - * skipA is skipped. - */ -[[nodiscard]] static bool otherQubitBitsMatch(const std::size_t row, - const std::size_t col, - const std::size_t numQubits, - const std::size_t skipA, - const std::size_t skipB) { - for (std::size_t q = 0; q < numQubits; ++q) { - if (q == skipA || q == skipB) { - continue; - } - if (qubitBitAt(row, numQubits, q) != qubitBitAt(col, numQubits, q)) { - return false; - } - } - return true; -} - -Matrix4x4 kron(const Matrix2x2& lhs, const Matrix2x2& rhs) { - const auto& a = lhs.data; - const auto& b = rhs.data; - return Matrix4x4::fromElements( - a[0] * b[0], a[0] * b[1], a[1] * b[0], a[1] * b[1], a[0] * b[2], - a[0] * b[3], a[1] * b[2], a[1] * b[3], a[2] * b[0], a[2] * b[1], - a[3] * b[0], a[3] * b[1], a[2] * b[2], a[2] * b[3], a[3] * b[2], - a[3] * b[3]); -} - -Matrix2x2 operator*(const Complex& scalar, const Matrix2x2& matrix) { - return matrix * scalar; -} - -Matrix4x4 operator*(const Complex& scalar, const Matrix4x4& matrix) { - return matrix * scalar; -} - -DynamicMatrix operator*(const Complex& scalar, const DynamicMatrix& matrix) { - return matrix * scalar; -} - -// Adapted from John Burkardt's MIT-licensed EISPACK C port (`tred2` / `tql2`): -// https://people.sc.fsu.edu/~jburkardt/c_src/eispack/eispack.c -// Specialized to `n = 4`; input is row-major, accumulator `z` is column-major. - -/// EISPACK `tred2` for `n = 4` (column-major `z[row + col*n]`). -static void symmetricTred24(const std::array& input, - std::array& z, - std::array& diag, - std::array& subdiag) { - constexpr std::size_t n = 4; - const auto zAt = [&z](const std::size_t row, - const std::size_t col) -> double& { - return z[row + (col * n)]; - }; - double h = 0.0; - - for (std::size_t col = 0; col < n; ++col) { - for (std::size_t row = col; row < n; ++row) { - zAt(row, col) = input[(row * n) + col]; - } - diag[col] = input[((n - 1) * n) + col]; - } - - for (int i = static_cast(n) - 1; i >= 1; --i) { - const auto ui = static_cast(i); - const std::size_t l = ui - 1; - h = 0.0; - double scale = 0.0; - for (std::size_t k = 0; k <= l; ++k) { - scale += std::abs(diag[k]); - } - if (scale == 0.0) { - subdiag[ui] = diag[l]; - for (std::size_t j = 0; j <= l; ++j) { - diag[j] = zAt(l, j); - zAt(ui, j) = 0.0; - zAt(j, ui) = 0.0; - } - diag[ui] = 0.0; - continue; - } - for (std::size_t k = 0; k <= l; ++k) { - diag[k] /= scale; - } - for (std::size_t k = 0; k <= l; ++k) { - h += diag[k] * diag[k]; - } - const double f = diag[l]; - const double g = -std::sqrt(h) * std::copysign(1.0, f); - subdiag[ui] = scale * g; - h -= f * g; - diag[l] = f - g; - - for (std::size_t k = 0; k <= l; ++k) { - subdiag[k] = 0.0; - } - for (std::size_t j = 0; j <= l; ++j) { - const double fj = diag[j]; - zAt(j, ui) = fj; - double gj = subdiag[j] + (zAt(j, j) * fj); - for (std::size_t k = j + 1; k <= l; ++k) { - gj += zAt(k, j) * diag[k]; - subdiag[k] += zAt(k, j) * fj; - } - subdiag[j] = gj; - } - double ff = 0.0; - for (std::size_t k = 0; k <= l; ++k) { - subdiag[k] /= h; - ff += subdiag[k] * diag[k]; - } - const double hh = 0.5 * ff / h; - for (std::size_t k = 0; k <= l; ++k) { - subdiag[k] -= hh * diag[k]; - } - for (std::size_t j = 0; j <= l; ++j) { - const double fj = diag[j]; - const double gj = subdiag[j]; - for (std::size_t k = j; k <= l; ++k) { - zAt(k, j) -= (fj * subdiag[k]) + (gj * diag[k]); - } - diag[j] = zAt(l, j); - zAt(ui, j) = 0.0; - } - diag[ui] = h; - } - - for (std::size_t i = 1; i < n; ++i) { - const std::size_t l = i - 1; - zAt(n - 1, l) = zAt(l, l); - zAt(l, l) = 1.0; - h = diag[i]; - if (h != 0.0) { - for (std::size_t k = 0; k <= l; ++k) { - diag[k] = zAt(k, i) / h; - } - for (std::size_t j = 0; j <= l; ++j) { - double g = 0.0; - for (std::size_t k = 0; k <= l; ++k) { - g += zAt(k, i) * zAt(k, j); - } - for (std::size_t k = 0; k <= l; ++k) { - zAt(k, j) -= g * diag[k]; - } - } - } - for (std::size_t k = 0; k <= l; ++k) { - zAt(k, i) = 0.0; - } - } - - for (std::size_t j = 0; j < n; ++j) { - diag[j] = zAt(n - 1, j); - } - for (std::size_t j = 0; j < n - 1; ++j) { - zAt(n - 1, j) = 0.0; - } - zAt(n - 1, n - 1) = 1.0; - subdiag[0] = 0.0; + for (std::size_t j = 0; j < n; ++j) { + diag[j] = zAt(n - 1, j); + } + for (std::size_t j = 0; j < n - 1; ++j) { + zAt(n - 1, j) = 0.0; + } + zAt(n - 1, n - 1) = 1.0; + subdiag[0] = 0.0; } /// EISPACK `tql2` for `n = 4` (column-major `z[row + col*n]`). @@ -907,18 +788,17 @@ static void symmetricTql24(std::array& diag, } } -SymmetricEigen4 symmetricEigen4(const std::array& symmetric) { +SymmetricEigen4 +Matrix4x4::symmetricEigen4(const std::array& symmetric) { constexpr std::size_t n = 4; + SymmetricEigen4 result; std::array z{}; - std::array diag{}; std::array subdiag{}; - symmetricTred24(symmetric, z, diag, subdiag); - symmetricTql24(diag, subdiag, z); + symmetricTred24(symmetric, z, result.eigenvalues, subdiag); + symmetricTql24(result.eigenvalues, subdiag, z); - SymmetricEigen4 result; for (std::size_t col = 0; col < n; ++col) { - result.eigenvalues[col] = diag[col]; for (std::size_t row = 0; row < n; ++row) { result.eigenvectors(row, col) = Complex{z[row + (col * n)], 0.0}; } @@ -926,89 +806,204 @@ SymmetricEigen4 symmetricEigen4(const std::array& symmetric) { return result; } -SymmetricEigen4 symmetricEigen4(const Matrix4x4& symmetric) { - return symmetricEigen4(symmetric.realPart()); +struct DynamicMatrix::Impl { + std::int64_t dim = 0; + SmallVector data; +}; + +DynamicMatrix::DynamicMatrix() : impl_(std::make_unique()) {} + +DynamicMatrix::DynamicMatrix(const std::int64_t dim) + : impl_(std::make_unique()) { + impl_->dim = dim; + impl_->data.assign(checkedStorageSize(dim), Complex{}); +} + +DynamicMatrix::DynamicMatrix(const Matrix2x2& src) + : impl_(std::make_unique()) { + assignFrom(src); +} + +DynamicMatrix::DynamicMatrix(const Matrix4x4& src) + : impl_(std::make_unique()) { + assignFrom(src); } -static Matrix4x4 embedSingleQubitInTwoQubit(const Matrix2x2& matrix, - const std::size_t qubitIndex) { - if (qubitIndex == 0) { - return kron(matrix, Matrix2x2::identity()); +DynamicMatrix::DynamicMatrix(const DynamicMatrix& other) + : impl_(std::make_unique(*other.impl_)) {} + +DynamicMatrix::DynamicMatrix(DynamicMatrix&& other) noexcept = default; + +DynamicMatrix& DynamicMatrix::operator=(const DynamicMatrix& other) { + if (this != &other) { + *impl_ = *other.impl_; } - if (qubitIndex == 1) { - return kron(Matrix2x2::identity(), matrix); + return *this; +} + +DynamicMatrix& +DynamicMatrix::operator=(DynamicMatrix&& other) noexcept = default; + +DynamicMatrix::~DynamicMatrix() = default; + +std::int64_t DynamicMatrix::rows() const { return impl_->dim; } + +std::int64_t DynamicMatrix::cols() const { return impl_->dim; } + +DynamicMatrix DynamicMatrix::identity(const std::int64_t dim) { + DynamicMatrix matrix(dim); + const auto udim = checkedDim(dim); + for (std::size_t i = 0; i < udim; ++i) { + matrix.impl_->data[(i * udim) + i] = 1.0; } - assert(false && "Invalid qubit index for single-qubit embed"); + return matrix; } -Matrix4x4 reorderTwoQubitMatrix(const Matrix4x4& matrix, - const std::size_t q0Index, - const std::size_t q1Index) { - if (q0Index == 0 && q1Index == 1) { - return matrix; +DynamicMatrix DynamicMatrix::fromAdjoint(const Matrix2x2& src) { + return DynamicMatrix(src.adjoint()); +} + +Complex& DynamicMatrix::operator()(const std::int64_t row, + const std::int64_t col) { + return impl_->data[static_cast((row * impl_->dim) + col)]; +} + +Complex DynamicMatrix::operator()(const std::int64_t row, + const std::int64_t col) const { + return impl_->data[static_cast((row * impl_->dim) + col)]; +} + +void DynamicMatrix::setBottomRightCorner(const Matrix2x2& block) { + copyBottomRightCorner(impl_->dim, impl_->data, + static_cast(Matrix2x2::K_ROWS), + block.data); +} + +void DynamicMatrix::setBottomRightCorner(const Matrix4x4& block) { + copyBottomRightCorner(impl_->dim, impl_->data, + static_cast(Matrix4x4::K_ROWS), + block.data); +} + +void DynamicMatrix::setBottomRightCorner(const DynamicMatrix& block) { + copyBottomRightCorner(impl_->dim, impl_->data, block.impl_->dim, + block.impl_->data); +} + +DynamicMatrix DynamicMatrix::adjoint() const { + DynamicMatrix out(impl_->dim); + adjointInto(impl_->data, out.impl_->data, checkedDim(impl_->dim)); + return out; +} + +void DynamicMatrix::assignFrom(const Matrix1x1& src) { + impl_->dim = 1; + impl_->data.assign({src.value}); +} + +void DynamicMatrix::assignFrom(const Matrix2x2& src) { + assignFixedImpl( + impl_->dim, impl_->data, src.data); +} + +void DynamicMatrix::assignFrom(const Matrix4x4& src) { + assignFixedImpl( + impl_->dim, impl_->data, src.data); +} + +void DynamicMatrix::assignFrom(const DynamicMatrix& src) { + *impl_ = *src.impl_; +} + +bool DynamicMatrix::isApprox(const Matrix1x1& other, const double tol) const { + if (impl_->dim != 1) { + return false; } - if (q0Index == 1 && q1Index == 0) { - // Conjugate by SWAP: out[i, j] = matrix[pi(i), pi(j)] with pi swapping |01> - // and |10> (basis indices 1 and 2). - const auto& m = matrix.data; - return Matrix4x4::fromElements(m[0], m[2], m[1], m[3], m[8], m[10], m[9], - m[11], m[4], m[6], m[5], m[7], m[12], m[14], - m[13], m[15]); + return entryIsApprox(impl_->data[0], other.value, tol); +} + +bool DynamicMatrix::isApprox(const Matrix2x2& other, const double tol) const { + return isApproxFixedImpl( + impl_->dim, impl_->data, other.data, tol); +} + +bool DynamicMatrix::isApprox(const Matrix4x4& other, const double tol) const { + return isApproxFixedImpl( + impl_->dim, impl_->data, other.data, tol); +} + +bool DynamicMatrix::isApprox(const DynamicMatrix& other, + const double tol) const { + return entriesAreApprox(impl_->data, other.impl_->data, tol); +} + +Complex DynamicMatrix::trace() const { + Complex sum{0.0, 0.0}; + const auto udim = checkedDim(impl_->dim); + for (std::size_t i = 0; i < udim; ++i) { + sum += impl_->data[(i * udim) + i]; } - assert(false && "Invalid qubit indices for two-qubit reorder"); + return sum; } -DynamicMatrix embedSingleQubitInNqubit(const Matrix2x2& matrix, - const std::size_t numQubits, - const std::size_t qubitIndex) { - assert(qubitIndex < numQubits && - "Invalid qubit index for single-qubit embed"); - if (numQubits == 2) { - return DynamicMatrix(embedSingleQubitInTwoQubit(matrix, qubitIndex)); +DynamicMatrix DynamicMatrix::operator*(const DynamicMatrix& rhs) const { + assert(impl_->dim == rhs.impl_->dim && + "DynamicMatrix multiply requires matching dimensions"); + DynamicMatrix out(impl_->dim); + if (std::cmp_equal(impl_->dim, Matrix2x2::K_ROWS)) { + multiply2x2(impl_->data, rhs.impl_->data, out.impl_->data); + return out; } - const auto dim = static_cast(1ULL << numQubits); - DynamicMatrix out(dim); - const auto udim = static_cast(dim); + if (std::cmp_equal(impl_->dim, Matrix4x4::K_ROWS)) { + multiply4x4(impl_->data, rhs.impl_->data, out.impl_->data); + return out; + } + + const auto udim = checkedDim(impl_->dim); for (std::size_t row = 0; row < udim; ++row) { for (std::size_t col = 0; col < udim; ++col) { - if (!otherQubitBitsMatch(row, col, numQubits, qubitIndex, numQubits)) { - continue; + Complex sum{0.0, 0.0}; + for (std::size_t k = 0; k < udim; ++k) { + sum += + impl_->data[(row * udim) + k] * rhs.impl_->data[(k * udim) + col]; } - const std::size_t rowBit = qubitBitAt(row, numQubits, qubitIndex); - const std::size_t colBit = qubitBitAt(col, numQubits, qubitIndex); - out(static_cast(row), static_cast(col)) = - matrix(rowBit, colBit); + out.impl_->data[(row * udim) + col] = sum; } } return out; } -DynamicMatrix embedTwoQubitInNqubit(const Matrix4x4& matrix, - const std::size_t numQubits, - const std::size_t q0Index, - const std::size_t q1Index) { - assert(q0Index < numQubits && q1Index < numQubits && q0Index != q1Index && - "Invalid qubit indices for two-qubit embed"); - if (numQubits == 2) { - return DynamicMatrix(reorderTwoQubitMatrix(matrix, q0Index, q1Index)); - } - const auto dim = static_cast(1ULL << numQubits); - DynamicMatrix out(dim); - const auto udim = static_cast(dim); - for (std::size_t row = 0; row < udim; ++row) { - for (std::size_t col = 0; col < udim; ++col) { - if (!otherQubitBitsMatch(row, col, numQubits, q0Index, q1Index)) { - continue; - } - const std::size_t rowPair = (qubitBitAt(row, numQubits, q0Index) << 1) | - qubitBitAt(row, numQubits, q1Index); - const std::size_t colPair = (qubitBitAt(col, numQubits, q0Index) << 1) | - qubitBitAt(col, numQubits, q1Index); - out(static_cast(row), static_cast(col)) = - matrix(rowPair, colPair); - } +DynamicMatrix DynamicMatrix::operator*(const Complex& scalar) const { + DynamicMatrix out(impl_->dim); + for (std::size_t i = 0; i < impl_->data.size(); ++i) { + out.impl_->data[i] = impl_->data[i] * scalar; } return out; } +DynamicMatrix& DynamicMatrix::operator*=(const Complex& scalar) { + for (Complex& entry : impl_->data) { + entry *= scalar; + } + return *this; +} + +bool DynamicMatrix::isIdentity(const double tol) const { + return isIdentityEntries(impl_->data, checkedDim(impl_->dim), tol); +} + +Matrix2x2 operator*(const Complex& scalar, const Matrix2x2& matrix) { + return matrix * scalar; +} + +Matrix4x4 operator*(const Complex& scalar, const Matrix4x4& matrix) { + return matrix * scalar; +} + +DynamicMatrix operator*(const Complex& scalar, const DynamicMatrix& matrix) { + return matrix * scalar; +} + } // namespace mlir::qco diff --git a/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp b/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp index 946dc1956e..c23426918b 100644 --- a/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp +++ b/mlir/unittests/Dialect/QCO/Utils/test_unitary_matrix.cpp @@ -370,7 +370,7 @@ TEST(UnitaryMatrix4x4, TransposeAndIsIdentity) { EXPECT_FALSE(swapMatrix().isIdentity()); } -TEST(UnitaryMatrix4x4, DiagonalColumnsAndParts) { +TEST(UnitaryMatrix4x4, DiagonalRowsColumnsAndParts) { Matrix4x4 m = Matrix4x4::fromElements(Complex{1, 1}, 0, 0, 0, 0, Complex{2, 2}, 0, 0, 0, 0, Complex{3, 3}, 0, 0, 0, 0, Complex{4, 4}); @@ -386,6 +386,19 @@ TEST(UnitaryMatrix4x4, DiagonalColumnsAndParts) { EXPECT_EQ(n(0, 2), 1i); EXPECT_EQ(n(3, 2), 4i); + const auto row1 = m.row(1); + ASSERT_EQ(row1.size(), Matrix4x4::K_COLS); + EXPECT_EQ(row1[0], (Complex{0, 0})); + EXPECT_EQ(row1[1], (Complex{2, 2})); + EXPECT_EQ(row1[0], m(1, 0)); + EXPECT_EQ(row1[3], m(1, 3)); + + Matrix4x4 r = Matrix4x4::identity(); + r.setRow(2, {1.0, 2.0, 3.0, 4.0}); + EXPECT_EQ(r(2, 0), 1.0); + EXPECT_EQ(r(2, 3), 4.0); + EXPECT_EQ(r.row(2)[1], 2.0); + const auto re = m.realPart(); const auto im = m.imagPart(); EXPECT_EQ(re[0], 1.0); @@ -397,13 +410,13 @@ TEST(UnitaryMatrix4x4, DiagonalColumnsAndParts) { TEST(UnitaryMatrix4x4, KroneckerProduct) { const Matrix2x2 x = pauliX(); // X (x) I should swap the high bit. - const Matrix4x4 xi = kron(x, Matrix2x2::identity()); + const Matrix4x4 xi = Matrix4x4::kron(x, Matrix2x2::identity()); EXPECT_TRUE(xi.isApprox(Matrix4x4::fromElements(0, 0, 1, 0, // row 0 0, 0, 0, 1, // row 1 1, 0, 0, 0, // row 2 0, 1, 0, 0))); // I (x) X swaps the low bit. - const Matrix4x4 ix = kron(Matrix2x2::identity(), x); + const Matrix4x4 ix = Matrix4x4::kron(Matrix2x2::identity(), x); EXPECT_TRUE(ix.isApprox(Matrix4x4::fromElements(0, 1, 0, 0, // row 0 1, 0, 0, 0, // row 1 0, 0, 0, 1, // row 2 @@ -412,12 +425,12 @@ TEST(UnitaryMatrix4x4, KroneckerProduct) { TEST(UnitaryMatrix4x4, ReorderTwoQubitMatrix) { const Matrix2x2 x = pauliX(); - const Matrix4x4 onHigh = kron(x, Matrix2x2::identity()); - const Matrix4x4 onLow = kron(Matrix2x2::identity(), x); + const Matrix4x4 onHigh = Matrix4x4::kron(x, Matrix2x2::identity()); + const Matrix4x4 onLow = Matrix4x4::kron(Matrix2x2::identity(), x); - EXPECT_TRUE(reorderTwoQubitMatrix(onHigh, 0, 1).isApprox(onHigh)); - EXPECT_TRUE(reorderTwoQubitMatrix(onHigh, 1, 0).isApprox(onLow)); - EXPECT_TRUE(reorderTwoQubitMatrix(onLow, 1, 0).isApprox(onHigh)); + EXPECT_TRUE(onHigh.reorderForQubits(0, 1).isApprox(onHigh)); + EXPECT_TRUE(onHigh.reorderForQubits(1, 0).isApprox(onLow)); + EXPECT_TRUE(onLow.reorderForQubits(1, 0).isApprox(onHigh)); } TEST(UnitaryDynamicMatrix, NQubitEmbedMatchesTwoQubitSpecialization) { @@ -426,14 +439,13 @@ TEST(UnitaryDynamicMatrix, NQubitEmbedMatchesTwoQubitSpecialization) { 0, 1, 0, 0, // 0, 0, 0, 1, // 0, 0, 1, 0); - EXPECT_TRUE(embedSingleQubitInNqubit(x, 2, 0).isApprox( - kron(x, Matrix2x2::identity()))); - EXPECT_TRUE(embedSingleQubitInNqubit(x, 2, 1).isApprox( - kron(Matrix2x2::identity(), x))); - EXPECT_TRUE(embedTwoQubitInNqubit(cx, 2, 0, 1) - .isApprox(reorderTwoQubitMatrix(cx, 0, 1))); - const DynamicMatrix cxOn01 = embedTwoQubitInNqubit(cx, 3, 0, 1); - const DynamicMatrix cxOn12 = embedTwoQubitInNqubit(cx, 3, 1, 2); + EXPECT_TRUE(x.embedInNqubit(2, 0).isApprox( + Matrix4x4::kron(x, Matrix2x2::identity()))); + EXPECT_TRUE(x.embedInNqubit(2, 1).isApprox( + Matrix4x4::kron(Matrix2x2::identity(), x))); + EXPECT_TRUE(cx.embedInNqubit(2, 0, 1).isApprox(cx.reorderForQubits(0, 1))); + const DynamicMatrix cxOn01 = cx.embedInNqubit(3, 0, 1); + const DynamicMatrix cxOn12 = cx.embedInNqubit(3, 1, 2); EXPECT_EQ(cxOn01.rows(), 8); EXPECT_EQ(cxOn12.rows(), 8); EXPECT_FALSE(cxOn01.isApprox(cxOn12)); @@ -441,7 +453,7 @@ TEST(UnitaryDynamicMatrix, NQubitEmbedMatchesTwoQubitSpecialization) { TEST(UnitaryDynamicMatrix, EmbedSingleQubitOnMiddleWire) { const Matrix2x2 x = pauliX(); - const DynamicMatrix embedded = embedSingleQubitInNqubit(x, 3, 1); + const DynamicMatrix embedded = x.embedInNqubit(3, 1); EXPECT_EQ(embedded.rows(), 8); EXPECT_FALSE(embedded.isIdentity()); @@ -452,7 +464,7 @@ TEST(UnitaryDynamicMatrix, EmbedSingleQubitOnMiddleWire) { TEST(UnitaryDynamicMatrix, MultiplyTraceAndScalar) { const Matrix2x2 x = pauliX(); - const DynamicMatrix embedded = embedSingleQubitInNqubit(x, 2, 0); + const DynamicMatrix embedded = x.embedInNqubit(2, 0); EXPECT_FALSE(embedded.isIdentity()); const Complex scalar = std::exp(1i * 0.3); EXPECT_TRUE((scalar * embedded).isApprox(embedded * scalar)); @@ -506,7 +518,7 @@ TEST(SymmetricEigensolver, DiagonalMatrix) { a[5] = 1.0; a[10] = 4.0; a[15] = 2.0; - const SymmetricEigen4 result = symmetricEigen4(a); + const SymmetricEigen4 result = Matrix4x4::symmetricEigen4(a); EXPECT_NEAR(result.eigenvalues[0], 1.0, MATRIX_TOLERANCE); EXPECT_NEAR(result.eigenvalues[1], 2.0, MATRIX_TOLERANCE); EXPECT_NEAR(result.eigenvalues[2], 3.0, MATRIX_TOLERANCE); @@ -523,8 +535,8 @@ TEST(SymmetricEigensolver, Matrix4x4Overload) { for (std::size_t k = 0; k < 16; ++k) { matrix(k / 4, k % 4) = a[k]; } - const SymmetricEigen4 fromArray = symmetricEigen4(a); - const SymmetricEigen4 fromMatrix = symmetricEigen4(matrix); + const SymmetricEigen4 fromArray = Matrix4x4::symmetricEigen4(a); + const SymmetricEigen4 fromMatrix = matrix.symmetricEigen4(); for (std::size_t i = 0; i < 4; ++i) { EXPECT_NEAR(fromMatrix.eigenvalues[i], fromArray.eigenvalues[i], MATRIX_TOLERANCE); @@ -544,7 +556,7 @@ TEST(SymmetricEigensolver, ReconstructsRandomSymmetric) { a[(j * 4) + i] = value; } } - const SymmetricEigen4 result = symmetricEigen4(a); + const SymmetricEigen4 result = Matrix4x4::symmetricEigen4(a); // Eigenvalues are ascending. for (std::size_t i = 0; i + 1 < 4; ++i) { @@ -576,7 +588,7 @@ TEST(SymmetricEigensolver, HandlesDegenerateSpectrum) { for (std::size_t i = 0; i < 4; ++i) { a[(i * 4) + i] = 2.5; } - const SymmetricEigen4 result = symmetricEigen4(a); + const SymmetricEigen4 result = Matrix4x4::symmetricEigen4(a); for (const double value : result.eigenvalues) { EXPECT_NEAR(value, 2.5, MATRIX_TOLERANCE); } From 01cfccbcaa28ed7c232bc0ae88af7731aa7e0a64 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 13:43:31 +0200 Subject: [PATCH 05/10] =?UTF-8?q?=F0=9F=93=9D=20Update=20CHANGELOG=20to=20?= =?UTF-8?q?include=20PR=20#1802?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eb58acb8be..614fe3c174 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -42,7 +42,7 @@ with the exception that minor releases may include breaking changes. [#1569], [#1570], [#1572], [#1573], [#1580], [#1602], [#1620], [#1623], [#1624], [#1626], [#1627], [#1635], [#1638], [#1673], [#1675], [#1700], [#1710], [#1717], [#1728], [#1730], [#1749], [#1751], [#1762], [#1765], - [#1774], [#1780], [#1781], [#1782], [#1787]) + [#1774], [#1780], [#1781], [#1782], [#1787], [#1802]) ([**@burgholzer**], [**@denialhaag**], [**@taminob**], [**@DRovara**], [**@li-mingbao**], [**@Ectras**], [**@MatthiasReumann**], [**@simon1hofmann**]) @@ -598,6 +598,7 @@ changelogs._ +[#1802]: https://github.com/munich-quantum-toolkit/core/pull/1802 [#1787]: https://github.com/munich-quantum-toolkit/core/pull/1787 [#1782]: https://github.com/munich-quantum-toolkit/core/pull/1782 [#1781]: https://github.com/munich-quantum-toolkit/core/pull/1781 From aa593accc2f8245e0422c109a678767a4e6587bc Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 13:54:18 +0200 Subject: [PATCH 06/10] =?UTF-8?q?=F0=9F=9A=A8=20Fix=20linter=20warnings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 806aaa51ae..aeea248840 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -11,6 +11,7 @@ #include "mlir/Dialect/QCO/Utils/Matrix.h" #include +#include #include #include @@ -271,10 +272,7 @@ Matrix2x2 Matrix2x2::operator*(const Matrix2x2& rhs) const { return out; } -void Matrix2x2::premultiplyBy(const Matrix2x2& lhs) { - const std::array rhs = data; - multiply2x2(lhs.data, rhs, data); -} +void Matrix2x2::premultiplyBy(const Matrix2x2& lhs) { *this = lhs * *this; } Matrix2x2 Matrix2x2::operator*(const Complex& scalar) const { Matrix2x2 out = *this; @@ -377,10 +375,7 @@ Matrix4x4 Matrix4x4::operator*(const Matrix4x4& rhs) const { return out; } -void Matrix4x4::premultiplyBy(const Matrix4x4& lhs) { - const std::array rhs = data; - multiply4x4(lhs.data, rhs, data); -} +void Matrix4x4::premultiplyBy(const Matrix4x4& lhs) { *this = lhs * *this; } Matrix4x4 Matrix4x4::operator*(const Complex& scalar) const { Matrix4x4 out = *this; From 3f38d9c34e1754bfa155fa0232a39b4a066d68f3 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 14:11:05 +0200 Subject: [PATCH 07/10] =?UTF-8?q?=F0=9F=90=87=20Address=20rabbit's=20comme?= =?UTF-8?q?nts?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/include/mlir/Dialect/QCO/Utils/Matrix.h | 2 ++ mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 20 +++++++++++++++----- 2 files changed, 17 insertions(+), 5 deletions(-) diff --git a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h index 369cd2209b..cdea703ed4 100644 --- a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h +++ b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h @@ -789,7 +789,9 @@ inline constexpr bool * eigenvector for `eigenvalues[j]`). */ struct SymmetricEigen4 { + /// Eigenvalues in ascending order. std::array eigenvalues{}; + /// Orthonormal eigenvectors as columns (column `j` matches `eigenvalues[j]`). Matrix4x4 eigenvectors{}; }; diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index aeea248840..00b590271b 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -155,6 +155,14 @@ static void multiply4x4(const ArrayRef lhs, return udim * udim; } +/// Returns `2^numQubits` as `int64_t` after checking it fits. +[[nodiscard]] static std::int64_t +checkedHilbertDim(const std::size_t numQubits) { + assert(numQubits < std::numeric_limits::digits && + "Hilbert-space dimension must fit in int64_t"); + return std::int64_t{static_cast(std::uint64_t{1} << numQubits)}; +} + static void validateCornerDims(const std::int64_t matrixDim, const std::int64_t blockDim) { assert(matrixDim >= 0 && blockDim >= 0 && blockDim <= matrixDim && @@ -321,7 +329,7 @@ DynamicMatrix Matrix2x2::embedInNqubit(const std::size_t numQubits, if (numQubits == 2) { return DynamicMatrix(embedInTwoQubit(qubitIndex)); } - const auto dim = static_cast(1ULL << numQubits); + const auto dim = checkedHilbertDim(numQubits); DynamicMatrix out(dim); const auto udim = static_cast(dim); for (std::size_t row = 0; row < udim; ++row) { @@ -345,7 +353,7 @@ Matrix4x4 Matrix2x2::embedInTwoQubit(const std::size_t qubitIndex) const { if (qubitIndex == 1) { return Matrix4x4::kron(Matrix2x2::identity(), *this); } - assert(false && "Invalid qubit index for single-qubit embed"); + llvm::reportFatalInternalError("Invalid qubit index for single-qubit embed"); } Matrix4x4 Matrix4x4::fromElements(const Complex& m00, const Complex& m01, @@ -452,12 +460,14 @@ Matrix4x4 Matrix4x4::kron(const Matrix2x2& lhs, const Matrix2x2& rhs) { std::array Matrix4x4::column(const std::size_t col) const { + assert(col < K_COLS); return {data[col], data[K_COLS + col], data[(2 * K_COLS) + col], data[(3 * K_COLS) + col]}; } void Matrix4x4::setColumn(const std::size_t col, const ArrayRef values) { + assert(col < K_COLS); assert(values.size() == K_ROWS && "setColumn requires exactly K_ROWS entries"); for (std::size_t row = 0; row < K_ROWS; ++row) { @@ -467,7 +477,7 @@ void Matrix4x4::setColumn(const std::size_t col, ArrayRef Matrix4x4::row(const std::size_t row) const { assert(row < K_ROWS); - return ArrayRef(data).slice(row * K_COLS, K_COLS); + return ArrayRef(data).slice(row * K_COLS, K_COLS); } void Matrix4x4::setRow(const std::size_t row, const ArrayRef values) { @@ -510,7 +520,7 @@ DynamicMatrix Matrix4x4::embedInNqubit(const std::size_t numQubits, if (numQubits == 2) { return DynamicMatrix(reorderForQubits(q0Index, q1Index)); } - const auto dim = static_cast(1ULL << numQubits); + const auto dim = checkedHilbertDim(numQubits); DynamicMatrix out(dim); const auto udim = static_cast(dim); for (std::size_t row = 0; row < udim; ++row) { @@ -541,7 +551,7 @@ Matrix4x4 Matrix4x4::reorderForQubits(const std::size_t q0Index, return fromElements(m[0], m[2], m[1], m[3], m[8], m[10], m[9], m[11], m[4], m[6], m[5], m[7], m[12], m[14], m[13], m[15]); } - assert(false && "Invalid qubit indices for two-qubit reorder"); + llvm::reportFatalInternalError("Invalid qubit indices for two-qubit reorder"); } SymmetricEigen4 Matrix4x4::symmetricEigen4() const { From e9c6141feececa03079e0453670f2bdb5934ff39 Mon Sep 17 00:00:00 2001 From: simon1hofmann <119581649+simon1hofmann@users.noreply.github.com> Date: Mon, 22 Jun 2026 15:18:56 +0200 Subject: [PATCH 08/10] Update mlir/lib/Dialect/QCO/Utils/Matrix.cpp Co-authored-by: Lukas Burgholzer Signed-off-by: simon1hofmann <119581649+simon1hofmann@users.noreply.github.com> --- mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 00b590271b..8960aa4dca 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -805,7 +805,7 @@ Matrix4x4::symmetricEigen4(const std::array& symmetric) { for (std::size_t col = 0; col < n; ++col) { for (std::size_t row = 0; row < n; ++row) { - result.eigenvectors(row, col) = Complex{z[row + (col * n)], 0.0}; + result.eigenvectors(row, col) = z[row + (col * n); } } return result; From 88532b9cc2352e6c71f590b0c395b7e0199326d8 Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 15:26:47 +0200 Subject: [PATCH 09/10] =?UTF-8?q?=F0=9F=8E=A8=20Unroll=20more=20loops?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 8960aa4dca..d722981dcc 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -443,9 +443,10 @@ Matrix4x4 Matrix4x4::fromDiagonal(const ArrayRef diagonalEntries) { assert(diagonalEntries.size() == K_ROWS && "fromDiagonal requires exactly K_ROWS entries"); Matrix4x4 out{}; - for (std::size_t i = 0; i < K_ROWS; ++i) { - out.data[(i * K_COLS) + i] = diagonalEntries[i]; - } + out.data[0] = diagonalEntries[0]; + out.data[5] = diagonalEntries[1]; + out.data[10] = diagonalEntries[2]; + out.data[15] = diagonalEntries[3]; return out; } @@ -470,9 +471,10 @@ void Matrix4x4::setColumn(const std::size_t col, assert(col < K_COLS); assert(values.size() == K_ROWS && "setColumn requires exactly K_ROWS entries"); - for (std::size_t row = 0; row < K_ROWS; ++row) { - data[(row * K_COLS) + col] = values[row]; - } + data[col] = values[0]; + data[K_COLS + col] = values[1]; + data[(2 * K_COLS) + col] = values[2]; + data[(3 * K_COLS) + col] = values[3]; } ArrayRef Matrix4x4::row(const std::size_t row) const { @@ -483,9 +485,11 @@ ArrayRef Matrix4x4::row(const std::size_t row) const { void Matrix4x4::setRow(const std::size_t row, const ArrayRef values) { assert(row < K_ROWS); assert(values.size() == K_COLS && "setRow requires exactly K_COLS entries"); - for (std::size_t col = 0; col < K_COLS; ++col) { - data[(row * K_COLS) + col] = values[col]; - } + const std::size_t rowBase = row * K_COLS; + data[rowBase + 0] = values[0]; + data[rowBase + 1] = values[1]; + data[rowBase + 2] = values[2]; + data[rowBase + 3] = values[3]; } std::array From 585646df4ac64bb4e6549b6f2b0e40fc8eaf3cfb Mon Sep 17 00:00:00 2001 From: Simon Hofmann Date: Mon, 22 Jun 2026 15:31:51 +0200 Subject: [PATCH 10/10] =?UTF-8?q?=F0=9F=9A=A8=20Fix=20syntax=20error=20in?= =?UTF-8?q?=20eigenvector=20assignment=20in=20Matrix.cpp?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- mlir/lib/Dialect/QCO/Utils/Matrix.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index d722981dcc..50d524daa1 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -809,7 +809,7 @@ Matrix4x4::symmetricEigen4(const std::array& symmetric) { for (std::size_t col = 0; col < n; ++col) { for (std::size_t row = 0; row < n; ++row) { - result.eigenvectors(row, col) = z[row + (col * n); + result.eigenvectors(row, col) = z[row + (col * n)]; } } return result;