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 diff --git a/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h b/mlir/include/mlir/Dialect/QCO/Utils/Matrix.h index 8f49d372d7..cdea703ed4 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. @@ -183,6 +188,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 +206,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. * @@ -215,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; }; /** @@ -322,6 +364,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 +382,81 @@ 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)`; must have + * length `K_ROWS`. + * @return Diagonal matrix with the given entries. + */ + [[nodiscard]] static Matrix4x4 + 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. + * @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; 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 setRow(std::size_t row, ArrayRef 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. * @@ -354,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); }; /** @@ -539,6 +714,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 +768,31 @@ inline constexpr bool std::disjunction_v, std::is_same, std::is_same, std::is_same>; + +/// 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 { + /// Eigenvalues in ascending order. + std::array eigenvalues{}; + /// Orthonormal eigenvectors as columns (column `j` matches `eigenvalues[j]`). + Matrix4x4 eigenvectors{}; +}; + } // namespace mlir::qco diff --git a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp index 8df45840a3..50d524daa1 100644 --- a/mlir/lib/Dialect/QCO/Utils/Matrix.cpp +++ b/mlir/lib/Dialect/QCO/Utils/Matrix.cpp @@ -11,6 +11,8 @@ #include "mlir/Dialect/QCO/Utils/Matrix.h" #include +#include +#include #include #include @@ -26,14 +28,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). @@ -81,10 +89,12 @@ 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]; @@ -92,10 +102,12 @@ multiply2x2(const std::array& lhs, } /// 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) { +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,27 @@ multiply4x4(const std::array& lhs, } } +/// Returns true if @p data is approximately the @p dim x @p dim identity +/// matrix. +[[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 (!entryIsApprox(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"); @@ -122,6 +155,14 @@ multiply4x4(const std::array& 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 && @@ -147,10 +188,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}; } @@ -166,7 +238,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) { @@ -208,10 +280,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; @@ -231,12 +300,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, K_ROWS, tol); +} + bool Matrix2x2::isApprox(const Matrix2x2& other, const double tol) const { return entriesAreApprox(data, other.data, tol); } @@ -245,6 +322,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 = checkedHilbertDim(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); + } + llvm::reportFatalInternalError("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, @@ -272,10 +383,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; @@ -296,6 +404,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 +431,83 @@ Complex Matrix4x4::determinant() const { data[12], data[13], data[14]); } +bool Matrix4x4::isIdentity(const double tol) const { + return isIdentityEntries(data, K_ROWS, tol); +} + +std::array Matrix4x4::diagonal() const { + return {data[0], data[5], data[10], data[15]}; +} + +Matrix4x4 Matrix4x4::fromDiagonal(const ArrayRef diagonalEntries) { + assert(diagonalEntries.size() == K_ROWS && + "fromDiagonal requires exactly K_ROWS entries"); + Matrix4x4 out{}; + out.data[0] = diagonalEntries[0]; + out.data[5] = diagonalEntries[1]; + out.data[10] = diagonalEntries[2]; + out.data[15] = diagonalEntries[3]; + 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 { + 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"); + 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 { + 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"); + 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 +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); } @@ -325,6 +516,310 @@ bool Matrix4x4::assignFrom(const DynamicMatrix& src) { return assignFromDynamicImpl(src, data); } +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 = checkedHilbertDim(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; +} + +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]); + } + llvm::reportFatalInternalError("Invalid qubit indices for two-qubit reorder"); +} + +SymmetricEigen4 Matrix4x4::symmetricEigen4() const { + return symmetricEigen4(realPart()); +} + +// 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; +} + +/// 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 = 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) { + 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 = std::hypot(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 +Matrix4x4::symmetricEigen4(const std::array& symmetric) { + constexpr std::size_t n = 4; + + SymmetricEigen4 result; + std::array z{}; + std::array subdiag{}; + symmetricTred24(symmetric, z, result.eigenvalues, subdiag); + symmetricTql24(result.eigenvalues, subdiag, z); + + 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)]; + } + } + return result; +} + +struct DynamicMatrix::Impl { + std::int64_t dim = 0; + SmallVector data; +}; + DynamicMatrix::DynamicMatrix() : impl_(std::make_unique()) {} DynamicMatrix::DynamicMatrix(const std::int64_t dim) @@ -433,7 +928,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 { @@ -453,4 +948,71 @@ 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 { + 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); +} + +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 afa0792415..c23426918b 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,251 @@ 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()); + 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) { + 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, 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}); + 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 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); + 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 = 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 = 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 + 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 = Matrix4x4::kron(x, Matrix2x2::identity()); + const Matrix4x4 onLow = Matrix4x4::kron(Matrix2x2::identity(), x); - DynamicMatrix as2x2; - as2x2.assignFrom(x); - EXPECT_TRUE(as2x2.isApprox(x)); - EXPECT_FALSE(as2x2.isApprox(Matrix2x2::identity())); + EXPECT_TRUE(onHigh.reorderForQubits(0, 1).isApprox(onHigh)); + EXPECT_TRUE(onHigh.reorderForQubits(1, 0).isApprox(onLow)); + EXPECT_TRUE(onLow.reorderForQubits(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(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)); +} - 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 = x.embedInNqubit(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 = x.embedInNqubit(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(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); + 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 = 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); + 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 = 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); + } + 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 = Matrix4x4::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 = Matrix4x4::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()); }