// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// eigenvalues.hpp
//
// Eigenvalue problem algorithms
//
// This file implements algorithms for computing the eigenvalues and
// eigenvectors of matrices.
//
// Main features:
// - eigen()                   — eigenvalues/eigenvectors of a general matrix (QR method, O(n^3))
// - eigen_symmetric()         — for symmetric matrices (Householder + QR, O(n^3), real eigenvalues)
// - eigen_generalized_symmetric() — generalized eigenvalue problem Av=λBv (eigen_symmetric after Cholesky transform)
// - power_iteration()         — power method for the largest eigenvalue (for large sparse matrices, O(iter * n^2))
// - inverse_iteration()       — inverse power method near a specified eigenvalue
//
// References:
//   Golub, Van Loan "Matrix Computations" 4th ed. Ch.8 (QR method), Ch.7 (symmetric eigenvalues)
//   Trefethen, Bau "Numerical Linear Algebra" Lec.24-29
//
// Note: eigenvalues of large sparse matrices require ARPACK-equivalent Lanczos/Arnoldi methods
//       (currently not implemented; consider integration with the accel/ module if needed)

#ifndef SANGI_EIGENVALUES_HPP
#define SANGI_EIGENVALUES_HPP

#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/core/Complex.hpp>
#include <math/linalg/decomposition.hpp>
#include <math/linalg/balance.hpp>
#include <cassert>

#include <algorithm>
#include <limits>
#include <random>
#include <utility>
#include <vector>

// Forward declarations (mutual references within eigenvalues.hpp)
namespace sangi { namespace algorithms {
    template<typename MA>
        requires sangi::BaseMatrixLike<MA>
    std::pair<Vector<sangi::element_t<MA>>, Matrix<sangi::element_t<MA>>> eigen_symmetric(const MA& a);
    template<typename MA, typename MB>
        requires sangi::BaseMatrixLike<MA>
              && sangi::BaseMatrixLike<MB>
              && std::same_as<sangi::element_t<MA>, sangi::element_t<MB>>
    std::pair<Vector<sangi::element_t<MA>>, Matrix<sangi::element_t<MA>>>
    eigen_generalized_symmetric(const MA& a, const MB& b);
}}

namespace sangi {
    namespace algorithms {

        //-----------------------------------------------------------------------------
        // dgebal-aware decomposition: convenient wrapper bundling balance + Hessenberg / Schur
        //
        //   The existing hessenberg_decomposition / schur_decomposition return-value
        //   contract guarantees that Q is an orthogonal matrix, so they do not apply
        //   balance internally (LAPACK convention).
        //   When you want to perform eigenvalue analysis with balance on an
        //   ill-conditioned matrix, the following wrappers give you balance +
        //   decomposition + inverse-transform information bundled together.
        //
        //   The returned Q is the orthogonal Schur vectors / Hessenberg transform
        //   matrix of the balanced matrix.
        //   To map back to the invariant subspace of the original A, use info to call
        //   rebalance_invariant_subspace (orthogonality is lost but the invariant subspace is preserved).
        //   The eigenvalues can be taken directly from the diagonal blocks of H / T,
        //   and since balance is a similarity transform they match the eigenvalues of
        //   the original A.
        //-----------------------------------------------------------------------------

        /// Return value of hessenberg_decomposition_balanced
        template<typename T>
        struct BalancedHessenberg {
            Matrix<T> H;         ///< upper Hessenberg form of the balanced matrix
            Matrix<T> Q;         ///< orthogonal transform of the balanced matrix (A_bal = Q H Q^T)
            BalanceInfo<T> info; ///< inverse balance information (for rebalance_invariant_subspace)
        };

        /// Return value of schur_decomposition_balanced
        template<typename T>
        struct BalancedSchur {
            Matrix<T> T_sch;     ///< real Schur form of the balanced matrix (quasi-upper-triangular)
            Matrix<T> Q;         ///< orthogonal Schur vectors of the balanced matrix (A_bal = Q T Q^T)
            BalanceInfo<T> info; ///< inverse balance information
        };

        /**
         * @brief balance + Hessenberg decomposition
         *
         * @param A in square matrix
         * @return BalancedHessenberg<T>
         *
         * @note The eigenvalues of A can be extracted from the diagonal blocks of H
         *       (after Wilkinson shift QR); since balance is a similarity transform they are identical to those of the original A.
         * @note To convert to the eigenvectors of A, map back to the invariant
         *       subspace with rebalance_invariant_subspace(Q, info) (orthogonality is lost).
         */
        template<typename T>
        BalancedHessenberg<T> hessenberg_decomposition_balanced(const BaseMatrix<T>& A) {
            if (A.rows() != A.cols()) {
                throw DimensionError("hessenberg_decomposition_balanced: matrix must be square");
            }
            Matrix<T> A_bal(A);
            auto info = balance_matrix(A_bal);
            auto [H, Q] = hessenberg_decomposition(A_bal);
            return BalancedHessenberg<T>{std::move(H), std::move(Q), std::move(info)};
        }

        /**
         * @brief balance + real Schur decomposition
         *
         * @param A in square matrix
         * @return BalancedSchur<T>
         *
         * @note A_bal = Q T Q^T holds (Q is orthogonal); however the relation to the
         *       original A is A = (P D Q) T (P D Q)^{-1}, where Q_full = P D Q is not orthogonal.
         * @note For LQR Laub use cases, map the selected Schur vector columns back to
         *       the invariant subspace with rebalance_invariant_subspace(Q, info) (orthogonality is lost).
         */
        template<typename T>
        BalancedSchur<T> schur_decomposition_balanced(const BaseMatrix<T>& A) {
            if (A.rows() != A.cols()) {
                throw DimensionError("schur_decomposition_balanced: matrix must be square");
            }
            Matrix<T> A_bal(A);
            auto info = balance_matrix(A_bal);
            auto [T_sch, Q] = schur_decomposition(A_bal);
            return BalancedSchur<T>{std::move(T_sch), std::move(Q), std::move(info)};
        }

        //-----------------------------------------------------------------------------
        // General eigenvalue computation functions
        //-----------------------------------------------------------------------------

        // Computes the eigenvalues and eigenvectors of a matrix
        // Eigenvalues are returned as complex numbers (even for real symmetric matrices)
        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
        std::pair<Vector<Complex<sangi::element_t<MA>>>, Matrix<Complex<sangi::element_t<MA>>>>
            eigen(const MA& a_basemat) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);

            Matrix<T> a(a_basemat);  // 1 copy: body modifies a (Hessenberg, QR)
            const auto n = a.rows();

            if (n != a.cols()) {
                assert(false && "DimensionError: Eigenvalue calculation requires a square matrix");
                throw DimensionError("Eigenvalue calculation requires a square matrix");
            }

            // Check symmetry - unified function name (isSymmetric → is_symmetric)
            bool is_symmetric = true;
            for (std::size_t i = 0; i < n && is_symmetric; ++i) {
                for (std::size_t j = i + 1; j < n && is_symmetric; ++j) {
                    if (!approximately_equal(a(i, j), a(j, i))) {
                        is_symmetric = false;
                    }
                }
            }

            // For a real symmetric matrix, use the dedicated symmetric-matrix function
            if (is_symmetric) {
                // The eigenvalues of a real symmetric matrix are real, so conversion is needed
                auto [eigenvalues_real, eigenvectors_real] = eigen_symmetric(a);

                // Conversion from real values
                Vector<Complex<T>> eigenvalues(n);
                for (std::size_t i = 0; i < n; ++i) {
                    eigenvalues[i] = Complex<T>(eigenvalues_real[i], T(0));
                }

                // Conversion from real matrix to complex matrix
                Matrix<Complex<T>> eigenvectors(n, n);
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        eigenvectors(i, j) = Complex<T>(eigenvectors_real(i, j), T(0));
                    }
                }

                return { eigenvalues, eigenvectors };
            }

            // For a non-symmetric matrix

            // Computation in complex numbers is required
            Vector<Complex<T>> eigenvalues(n);
            Matrix<Complex<T>> eigenvectors(n, n);

#if SANGI_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                Matrix<T> Acopy = a;
                std::vector<T> wr(n), wi(n);
                Matrix<T> vr(n, n);  // right eigenvectors (real form)
                lapack_int info;
                if constexpr (std::is_same_v<T, float>)
                    info = LAPACKE_sgeev(LAPACK_ROW_MAJOR, 'N', 'V',
                                         (lapack_int)n, Acopy.data(), (lapack_int)n,
                                         wr.data(), wi.data(),
                                         nullptr, (lapack_int)n,
                                         vr.data(), (lapack_int)n);
                else
                    info = LAPACKE_dgeev(LAPACK_ROW_MAJOR, 'N', 'V',
                                         (lapack_int)n, Acopy.data(), (lapack_int)n,
                                         wr.data(), wi.data(),
                                         nullptr, (lapack_int)n,
                                         vr.data(), (lapack_int)n);
                if (info > 0)
                    throw MathError("eigen: LAPACKE_geev did not converge");

                // Convert the LAPACKE output to complex form
                for (std::size_t i = 0; i < n; ++i)
                    eigenvalues[i] = Complex<T>(wr[i], wi[i]);

                // Eigenvector conversion: complex conjugate pairs are stored as real/imaginary parts in adjacent columns
                for (std::size_t j = 0; j < n; ++j) {
                    if (wi[j] == T(0)) {
                        // real eigenvalue → real eigenvector
                        for (std::size_t i = 0; i < n; ++i)
                            eigenvectors(i, j) = Complex<T>(vr(i, j), T(0));
                    } else if (j + 1 < n && wi[j] > T(0)) {
                        // complex conjugate pair: v(:,j) ± i*v(:,j+1)
                        for (std::size_t i = 0; i < n; ++i) {
                            eigenvectors(i, j)     = Complex<T>(vr(i, j),  vr(i, j + 1));
                            eigenvectors(i, j + 1) = Complex<T>(vr(i, j), -vr(i, j + 1));
                        }
                        ++j;  // skip the next column
                    }
                }
                return { eigenvalues, eigenvectors };
            }
#endif

            // dgebal-equivalent balancing: applies A → P^T D^{-1} A D P in-place.
            // Eigenvalues are preserved by the similarity transform, but the matrix
            // norms become balanced, improving the numerical stability of the Schur
            // decomposition and the subsequent inverse-iteration LU. Well-conditioned
            // cases (= symmetric matrix norms, no isolated eigenvalues) are effectively a no-op (= regression 0).
            // The eigenvectors are mapped back to the original A at the end with unbalance_eigenvectors.
            // Reference: LAPACK dgeev calls dgebal internally. This is the equivalent treatment for the non-MKL fallback.
            auto bal_info = balance_matrix(a);

            // Real Schur decomposition based: A_balanced = Q T Q^T
            // T is quasi-upper-triangular (1×1 real eigenvalues + 2×2 complex conjugate pair blocks)
            auto [T_sch, Q_sch] = schur_decomposition(a);
            const T eps = numeric_traits<T>::epsilon();

            // Extract eigenvalues from the quasi-upper-triangular matrix T
            std::size_t idx = 0;
            while (idx < n) {
                if (idx + 1 < n && std::abs(T_sch(idx + 1, idx)) > eps *
                    (std::abs(T_sch(idx, idx)) + std::abs(T_sch(idx + 1, idx + 1)) + T(1))) {
                    // 2×2 block → complex conjugate pair
                    T a11 = T_sch(idx, idx), a12 = T_sch(idx, idx + 1);
                    T a21 = T_sch(idx + 1, idx), a22 = T_sch(idx + 1, idx + 1);
                    T tr = a11 + a22;
                    T det = a11 * a22 - a12 * a21;
                    T disc = tr * tr - T(4) * det;
                    if (disc < T(0)) {
                        T re = tr / T(2);
                        T im = std::sqrt(-disc) / T(2);
                        eigenvalues[idx]     = Complex<T>(re,  im);
                        eigenvalues[idx + 1] = Complex<T>(re, -im);
                    } else {
                        T sqd = std::sqrt(disc);
                        eigenvalues[idx]     = Complex<T>((tr + sqd) / T(2), T(0));
                        eigenvalues[idx + 1] = Complex<T>((tr - sqd) / T(2), T(0));
                    }
                    idx += 2;
                } else {
                    // 1×1 block → real eigenvalue
                    eigenvalues[idx] = Complex<T>(T_sch(idx, idx), T(0));
                    idx += 1;
                }
            }

            // Eigenvectors: inverse iteration
            // Repeatedly solve (A - λI)x = b to find the eigenvector for λ
            for (std::size_t k = 0; k < n; ++k) {
                Complex<T> lambda = eigenvalues[k];

                // For the second of a complex conjugate pair, set the conjugate
                if (k > 0 && std::abs(lambda.imag()) > eps &&
                    sangi::abs(lambda - sangi::conj(eigenvalues[k - 1])) < eps * T(100)) {
                    for (std::size_t i = 0; i < n; ++i)
                        eigenvectors(i, k) = sangi::conj(eigenvectors(i, k - 1));
                    continue;
                }

                // Build (A - λI) as a complex matrix
                Matrix<Complex<T>> B(n, n);
                for (std::size_t i = 0; i < n; ++i)
                    for (std::size_t j = 0; j < n; ++j)
                        B(i, j) = Complex<T>(a(i, j), T(0))
                                  - (i == j ? lambda : Complex<T>(T(0)));

                // Shift by a tiny perturbation (avoid singularity)
                for (std::size_t i = 0; i < n; ++i)
                    B(i, i) += Complex<T>(eps * T(10) * sangi::abs(lambda) + eps, T(0));

                // LU decomposition (partial pivoting)
                std::vector<std::size_t> piv(n);
                for (std::size_t i = 0; i < n; ++i) piv[i] = i;
                Matrix<Complex<T>> LU = B;
                for (std::size_t col = 0; col < n; ++col) {
                    // pivot selection
                    T maxval = T(0);
                    std::size_t maxrow = col;
                    for (std::size_t row = col; row < n; ++row) {
                        T val = sangi::abs(LU(row, col));
                        if (val > maxval) { maxval = val; maxrow = row; }
                    }
                    if (maxrow != col) {
                        std::swap(piv[col], piv[maxrow]);
                        for (std::size_t j = 0; j < n; ++j)
                            std::swap(LU(col, j), LU(maxrow, j));
                    }
                    if (sangi::abs(LU(col, col)) < eps * eps)
                        LU(col, col) = Complex<T>(eps * eps, T(0));
                    for (std::size_t row = col + 1; row < n; ++row) {
                        LU(row, col) /= LU(col, col);
                        for (std::size_t j = col + 1; j < n; ++j)
                            LU(row, j) -= LU(row, col) * LU(col, j);
                    }
                }

                // Inverse iteration: repeat x ← solve(B, x) 3 times
                Vector<Complex<T>> x(n);
                for (std::size_t i = 0; i < n; ++i)
                    x[i] = Complex<T>(T(1), T(0));

                for (int iter = 0; iter < 3; ++iter) {
                    // Build Pb
                    Vector<Complex<T>> b(n);
                    for (std::size_t i = 0; i < n; ++i) b[i] = x[piv[i]];
                    // forward substitution
                    for (std::size_t i = 1; i < n; ++i)
                        for (std::size_t j = 0; j < i; ++j)
                            b[i] -= LU(i, j) * b[j];
                    // backward substitution
                    for (int i = static_cast<int>(n) - 1; i >= 0; --i) {
                        for (std::size_t j = i + 1; j < n; ++j)
                            b[i] -= LU(i, j) * b[j];
                        b[i] /= LU(i, i);
                    }
                    x = b;
                    // normalization
                    T nrm = T(0);
                    for (std::size_t i = 0; i < n; ++i) nrm += sangi::normSq(x[i]);
                    nrm = std::sqrt(nrm);
                    if (nrm > T(0))
                        for (std::size_t i = 0; i < n; ++i) x[i] /= nrm;
                }

                for (std::size_t i = 0; i < n; ++i)
                    eigenvectors(i, k) = x[i];
            }

            // Map the eigenvectors of the balanced matrix back to the eigenvectors of the original A
            // (v = P D v_balanced, inverse scaling + permutation transform, dgebak equivalent)
            unbalance_eigenvectors(eigenvectors, bal_info);

            return { eigenvalues, eigenvectors };
        }

        // Computes the eigenvalues and eigenvectors of a real symmetric matrix
        // All eigenvalues of a real symmetric matrix are real
        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
        std::pair<Vector<sangi::element_t<MA>>, Matrix<sangi::element_t<MA>>>
            eigen_symmetric(const MA& a_basemat) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);

            Matrix<T> a(a_basemat);  // 1 copy: body modifies a (Jacobi)
            const auto n = a.rows();

            if (n != a.cols()) {
                assert(false && "DimensionError: Symmetric eigenvalue calculation requires a square matrix");
                throw DimensionError("Symmetric eigenvalue calculation requires a square matrix");
            }

            // Check symmetry
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = i + 1; j < n; ++j) {
                    if (!approximately_equal(a(i, j), a(j, i))) {
                        throw MathError("Matrix must be symmetric for eigen_symmetric");
                    }
                }
            }

            // Arrays storing the eigenvalues and eigenvectors
            Vector<T> eigenvalues(n);
            Matrix<T> eigenvectors(n, n);

#if SANGI_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                lapack_int info;
                // For large matrices, MRRR (Multiple Relatively Robust Representations,
                // syevr via DSTEMR) is faster than QR-based syev. Like syev, the
                // eigenvalues are returned in ascending order and the eigenvectors as
                // columns. syevr destroys the input matrix and writes the eigenvectors
                // to a separate buffer z, so its handling differs from syev's in-place.
                constexpr std::size_t SYEVR_CUTOFF = 1000;
                if (n >= SYEVR_CUTOFF) {
                    Matrix<T> acopy = a;            // working copy that gets destroyed
                    lapack_int m_found = 0;
                    std::vector<lapack_int> isuppz(2 * (n > 0 ? n : 1));
                    if constexpr (std::is_same_v<T, float>)
                        info = LAPACKE_ssyevr(LAPACK_ROW_MAJOR, 'V', 'A', 'U',
                                              (lapack_int)n, acopy.data(), (lapack_int)n,
                                              0.0f, 0.0f, 0, 0, 0.0f, &m_found,
                                              eigenvalues.data(), eigenvectors.data(),
                                              (lapack_int)n, isuppz.data());
                    else
                        info = LAPACKE_dsyevr(LAPACK_ROW_MAJOR, 'V', 'A', 'U',
                                              (lapack_int)n, acopy.data(), (lapack_int)n,
                                              0.0, 0.0, 0, 0, 0.0, &m_found,
                                              eigenvalues.data(), eigenvectors.data(),
                                              (lapack_int)n, isuppz.data());
                    if (info != 0)
                        throw MathError("eigen_symmetric: LAPACKE_syevr did not converge");
                    return { eigenvalues, eigenvectors };
                }

                // LAPACKE_syev: computes all eigenvalues and eigenvectors of a symmetric matrix
                eigenvectors = a;  // gets overwritten
                if constexpr (std::is_same_v<T, float>)
                    info = LAPACKE_ssyev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n, eigenvectors.data(), (lapack_int)n,
                                         eigenvalues.data());
                else
                    info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n, eigenvectors.data(), (lapack_int)n,
                                         eigenvalues.data());
                if (info > 0)
                    throw MathError("eigen_symmetric: LAPACKE_syev did not converge");
                // LAPACKE_syev stores eigenvalues in ascending order and eigenvectors as columns (same format as the existing code)
                return { eigenvalues, eigenvectors };
            }
#endif

            // ================================================================
            // Block tridiagonalization + implicit QL iteration + direct inverse transform
            // (LAPACK DSYEV approach, optimized version that does not form Q explicitly)
            // ================================================================
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                // --- Phase 1: block Householder tridiagonalization (compact form) ---
                auto tri = tridiagonalize_compact(a);
                const std::size_t nm1 = n > 1 ? n - 1 : 0;
                std::vector<T> d_arr(n), e_arr(n, T(0));
                for (std::size_t i = 0; i < n; ++i) d_arr[i] = tri.d[i];
                for (std::size_t i = 0; i < nm1; ++i) e_arr[i] = tri.e[i];

                // --- Phase 2: implicit QL iteration (column-major Z for better cache efficiency) ---
                // Store Z column-major: z_cm[j*n + k] = Z(k, j)
                std::vector<T> z_cm(n * n, T(0));
                for (std::size_t i = 0; i < n; ++i)
                    z_cm[i * n + i] = T(1);  // Z = I (column-major)

                const int max_iter = 30 * static_cast<int>(n);
                const T eps = numeric_traits<T>::epsilon();
                int total_iter = 0;

                for (std::size_t l = 0; l < n; ) {
                    std::size_t m = l;
                    while (m + 1 < n) {
                        T tst = std::abs(d_arr[m]) + std::abs(d_arr[m + 1]);
                        if (std::abs(e_arr[m]) <= eps * tst)
                            break;
                        ++m;
                    }

                    if (m == l) { ++l; continue; }

                    if (++total_iter > max_iter)
                        throw MathError("eigen_symmetric: QL iteration did not converge");

                    T g = (d_arr[l + 1] - d_arr[l]) / (T(2) * e_arr[l]);
                    T r = std::sqrt(g * g + T(1));
                    T sign_g = (g >= T(0)) ? T(1) : T(-1);
                    g = d_arr[m] - d_arr[l] + e_arr[l] / (g + sign_g * r);

                    T s = T(1), c = T(1), p = T(0);

                    bool lucky = false;
                    for (std::size_t ii = 0; ii < m - l; ++ii) {
                        std::size_t i = m - 1 - ii;
                        T f = s * e_arr[i];
                        T b = c * e_arr[i];

                        if (std::abs(f) >= std::abs(g)) {
                            c = g / f;
                            r = std::sqrt(c * c + T(1));
                            e_arr[i + 1] = f * r;
                            s = T(1) / r;
                            c *= s;
                        } else {
                            s = f / g;
                            r = std::sqrt(s * s + T(1));
                            e_arr[i + 1] = g * r;
                            c = T(1) / r;
                            s *= c;
                        }

                        if (e_arr[i + 1] == T(0)) {
                            d_arr[i + 1] -= p;
                            e_arr[m] = T(0);
                            lucky = true;
                            break;
                        }

                        g = d_arr[i + 1] - p;
                        r = (d_arr[i] - g) * s + T(2) * c * b;
                        p = s * r;
                        d_arr[i + 1] = g + p;
                        g = c * r - b;

                        // Column-major Z: columns i, i+1 are contiguous memory → cache-friendly
                        T* __restrict col_i  = z_cm.data() + i * n;
                        T* __restrict col_i1 = z_cm.data() + (i + 1) * n;
                        for (std::size_t k = 0; k < n; ++k) {
                            T zi  = col_i[k];
                            T zi1 = col_i1[k];
                            col_i1[k] = s * zi + c * zi1;
                            col_i[k]  = c * zi - s * zi1;
                        }
                    }

                    if (!lucky) {
                        d_arr[l] -= p;
                        e_arr[l] = g;
                        e_arr[m] = T(0);
                    }
                }

                // --- Phase 3: direct Householder inverse transform ---
                // Z ← Q * Z (apply the Householder columns directly to Z without forming Q)
                apply_householder_sequence<T>(
                    tri.refs, n, z_cm.data(), n);

                // Convert column-major Z → row-major eigenvectors + sort ascending
                std::vector<std::size_t> idx(n);
                for (std::size_t i = 0; i < n; ++i) { eigenvalues[i] = d_arr[i]; idx[i] = i; }
                std::sort(idx.begin(), idx.end(),
                          [&](std::size_t aa, std::size_t bb) { return eigenvalues[aa] < eigenvalues[bb]; });

                Vector<T> sorted_evals(n);
                Matrix<T> sorted_evecs(n, n);
                for (std::size_t i = 0; i < n; ++i) {
                    sorted_evals[i] = eigenvalues[idx[i]];
                    // Z column idx[i] (column-major) → sorted_evecs column i (row-major)
                    const T* src_col = z_cm.data() + idx[i] * n;
                    for (std::size_t j = 0; j < n; ++j)
                        sorted_evecs(j, i) = src_col[j];
                }
                return { std::move(sorted_evals), std::move(sorted_evecs) };
            }

            // Generic type: Schur-decomposition-based eigenvalue computation
            auto [Sch, Q_sch] = schur_decomposition(a);
            for (std::size_t i = 0; i < n; ++i) {
                eigenvalues[i] = Sch(i, i);
                for (std::size_t j = 0; j < n; ++j) {
                    eigenvectors(j, i) = Q_sch(j, i);
                }
            }
            // Sort eigenvalues in ascending order
            std::vector<std::size_t> idx(n);
            for (std::size_t i = 0; i < n; ++i) idx[i] = i;
            std::sort(idx.begin(), idx.end(),
                      [&](std::size_t aa, std::size_t bb) {
                          return eigenvalues[aa] < eigenvalues[bb];
                      });
            Vector<T> sorted_evals(n);
            Matrix<T> sorted_evecs(n, n);
            for (std::size_t i = 0; i < n; ++i) {
                sorted_evals[i] = eigenvalues[idx[i]];
                for (std::size_t j = 0; j < n; ++j) {
                    sorted_evecs(j, i) = eigenvectors(j, idx[i]);
                }
            }
            eigenvalues = sorted_evals;
            eigenvectors = sorted_evecs;

            return { eigenvalues, eigenvectors };
        }

        // ====================================================================
        // Symmetric eigenvalues with a convergence flag ─ returns results after *post-verification* (residual ‖A·v−λ·v‖)
        // ====================================================================
        //
        // Motivation (TODO A.1): eigen_symmetric throws an *exception* (MathError) on non-convergence, and on
        //   convergence merely returns the eigenpairs, without telling the caller "how accurate each eigenpair is".
        //   Here, without changing the existing API (non-destructive), we (1) catch the exception and convert it
        //   into a converged flag, and (2) actually compute the residual ‖A·vⱼ − λⱼ·vⱼ‖₂ of each returned eigenpair
        //   and return the maximum residual.
        //   ★ Because we verify the result rather than "trusting the solver", converged will not erroneously
        //     become true even if the solver is buggy. If you need exact (rational, proof-backed) verification,
        //     also use numeric_bridge::certifyEigenpair. For symmetric matrices (real eigenvalues / real eigenvectors ⇒ the residual is computed entirely in real arithmetic).
        template<typename T>
        struct SymEigenResult {
            Vector<T> values;                 // eigenvalues (ascending)
            Matrix<T> vectors;                // eigenvectors (per column)
            bool   converged   = false;       // solver success ∧ all residuals ≤ tol
            bool   solverOk    = false;       // whether the solver completed without exceptions
            T      maxResidual = T(0);        // max_j ‖A·vⱼ − λⱼ·vⱼ‖₂
        };

        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
        SymEigenResult<sangi::element_t<MA>>
        eigen_symmetric_checked(const MA& a,
                                sangi::element_t<MA> residualTol =
                                    std::numeric_limits<sangi::element_t<MA>>::epsilon() * sangi::element_t<MA>(1e4)) {
            using T = sangi::element_t<MA>;
            SymEigenResult<T> r;
            try {
                auto pr   = eigen_symmetric(a);
                r.values  = pr.first;
                r.vectors = pr.second;
                r.solverOk = true;
            } catch (...) {
                r.solverOk = false;           // non-convergence/non-symmetric → converged=false (fail-closed)
                return r;
            }
            Matrix<T> A(a);
            const std::size_t n = A.rows();
            using std::sqrt;
            T maxres = T(0);
            for (std::size_t j = 0; j < n; ++j) {           // each eigenpair (λⱼ, j-th column vector)
                T sumsq = T(0);
                for (std::size_t i = 0; i < n; ++i) {
                    T av = T(0);
                    for (std::size_t k = 0; k < n; ++k) av += A(i, k) * r.vectors(k, j);
                    const T resi = av - r.values[j] * r.vectors(i, j);
                    sumsq += resi * resi;
                }
                const T nrm = sqrt(sumsq);
                if (nrm > maxres) maxres = nrm;
            }
            r.maxResidual = maxres;
            r.converged   = r.solverOk && (n == 0 || maxres <= residualTol);
            return r;
        }

        // Computes the eigenvalues and eigenvectors of a Hermitian matrix (complex symmetric matrix)
        // All eigenvalues of a Hermitian matrix are real
        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
                  && requires { typename numeric_traits<sangi::element_t<MA>>::real_type; }
        std::pair<Vector<typename numeric_traits<sangi::element_t<MA>>::real_type>,
                  Matrix<sangi::element_t<MA>>>
            eigen_hermitian(const MA& a) {
            using EltT = sangi::element_t<MA>;                       // Complex<T>
            using T = typename numeric_traits<EltT>::real_type;      // T (real)
            SANGI_STATIC_ASSERT_SQUARE(MA);

            const auto n = a.rows();

            if (n != a.cols()) {
                assert(false && "DimensionError: Hermitian eigenvalue calculation requires a square matrix");
                throw DimensionError("Hermitian eigenvalue calculation requires a square matrix");
            }

            // Check Hermitian property (A^H = A)
            {
                T eps = numeric_traits<T>::epsilon() * T(100);
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = i + 1; j < n; ++j) {
                        if (sangi::abs(a(i, j) - sangi::conj(a(j, i))) > eps)
                            throw MathError("Matrix must be Hermitian for eigen_hermitian");
                    }
                }
            }

            // Arrays storing the eigenvalues and eigenvectors
            Vector<T> eigenvalues(n);
            Matrix<Complex<T>> eigenvectors(n, n);

#if SANGI_HAS_MKL
            if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double>) {
                eigenvectors = a;
                lapack_int info;
                if constexpr (std::is_same_v<T, float>)
                    info = LAPACKE_cheev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n,
                                         reinterpret_cast<lapack_complex_float*>(eigenvectors.data()),
                                         (lapack_int)n, eigenvalues.data());
                else
                    info = LAPACKE_zheev(LAPACK_ROW_MAJOR, 'V', 'U',
                                         (lapack_int)n,
                                         reinterpret_cast<lapack_complex_double*>(eigenvectors.data()),
                                         (lapack_int)n, eigenvalues.data());
                if (info > 0)
                    throw MathError("eigen_hermitian: LAPACKE_heev did not converge");
                return { eigenvalues, eigenvectors };
            }
#endif
            // Solve by converting the Hermitian matrix A = X + iY into a 2n×2n real symmetric matrix
            // A(u+iv) = λ(u+iv) ⇔ Xu-Yv = λu, Yu+Xv = λv
            // The eigenvalues of M = [[X, -Y], [Y, X]] are a double copy of the eigenvalues of A
            // (X is symmetric, Y is anti-symmetric, so M is symmetric)
            Matrix<T> M(2 * n, 2 * n, T(0));
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = 0; j < n; ++j) {
                    T re = a(i, j).real();
                    T im = a(i, j).imag();
                    M(i, j) = re;           // X
                    M(i, n + j) = -im;      // -Y
                    M(n + i, j) = im;       // Y
                    M(n + i, n + j) = re;   // X
                }
            }

            auto [evals_2n, evecs_2n] = eigen_symmetric(M);

            // Eigenvalues appear in duplicate, so take every other one (ascending)
            for (std::size_t i = 0; i < n; ++i)
                eigenvalues[i] = evals_2n[2 * i];

            // Convert eigenvectors to Complex
            // Construct the eigenvector u + iv of A from the eigenvector [u; v] of M
            for (std::size_t k = 0; k < n; ++k) {
                std::size_t k2 = 2 * k; // corresponding index on the 2n side
                T norm_sq = T(0);
                for (std::size_t i = 0; i < n; ++i) {
                    T re = evecs_2n(i, k2);
                    T im = evecs_2n(n + i, k2);
                    eigenvectors(i, k) = Complex<T>(re, im);
                    norm_sq += re * re + im * im;
                }
                // normalization
                T inv_norm = T(1) / std::sqrt(norm_sq);
                for (std::size_t i = 0; i < n; ++i)
                    eigenvectors(i, k) = eigenvectors(i, k) * Complex<T>(inv_norm);
            }

            return { eigenvalues, eigenvectors };
        }

        //-----------------------------------------------------------------------------
        // Generalized eigenvalue problem
        //-----------------------------------------------------------------------------

        // Solution of the generalized eigenvalue problem Ax = λBx
        // (A,B) is called a generalized eigenpair
        template<typename MA, typename MB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseMatrixLike<MB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<MB>>
        std::pair<Vector<Complex<sangi::element_t<MA>>>, Matrix<Complex<sangi::element_t<MA>>>>
            eigen_generalized(const MA& a_basemat, const MB& b_basemat) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_SQUARE(MB);

            Matrix<T> a(a_basemat);  // 1 copy: body modifies a
            Matrix<T> b(b_basemat);  // 1 copy: body modifies b
            const auto n = a.rows();

            if (n != a.cols() || b.rows() != n || b.cols() != n) {
                assert(false && "DimensionError: Generalized eigenvalue problem requires square matrices of same size");
                throw DimensionError("Generalized eigenvalue problem requires square matrices of same size");
            }

            // Check symmetry
            bool is_symmetric_a = true;
            bool is_symmetric_b = true;

            for (std::size_t i = 0; i < n && (is_symmetric_a || is_symmetric_b); ++i) {
                for (std::size_t j = i + 1; j < n && (is_symmetric_a || is_symmetric_b); ++j) {
                    if (is_symmetric_a && !approximately_equal(a(i, j), a(j, i))) {
                        is_symmetric_a = false;
                    }
                    if (is_symmetric_b && !approximately_equal(b(i, j), b(j, i))) {
                        is_symmetric_b = false;
                    }
                }
            }

            // If both are symmetric matrices, it is a symmetric generalized eigenvalue problem
            if (is_symmetric_a && is_symmetric_b) {
                auto [eigenvalues_real, eigenvectors_real] = eigen_generalized_symmetric(a, b);

                // Convert real values to complex
                Vector<Complex<T>> eigenvalues(n);
                Matrix<Complex<T>> eigenvectors(n, n);

                for (std::size_t i = 0; i < n; ++i) {
                    eigenvalues[i] = Complex<T>(eigenvalues_real[i], T(0));
                    for (std::size_t j = 0; j < n; ++j) {
                        eigenvectors(i, j) = Complex<T>(eigenvectors_real(i, j), T(0));
                    }
                }

                return { eigenvalues, eigenvectors };
            }

            // Non-symmetric generalized eigenvalue problem

            // Arrays storing the eigenvalues and eigenvectors
            Vector<Complex<T>> eigenvalues(n);
            Matrix<Complex<T>> eigenvectors(n, n);

            // Solution of the generalized eigenvalue problem via QZ decomposition
            // (The MKL LAPACKE path is planned for the future since Matrix::data() is not yet supported)
            auto [S, Tri, Qm, Zm] = qz_decomposition(a, b);

            const T eps_eig = numeric_traits<T>::epsilon();

            // Extract generalized eigenvalues from the block diagonal of (S, Tri)
            std::size_t i = 0;
            while (i < n) {
                if (i + 1 < n && std::abs(S(i + 1, i)) > eps_eig *
                    (std::abs(S(i, i)) + std::abs(S(i + 1, i + 1)) + T(1))) {
                    // 2×2 block → complex conjugate pair
                    T a11 = S(i, i), a12 = S(i, i + 1);
                    T a21 = S(i + 1, i), a22 = S(i + 1, i + 1);
                    T b11 = Tri(i, i), b12 = Tri(i, i + 1), b22 = Tri(i + 1, i + 1);

                    // 2×2 generalized eigenvalue: det(S22 - λ T22) = 0
                    // (a11-λb11)(a22-λb22) - (a12-λb12)*a21 = 0
                    // coefficients: b11*b22*λ² - (a11*b22+a22*b11-a21*b12)*λ + (a11*a22-a12*a21) = 0
                    T qa = b11 * b22;
                    T qb = -(a11 * b22 + a22 * b11 - a21 * b12);
                    T qc = a11 * a22 - a12 * a21;

                    if (std::abs(qa) < eps_eig) {
                        // degenerate case
                        eigenvalues[i] = Complex<T>(std::numeric_limits<T>::infinity(), T(0));
                        eigenvalues[i + 1] = Complex<T>(std::numeric_limits<T>::infinity(), T(0));
                    } else {
                        T disc = qb * qb - T(4) * qa * qc;
                        T real_part = -qb / (T(2) * qa);
                        if (disc < T(0)) {
                            T imag_part = std::sqrt(-disc) / (T(2) * qa);
                            eigenvalues[i] = Complex<T>(real_part, imag_part);
                            eigenvalues[i + 1] = Complex<T>(real_part, -imag_part);
                        } else {
                            T sq = std::sqrt(disc);
                            eigenvalues[i] = Complex<T>((-qb + sq) / (T(2) * qa), T(0));
                            eigenvalues[i + 1] = Complex<T>((-qb - sq) / (T(2) * qa), T(0));
                        }
                    }
                    i += 2;
                } else {
                    // 1×1 block
                    T alpha = S(i, i);
                    T beta_val = Tri(i, i);
                    if (std::abs(beta_val) < eps_eig) {
                        eigenvalues[i] = Complex<T>(
                            std::numeric_limits<T>::infinity(), T(0));
                    } else {
                        eigenvalues[i] = Complex<T>(alpha / beta_val, T(0));
                    }
                    ++i;
                }
            }

            // Eigenvectors: obtained approximately by inverse transform
            // Solve (S - λ_i T) z = 0 and transform to the original space with Zm * z
            for (std::size_t col = 0; col < n; ++col) {
                if (std::abs(eigenvalues[col].imag()) > eps_eig) {
                    // complex eigenvalue: approximately use the corresponding column of Z
                    for (std::size_t row = 0; row < n; ++row) {
                        eigenvectors(row, col) = Complex<T>(Zm(row, col), T(0));
                    }
                    continue;
                }
                // real eigenvalue: refine the eigenvector by inverse iteration
                T lam = eigenvalues[col].real();
                // Find an approximate solution of (A - λB)v = 0 by inverse iteration
                // initial vector: column col of Z
                Matrix<T> C(n, n);
                for (std::size_t r = 0; r < n; ++r) {
                    for (std::size_t c = 0; c < n; ++c) {
                        C(r, c) = a(r, c) - lam * b(r, c);
                    }
                }
                // Add a tiny perturbation to C (avoid singularity)
                for (std::size_t r = 0; r < n; ++r) {
                    C(r, r) += eps_eig * (std::abs(C(r, r)) + T(1));
                }
                // Inverse iteration (up to 5 steps)
                Vector<T> v(n);
                for (std::size_t r = 0; r < n; ++r) v[r] = Zm(r, col);
                try {
                    auto [lu_mat, pivots] = lu_decomposition(C);
                    for (int iter = 0; iter < 5; ++iter) {
                        // pivot permutation
                        Vector<T> rhs = v;
                        for (std::size_t ii = 0; ii < n; ++ii)
                            if (pivots[ii] != ii) std::swap(rhs[ii], rhs[pivots[ii]]);
                        // forward substitution
                        for (std::size_t ii = 1; ii < n; ++ii)
                            for (std::size_t jj = 0; jj < ii; ++jj)
                                rhs[ii] -= lu_mat(ii, jj) * rhs[jj];
                        // backward substitution
                        for (std::size_t ii = n; ii-- > 0; ) {
                            for (std::size_t jj = ii + 1; jj < n; ++jj)
                                rhs[ii] -= lu_mat(ii, jj) * rhs[jj];
                            rhs[ii] /= lu_mat(ii, ii);
                        }
                        v = rhs;
                        // normalization
                        T nrm = T(0);
                        for (std::size_t r = 0; r < n; ++r) nrm += v[r] * v[r];
                        nrm = std::sqrt(nrm);
                        if (nrm > eps_eig)
                            for (std::size_t r = 0; r < n; ++r) v[r] /= nrm;
                    }
                    for (std::size_t r = 0; r < n; ++r)
                        eigenvectors(r, col) = Complex<T>(v[r], T(0));
                } catch (...) {
                    for (std::size_t r = 0; r < n; ++r) {
                        eigenvectors(r, col) = Complex<T>(Zm(r, col), T(0));
                    }
                }
            }

            return { eigenvalues, eigenvectors };
        }

        // Solution of the symmetric generalized eigenvalue problem Ax = λBx
        // (A,B) is a pair of symmetric matrices, with B positive definite
        template<typename MA, typename MB>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseMatrixLike<MB>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<MB>>
        std::pair<Vector<sangi::element_t<MA>>, Matrix<sangi::element_t<MA>>>
            eigen_generalized_symmetric(const MA& a_basemat, const MB& b_basemat) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_SQUARE(MB);

            Matrix<T> a(a_basemat);  // 1 copy: body modifies a
            Matrix<T> b(b_basemat);  // 1 copy: body modifies b
            const auto n = a.rows();

            if (n != a.cols() || b.rows() != n || b.cols() != n) {
                assert(false && "DimensionError: Symmetric generalized eigenvalue problem requires square matrices of same size");
                throw DimensionError("Symmetric generalized eigenvalue problem requires square matrices of same size");
            }

            // Check symmetry
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = i + 1; j < n; ++j) {
                    if (!approximately_equal(a(i, j), a(j, i)) || !approximately_equal(b(i, j), b(j, i))) {
                        throw MathError("Matrices must be symmetric for eigen_generalized_symmetric");
                    }
                }
            }

            // Arrays storing the eigenvalues and eigenvectors
            Vector<T> eigenvalues(n);
            Matrix<T> eigenvectors(n, n);

            // Solution of the symmetric generalized eigenvalue problem via Cholesky + standard eigenvalue transform
            // (The MKL LAPACKE path is planned for the future since Matrix::data() is not yet supported)

            // Method: use Cholesky decomposition to set L^T * L = B, and
            // convert to the eigenvalue problem of C = L^(-1) * A * L^(-T)

            try {
                // Cholesky decomposition of B: B = L * L^T
                Matrix<T> L = cholesky_decomposition(b);

                // Build C = L^{-1} A L^{-T} via triangular solves (without explicitly forming the inverse)
                // Step 1: Y = L^{-1} A  (solve each column by forward substitution)
                Matrix<T> Y(n, n);
                for (std::size_t col = 0; col < n; ++col) {
                    // Solve L * y_col = A_col
                    for (std::size_t i = 0; i < n; ++i) {
                        T sum = a(i, col);
                        for (std::size_t k = 0; k < i; ++k)
                            sum -= L(i, k) * Y(k, col);
                        Y(i, col) = sum / L(i, i);
                    }
                }

                // Step 2: C = Y * L^{-T}
                // Since C^T = L^{-1} Y^T, for each column j of C^T (= each row j of C),
                // solve L * z_j = column j of Y^T (= row j of Y) by forward substitution
                // → z_j = column j of C^T = row j of C
                Matrix<T> C(n, n);
                for (std::size_t row = 0; row < n; ++row) {
                    // Solve L * z = Y(row, :)^T by forward substitution → z = C(row, :)^T
                    for (std::size_t i = 0; i < n; ++i) {
                        T sum = Y(row, i);
                        for (std::size_t k = 0; k < i; ++k)
                            sum -= L(i, k) * C(row, k);
                        C(row, i) = sum / L(i, i);
                    }
                }

                // Compute the eigenvalues and eigenvectors of C
                auto [c_eigenvalues, c_eigenvectors] = eigen_symmetric(C);

                eigenvalues = c_eigenvalues;

                // Transform the eigenvectors: v = L^{-T} * c_v
                // Solve L^T * v = c_v by backward substitution
                for (std::size_t j = 0; j < n; ++j) {
                    for (std::size_t ii = 0; ii < n; ++ii) {
                        std::size_t i = n - 1 - ii;
                        T sum = c_eigenvectors(i, j);
                        for (std::size_t k = i + 1; k < n; ++k)
                            sum -= L(k, i) * eigenvectors(k, j);
                        eigenvectors(i, j) = sum / L(i, i);
                    }
                }
            }
            catch (const std::exception&) {
                throw MathError("Symmetric generalized eigenvalue calculation failed: B matrix is not positive definite");
            }

            return { eigenvalues, eigenvectors };
        }

        //-----------------------------------------------------------------------------
        // Special eigenvalue algorithms
        //-----------------------------------------------------------------------------

        // Power method (computes the eigenvalue of largest absolute value and its eigenvector)
        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
        std::pair<sangi::element_t<MA>, Vector<sangi::element_t<MA>>> power_method(const MA& a,
            sangi::element_t<MA> tolerance = sangi::element_t<MA>(1e-10),
            std::size_t max_iterations = 1000) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);

            const auto n = a.rows();

            if (n != a.cols()) {
                assert(false && "DimensionError: Power method requires a square matrix");
                throw DimensionError("Power method requires a square matrix");
            }

            // dgebal-equivalent balancing: A → P^T D^{-1} A D P resolves row/column norm imbalance.
            // Eigenvalues are preserved by the similarity transform and the dominant eigenvalue is unchanged.
            // Improves the convergence rate of the power method on ill-conditioned systems (the convergence
            // rate depends on the |λ_2 / λ_1| ratio, but balance suppresses round-off-induced amplitude drift).
            // Well-conditioned systems are a no-op.
            Matrix<T> A_bal(a);
            auto bal_info = balance_matrix(A_bal);

            // Initial vector is random or a unit vector
            Vector<T> x(n, T(1));
            x[0] = T(1);  // set the first element to 1

            // Normalize the initial vector
            const T x_norm = std::sqrt(dot(x, x));
            for (std::size_t i = 0; i < n; ++i) {
                x[i] /= x_norm;
            }

            // Initial estimate of the eigenvalue
            T lambda = T(0);
            T lambda_prev;

            // Power method iteration (performed on the balanced matrix A_bal)
            for (std::size_t iter = 0; iter < max_iterations; ++iter) {
                // Save the previous eigenvalue
                lambda_prev = lambda;

                // y = A_bal * x
                Vector<T> y(n, T(0));
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        y[i] += A_bal(i, j) * x[j];
                    }
                }

                // Rayleigh quotient λ = (x^T * A * x) / (x^T * x)
                lambda = dot(x, y);

                // Normalize the new vector
                const T y_norm = std::sqrt(dot(y, y));
                for (std::size_t i = 0; i < n; ++i) {
                    x[i] = y[i] / y_norm;
                }

                // Convergence test
                if (std::abs(lambda - lambda_prev) < tolerance * std::abs(lambda)) {
                    break;
                }

                // Check the maximum number of iterations
                if (iter == max_iterations - 1) {
                    throw std::runtime_error("Power method did not converge within "
                        + std::to_string(max_iterations) + " iterations");
                }
            }

            // Map the eigenvector of the balanced matrix back to the eigenvector of the original A, and re-normalize
            unbalance_vector(x, bal_info);
            const T xn = std::sqrt(dot(x, x));
            if (xn > T(0)) {
                for (std::size_t i = 0; i < n; ++i) x[i] /= xn;
            }

            return { lambda, x };
        }

        // Inverse power method (computes the eigenvalue closest to a specified value and its eigenvector)
        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
        std::pair<sangi::element_t<MA>, Vector<sangi::element_t<MA>>> inverse_power_method(const MA& a,
            sangi::element_t<MA> sigma = sangi::element_t<MA>(0),  // shift (target value)
            sangi::element_t<MA> tolerance = sangi::element_t<MA>(1e-10),
            std::size_t max_iterations = 1000) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);

            const auto n = a.rows();

            if (n != a.cols()) {
                assert(false && "DimensionError: Inverse power method requires a square matrix");
                throw DimensionError("Inverse power method requires a square matrix");
            }

            // dgebal-equivalent balancing: A → P^T D^{-1} A D P.
            // Because it is a similarity transform, the (A - σI) shift and eigenvalues are preserved.
            // The numerical stability of the LU decomposition improves, and the convergence
            // reliability of inverse iteration on ill-conditioned systems is improved.
            Matrix<T> a_balanced(a);
            auto bal_info = balance_matrix(a_balanced);

            // Computation of (A_balanced - σI)
            Matrix<T> a_shifted = a_balanced;
            for (std::size_t i = 0; i < n; ++i) {
                a_shifted(i, i) -= sigma;
            }

            // LU decomposition of (A - σI)
            auto [lu, pivots] = lu_decomposition(a_shifted);

            // Initial vector: all components non-zero (so that even a diagonal matrix can access all eigenspaces)
            Vector<T> x(n, T(1));

            // Normalize the initial vector
            {
                const T x_norm = std::sqrt(dot(x, x));
                for (std::size_t i = 0; i < n; ++i) {
                    x[i] /= x_norm;
                }
            }

            // Initial estimate of the eigenvalue
            T mu = T(0);
            T mu_prev;

            // Inverse power method iteration
            for (std::size_t iter = 0; iter < max_iterations; ++iter) {
                // Save the previous eigenvalue
                mu_prev = mu;

                // Solve (A - σI)y = x (using the already LU-decomposed matrix)
                Vector<T> y = x;
                // pivot permutation
                for (typename Matrix<T>::size_type pi = 0; pi < n; ++pi) {
                    if (pivots[pi] != pi) {
                        std::swap(y[pi], y[pivots[pi]]);
                    }
                }
                // forward substitution
                for (typename Matrix<T>::size_type fi = 1; fi < n; ++fi) {
                    for (typename Matrix<T>::size_type fj = 0; fj < fi; ++fj) {
                        y[fi] -= lu(fi, fj) * y[fj];
                    }
                }
                // backward substitution
                for (typename Matrix<T>::size_type bi = n; bi-- > 0; ) {
                    for (typename Matrix<T>::size_type bj = bi + 1; bj < n; ++bj) {
                        y[bi] -= lu(bi, bj) * y[bj];
                    }
                    y[bi] /= lu(bi, bi);
                }

                // Rayleigh quotient of (A-σI)^{-1}: μ = x^T * y (x is already normalized)
                mu = dot(x, y);

                // Normalize the new vector
                const T y_norm = std::sqrt(dot(y, y));
                for (std::size_t i = 0; i < n; ++i) {
                    x[i] = y[i] / y_norm;
                }

                // Convergence test
                if (std::abs(mu - mu_prev) < tolerance * std::abs(mu)) {
                    break;
                }

                // Check the maximum number of iterations
                if (iter == max_iterations - 1) {
                    throw std::runtime_error("Inverse power method did not converge within "
                        + std::to_string(max_iterations) + " iterations");
                }
            }

            // Compute the actual eigenvalue (undo the shift)
            const T lambda = sigma + T(1) / mu;

            // Map the eigenvector of the balanced matrix back to the eigenvector of the original A, and re-normalize
            unbalance_vector(x, bal_info);
            const T xn = std::sqrt(dot(x, x));
            if (xn > T(0)) {
                for (std::size_t i = 0; i < n; ++i) x[i] /= xn;
            }

            return { lambda, x };
        }

        // Rayleigh quotient iteration (eigenvalue computation with fast convergence)
        template<typename MA, typename V0>
            requires sangi::BaseMatrixLike<MA>
                  && sangi::BaseVectorLike<V0>
                  && std::same_as<sangi::element_t<MA>, sangi::element_t<V0>>
        std::pair<sangi::element_t<MA>, Vector<sangi::element_t<MA>>> rayleigh_quotient_iteration(const MA& a,
            const V0& initial_guess,
            sangi::element_t<MA> tolerance = sangi::element_t<MA>(1e-10),
            std::size_t max_iterations = 100) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);
            SANGI_STATIC_ASSERT_MATRIX_ROWS_EQ_VECTOR_SIZE(MA, V0);

            const auto n = a.rows();

            if (n != a.cols()) {
                assert(false && "DimensionError: Rayleigh quotient iteration requires a square matrix");
                throw DimensionError("Rayleigh quotient iteration requires a square matrix");
            }

            if (initial_guess.size() != n) {
                assert(false && "DimensionError: Initial guess vector size must match matrix dimension");
                throw DimensionError("Initial guess vector size must match matrix dimension");
            }

            // dgebal-equivalent balancing: A → P^T D^{-1} A D P balances the row/col norms.
            // Eigenvalues are unchanged (similarity transform). The numerical stability of the
            // LU decomposition improves, suppressing the divergence risk of Rayleigh quotient
            // iteration on ill-conditioned systems.
            // The user-specified initial_guess is in the original coordinate system, so convert
            // it to balanced coordinates with balance_vector, and map back to the original
            // coordinates with unbalance_vector at the end.
            Matrix<T> A_bal(a);
            auto bal_info = balance_matrix(A_bal);

            // Copy the initial vector + transform to balanced coordinates + normalize
            Vector<T> x = initial_guess;
            balance_vector(x, bal_info);
            const T x_norm = std::sqrt(dot(x, x));
            if (x_norm == T(0)) {
                throw DimensionError("rayleigh_quotient_iteration: initial guess maps to zero in balanced coordinates");
            }
            for (std::size_t i = 0; i < n; ++i) {
                x[i] /= x_norm;
            }

            // Initial Rayleigh quotient (computed on the balanced matrix; the eigenvalues are identical)
            T rho = T(0);
            for (std::size_t i = 0; i < n; ++i) {
                for (std::size_t j = 0; j < n; ++j) {
                    rho += x[i] * A_bal(i, j) * x[j];
                }
            }

            // Rayleigh quotient iteration (= performed on the balanced matrix A_bal)
            for (std::size_t iter = 0; iter < max_iterations; ++iter) {
                // Computation of (A_bal - ρI)
                Matrix<T> a_shifted = A_bal;
                for (std::size_t i = 0; i < n; ++i) {
                    a_shifted(i, i) -= rho;
                }

                // LU decomposition of (A - ρI)
                Matrix<T> lu;
                std::vector<typename Matrix<T>::size_type> pivots;

                try {
                    auto decomp = lu_decomposition(a_shifted);
                    lu = std::move(decomp.first);
                    pivots = std::move(decomp.second);
                }
                catch (const std::exception&) {
                    // When the LU decomposition fails (singular matrix)
                    // This can occur when ρ is very close to an eigenvalue
                    // Add a tiny perturbation and retry
                    const T epsilon = numeric_traits<T>::epsilon() * T(100);
                    rho += epsilon;
                    continue;
                }

                // Solve (A - ρI)y = x (using the already LU-decomposed matrix)
                Vector<T> y = x;
                for (typename Matrix<T>::size_type pi = 0; pi < n; ++pi) {
                    if (pivots[pi] != pi) {
                        std::swap(y[pi], y[pivots[pi]]);
                    }
                }
                for (typename Matrix<T>::size_type fi = 1; fi < n; ++fi) {
                    for (typename Matrix<T>::size_type fj = 0; fj < fi; ++fj) {
                        y[fi] -= lu(fi, fj) * y[fj];
                    }
                }
                for (typename Matrix<T>::size_type bi = n; bi-- > 0; ) {
                    for (typename Matrix<T>::size_type bj = bi + 1; bj < n; ++bj) {
                        y[bi] -= lu(bi, bj) * y[bj];
                    }
                    y[bi] /= lu(bi, bi);
                }

                // Normalize the new vector
                const T y_norm = std::sqrt(dot(y, y));
                for (std::size_t i = 0; i < n; ++i) {
                    y[i] /= y_norm;
                }

                // Compute the new Rayleigh quotient (on the balanced matrix)
                T rho_new = T(0);
                for (std::size_t i = 0; i < n; ++i) {
                    for (std::size_t j = 0; j < n; ++j) {
                        rho_new += y[i] * A_bal(i, j) * y[j];
                    }
                }

                // Convergence test
                if (std::abs(rho_new - rho) < tolerance * std::abs(rho_new)) {
                    rho = rho_new;
                    x = y;
                    break;
                }

                // Update
                rho = rho_new;
                x = y;

                // Check the maximum number of iterations
                if (iter == max_iterations - 1) {
                    throw std::runtime_error("Rayleigh quotient iteration did not converge within "
                        + std::to_string(max_iterations) + " iterations");
                }
            }

            // Map the eigenvector of the balanced matrix back to the eigenvector of the original A, and re-normalize
            unbalance_vector(x, bal_info);
            const T xn = std::sqrt(dot(x, x));
            if (xn > T(0)) {
                for (std::size_t i = 0; i < n; ++i) x[i] /= xn;
            }

            return { rho, x };
        }

        //-----------------------------------------------------------------------------
        // Gershgorin Circle Theorem
        //-----------------------------------------------------------------------------
        // Each eigenvalue lies in at least one Gershgorin disc D(a_ii, R_i).
        // D(a_ii, R_i) = { z ∈ C : |z - a_ii| ≤ R_i }
        // R_i = Σ_{j≠i} |a_ij|

        // Gershgorin disc: a (center, radius) pair
        template<typename T>
        struct GershgorinDisc {
            T center;   // a_ii (diagonal element)
            T radius;   // Σ_{j≠i} |a_ij|
        };

        // Returns all Gershgorin discs
        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
        std::vector<GershgorinDisc<sangi::element_t<MA>>> gershgorin_discs(const MA& A) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);

            const auto n = A.rows();
            if (n != A.cols()) {
                assert(false && "DimensionError: gershgorin_discs: matrix must be square");
                throw DimensionError("gershgorin_discs: matrix must be square");
            }

            std::vector<GershgorinDisc<T>> discs(n);
            for (std::size_t i = 0; i < n; ++i) {
                discs[i].center = A(i, i);
                T r = T(0);
                for (std::size_t j = 0; j < n; ++j)
                    if (j != i) r += std::abs(A(i, j));
                discs[i].radius = r;
            }
            return discs;
        }

        // Returns the eigenvalue existence range [min, max] (intersection of row-based and column-based)
        template<typename MA>
            requires sangi::BaseMatrixLike<MA>
        std::pair<sangi::element_t<MA>, sangi::element_t<MA>> gershgorin_bounds(const MA& A) {
            using T = sangi::element_t<MA>;
            SANGI_STATIC_ASSERT_SQUARE(MA);

            const auto n = A.rows();
            if (n != A.cols()) {
                assert(false && "DimensionError: gershgorin_bounds: matrix must be square");
                throw DimensionError("gershgorin_bounds: matrix must be square");
            }

            // Row-based discs (initial value set from the first row)
            T r0 = T(0);
            for (std::size_t j = 1; j < n; ++j) r0 += std::abs(A(0, j));
            T lo = A(0, 0) - r0;
            T hi = A(0, 0) + r0;
            for (std::size_t i = 0; i < n; ++i) {
                T r = T(0);
                for (std::size_t j = 0; j < n; ++j)
                    if (j != i) r += std::abs(A(i, j));
                lo = std::min(lo, A(i, i) - r);
                hi = std::max(hi, A(i, i) + r);
            }

            // Column-based discs (Gershgorin discs of A^T)
            for (std::size_t j = 0; j < n; ++j) {
                T r = T(0);
                for (std::size_t i = 0; i < n; ++i)
                    if (i != j) r += std::abs(A(i, j));
                lo = std::max(lo, A(j, j) - r);  // intersection → max of lower bounds
                hi = std::min(hi, A(j, j) + r);  // intersection → min of upper bounds
            }

            // If the intersection is empty, fall back to the row-based result
            if (lo > hi) {
                lo = A(0, 0) - r0;
                hi = A(0, 0) + r0;
                for (std::size_t i = 0; i < n; ++i) {
                    T r = T(0);
                    for (std::size_t j = 0; j < n; ++j)
                        if (j != i) r += std::abs(A(i, j));
                    lo = std::min(lo, A(i, i) - r);
                    hi = std::max(hi, A(i, i) + r);
                }
            }

            return {lo, hi};
        }

        //-----------------------------------------------------------------------------
        // Complex matrix eigenvalue solver (ComplexEigenSolver)
        //-----------------------------------------------------------------------------
        // Computes the eigenvalues and eigenvectors of a general complex matrix A ∈ C^{n×n}.
        // Method: complex Hessenberg decomposition → complex QR iteration (single shift)

        template<typename T>
        std::pair<Vector<Complex<T>>, Matrix<Complex<T>>>
        eigen_complex(const BaseMatrix<Complex<T>>& A) {
            using C = Complex<T>;
            const auto n = A.rows();
            if (n != A.cols()) {
                assert(false && "DimensionError: eigen_complex: matrix must be square");
                throw DimensionError("eigen_complex: matrix must be square");
            }

            if (n == 0) return {{}, {}};

            if (n == 1) {
                Vector<C> evals(1);
                evals[0] = A(0, 0);
                Matrix<C> evecs(1, 1, C(T(1), T(0)));
                return {evals, evecs};
            }

            // Hermitian check: if A = A* the eigenvalues are real (fast path)
            bool is_hermitian = true;
            for (std::size_t i = 0; i < n && is_hermitian; ++i)
                for (std::size_t j = i + 1; j < n && is_hermitian; ++j)
                    if (sangi::abs(A(i, j) - sangi::conj(A(j, i))) >
                        numeric_traits<T>::epsilon() * T(100))
                        is_hermitian = false;

            // Step 1: complex Hessenberg decomposition A = Q H Q*
            // Convert to upper Hessenberg form via Householder transforms
            Matrix<C> H = A;
            Matrix<C> Q(n, n, C(T(0)));
            for (std::size_t i = 0; i < n; ++i) Q(i, i) = C(T(1));

            for (std::size_t k = 0; k + 2 < n; ++k) {
                // Compute the Householder vector of v = H(k+1:n, k)
                std::size_t m = n - k - 1;
                std::vector<C> v(m);
                for (std::size_t i = 0; i < m; ++i)
                    v[i] = H(k + 1 + i, k);

                T nrm = T(0);
                for (auto& vi : v) nrm += sangi::normSq(vi);
                nrm = std::sqrt(nrm);

                if (nrm < numeric_traits<T>::epsilon()) continue;

                // Align the phase to v[0]
                T abs_v0 = sangi::abs(v[0]);
                C alpha = (abs_v0 > T(0)) ?
                    C(-nrm) * v[0] / C(abs_v0) : C(nrm, T(0));
                v[0] -= alpha;

                // Normalize v
                T vnrm = T(0);
                for (auto& vi : v) vnrm += sangi::normSq(vi);
                vnrm = std::sqrt(vnrm);
                if (vnrm < numeric_traits<T>::epsilon()) continue;
                for (auto& vi : v) vi /= vnrm;

                // H ← (I - 2vv*) H (I - 2vv*)
                // from the left: H(k+1:n, :) -= 2 v (v* H(k+1:n, :))
                for (std::size_t j = k; j < n; ++j) {
                    C dot(T(0));
                    for (std::size_t i = 0; i < m; ++i)
                        dot += sangi::conj(v[i]) * H(k + 1 + i, j);
                    dot *= C(T(2));
                    for (std::size_t i = 0; i < m; ++i)
                        H(k + 1 + i, j) -= v[i] * dot;
                }
                // from the right: H(:, k+1:n) -= 2 (H(:, k+1:n) v) v*
                for (std::size_t i = 0; i < n; ++i) {
                    C dot(T(0));
                    for (std::size_t j = 0; j < m; ++j)
                        dot += H(i, k + 1 + j) * v[j];
                    dot *= C(T(2));
                    for (std::size_t j = 0; j < m; ++j)
                        H(i, k + 1 + j) -= dot * sangi::conj(v[j]);
                }
                // Q accumulation: Q(:, k+1:n) -= 2 (Q(:, k+1:n) v) v*
                for (std::size_t i = 0; i < n; ++i) {
                    C dot(T(0));
                    for (std::size_t j = 0; j < m; ++j)
                        dot += Q(i, k + 1 + j) * v[j];
                    dot *= C(T(2));
                    for (std::size_t j = 0; j < m; ++j)
                        Q(i, k + 1 + j) -= dot * sangi::conj(v[j]);
                }
            }

            // Clean up the subdiagonal elements
            for (std::size_t i = 0; i + 2 < n; ++i)
                for (std::size_t j = 0; j < i; ++j)
                    H(i + 2, j) = C(T(0));

            // Step 2: complex QR iteration (single shift + Givens rotations)
            const int max_iter = 300 * static_cast<int>(n);
            const T eps = numeric_traits<T>::epsilon();
            int ihi = static_cast<int>(n);
            int total_iter = 0;

            while (ihi > 1 && total_iter < max_iter) {
                // deflation
                int ilo = ihi - 1;
                while (ilo > 0) {
                    T threshold = eps * (sangi::abs(H(ilo - 1, ilo - 1)) +
                                         sangi::abs(H(ilo, ilo)));
                    if (threshold == T(0)) threshold = eps;
                    if (sangi::abs(H(ilo, ilo - 1)) <= threshold) {
                        H(ilo, ilo - 1) = C(T(0));
                        break;
                    }
                    --ilo;
                }

                if (ilo == ihi - 1) {
                    --ihi;
                    continue;
                }

                // Wilkinson shift: eigenvalue close to H(ihi-1, ihi-1)
                C sigma = H(ihi - 1, ihi - 1);

                // QR step: H - σI = QR via Givens rotations → H ← RQ + σI
                for (int i = ilo; i < ihi - 1; ++i) {
                    C x = H(i, i) - sigma;
                    C y = H(i + 1, i);
                    T r = std::sqrt(sangi::normSq(x) + sangi::normSq(y));
                    if (r < eps) { ++total_iter; continue; }
                    C c = x / C(r);
                    C s = y / C(r);

                    // rotate from the left: update H([i, i+1], :)
                    for (int j = ilo; j < static_cast<int>(n); ++j) {
                        C t1 = H(i, j);
                        C t2 = H(i + 1, j);
                        H(i, j)     = sangi::conj(c) * t1 + sangi::conj(s) * t2;
                        H(i + 1, j) = -s * t1 + c * t2;
                    }
                    // rotate from the right: update H(:, [i, i+1])
                    int limit = std::min(i + 3, ihi);
                    for (int j = 0; j < limit; ++j) {
                        C t1 = H(j, i);
                        C t2 = H(j, i + 1);
                        H(j, i)     = c * t1 + s * t2;
                        H(j, i + 1) = -sangi::conj(s) * t1 + sangi::conj(c) * t2;
                    }
                    // Q accumulation
                    for (std::size_t j = 0; j < n; ++j) {
                        C t1 = Q(j, i);
                        C t2 = Q(j, i + 1);
                        Q(j, i)     = c * t1 + s * t2;
                        Q(j, i + 1) = -sangi::conj(s) * t1 + sangi::conj(c) * t2;
                    }
                }
                ++total_iter;
            }

            // Step 3: extract the eigenvalues from the diagonal
            Vector<C> eigenvalues(n);
            for (std::size_t i = 0; i < n; ++i)
                eigenvalues[i] = H(i, i);

            // Step 4: eigenvectors (backward substitution + Q transform)
            Matrix<C> eigenvectors(n, n, C(T(0)));
            for (std::size_t k = 0; k < n; ++k) {
                C lambda = eigenvalues[k];
                std::vector<C> v(n, C(T(0)));
                v[k] = C(T(1));

                for (int i = static_cast<int>(k) - 1; i >= 0; --i) {
                    C sum(T(0));
                    for (std::size_t j = i + 1; j <= k; ++j)
                        sum += H(i, j) * v[j];
                    C diag = H(i, i) - lambda;
                    if (sangi::abs(diag) > eps * T(100))
                        v[i] = -sum / diag;
                    else
                        v[i] = -sum / C(eps * T(100));
                }

                // Q * v
                for (std::size_t i = 0; i < n; ++i) {
                    C sum(T(0));
                    for (std::size_t j = 0; j < n; ++j)
                        sum += Q(i, j) * v[j];
                    eigenvectors(i, k) = sum;
                }

                // normalization
                T norm_sq = T(0);
                for (std::size_t i = 0; i < n; ++i)
                    norm_sq += sangi::normSq(eigenvectors(i, k));
                T nrm = std::sqrt(norm_sq);
                if (nrm > T(0))
                    for (std::size_t i = 0; i < n; ++i)
                        eigenvectors(i, k) = eigenvectors(i, k) / nrm;
            }

            return {eigenvalues, eigenvectors};
        }

    } // namespace algorithms
} // namespace sangi

#endif // SANGI_EIGENVALUES_HPP