From 92220ea5a483d6ece73bb6472af0773bc4106d73 Mon Sep 17 00:00:00 2001 From: y-jan137 Date: Sun, 3 May 2026 19:08:14 +0300 Subject: Add some QR algos --- README.md | 10 +- src/qr.cpp | 134 +++++++++++++++++++++++ src/qr_iteration.cpp | 258 ++++++++++++++++++++++++++++++++++++++++++++ tests/test_qr.cpp | 61 +++++++++++ tests/test_qr_iteration.cpp | 88 +++++++++++++++ 5 files changed, 547 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 0560b53..00105cd 100644 --- a/README.md +++ b/README.md @@ -57,18 +57,20 @@ ctest --test-dir build --output-on-failure - LU factorization with partial pivoting (`lu_factor`, `lu_solve`) - QR factorization — classical GS, modified GS, and Householder (`qr_classical_gs`, `qr_modified_gs`, `qr_householder`) +- Rank-revealing QR — Householder with column pivoting (`qr_colpiv`); reports numerical rank + and ensures |R(i,i)| ≥ |R(i+1,i+1)| - Eigenvalue computation via QR iteration: - Unshifted QR (`eigenvalues_unshifted`) — linear convergence, T&B Algorithm 28.1 - Wilkinson-shifted QR (`eigenvalues_shifted`) — typically cubic convergence, T&B Lecture 29 - Hessenberg + Givens QR (`eigenvalues_hessenberg`) — O(n²) per step after one O(n³) reduction; ~10–30× faster than `eigenvalues_shifted` for n ≥ 50 + - Francis double-shift QR (`eigenvalues_francis`) — implicit bulge chasing on Hessenberg form; + handles complex conjugate eigenvalue pairs without complex arithmetic; robust subdiagonal + + 2×2 block deflation with exceptional shifts (GVL §7.5) +- Cholesky factorization (`cholesky_factor`, `cholesky_solve`) — for symmetric positive definite systems ## TODO: -- [ ] Cholesky factorization (cholesky) — for symmetric positive definite systems -- [ ] Rank-revealing QR — Householder QR with column pivoting (qr_colpiv) - [ ] Symmetric tridiagonalization — Householder reduction before symmetric QR (tridiagonalize) -- [ ] Francis double-shift QR — bulge chasing for real matrices with complex conjugate eigenvalue pairs (eigenvalues_francis) -- [ ] Deflation — robust subdiagonal + 2×2 block deflation in Hessenberg QR - [ ] Eigenvectors via inverse iteration (eigenvectors_inverse_iteration) - [ ] SVD — Golub-Kahan bidiagonalization + QR (svd) - [ ] Conjugate Gradient (solve_cg) — for symmetric positive definite systems diff --git a/src/qr.cpp b/src/qr.cpp index 7e32ce3..a907bab 100644 --- a/src/qr.cpp +++ b/src/qr.cpp @@ -37,6 +37,24 @@ QRResult qr_modified_gs(const Matrix& A, double zero_tolerance = 1e-14); // Throws DimensionMismatchError if rows < cols. QRResult qr_householder(const Matrix& A); +struct QRColPivResult { + Matrix Q; + Matrix R; + std::vector perm; + std::size_t rank; +}; + +// Householder QR with column pivoting (rank-revealing). +// At each step, the column with largest remaining norm is selected as pivot. +// The numerical rank is determined by comparing diagonal entries of R to +// rank_tolerance * |R(0,0)|. +// +// Returns Q (m x n), R (n x n upper triangular), perm (column permutation), +// and rank (numerical rank estimate). +// +// Throws DimensionMismatchError if rows < cols. +QRColPivResult qr_colpiv(const Matrix& A, double rank_tolerance = 1e-12); + } // namespace linalgebra namespace { @@ -194,4 +212,120 @@ QRResult qr_householder(const Matrix& A) { return QRResult{std::move(Q), std::move(R)}; } +QRColPivResult qr_colpiv(const Matrix& A, double rank_tolerance) { + require_tall(A, "qr_colpiv"); + const std::size_t m = A.rows(); + const std::size_t n = A.cols(); + + Matrix work = A; + Matrix Q_full = Matrix::identity(m); + + std::vector perm(n); + std::iota(perm.begin(), perm.end(), std::size_t{0}); + + // Precompute column norms squared. + std::vector col_norms_sq(n); + for (std::size_t j = 0; j < n; ++j) { + double s = 0.0; + for (std::size_t i = 0; i < m; ++i) s += work(i, j) * work(i, j); + col_norms_sq[j] = s; + } + + std::size_t rank = n; + + for (std::size_t k = 0; k < n; ++k) { + // Find pivot: column with largest remaining norm. + std::size_t pivot = k; + double max_norm = col_norms_sq[k]; + for (std::size_t j = k + 1; j < n; ++j) { + if (col_norms_sq[j] > max_norm) { + max_norm = col_norms_sq[j]; + pivot = j; + } + } + + // Swap columns k and pivot. + if (pivot != k) { + for (std::size_t i = 0; i < m; ++i) { + std::swap(work(i, k), work(i, pivot)); + } + std::swap(col_norms_sq[k], col_norms_sq[pivot]); + std::swap(perm[k], perm[pivot]); + } + + // Householder reflector for column k. + const std::size_t p = m - k; + + std::vector u(p); + for (std::size_t i = 0; i < p; ++i) u[i] = work(k + i, k); + + const double x_norm = [&] { + double s = 0.0; + for (double v : u) s += v * v; + return std::sqrt(s); + }(); + + if (x_norm == 0.0) { + // Remaining columns are zero — rank determined. + rank = k; + break; + } + + // Check rank: if this pivot norm is small relative to R(0,0). + if (k > 0) { + const double r00 = std::abs(work(0, 0)); + if (x_norm <= rank_tolerance * r00) { + rank = k; + break; + } + } + + const double sigma = (u[0] >= 0.0 ? 1.0 : -1.0) * x_norm; + u[0] += sigma; + + const double utu = [&] { + double s = 0.0; + for (double v : u) s += v * v; + return s; + }(); + const double tau = 2.0 / utu; + + // Apply reflector to work columns k..n-1. + for (std::size_t j = k; j < n; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += u[i] * work(k + i, j); + const double coeff = tau * d; + for (std::size_t i = 0; i < p; ++i) work(k + i, j) -= coeff * u[i]; + } + + // Apply reflector to Q_full. + for (std::size_t j = 0; j < m; ++j) { + double d = 0.0; + for (std::size_t i = 0; i < p; ++i) d += u[i] * Q_full(k + i, j); + const double coeff = tau * d; + for (std::size_t i = 0; i < p; ++i) Q_full(k + i, j) -= coeff * u[i]; + } + + // Update column norms (downdate). + for (std::size_t j = k + 1; j < n; ++j) { + const double val = work(k, j); + col_norms_sq[j] -= val * val; + if (col_norms_sq[j] < 0.0) col_norms_sq[j] = 0.0; + } + } + + // Extract Q (m x n) and R (n x n). + Matrix Q(m, n); + for (std::size_t i = 0; i < m; ++i) + for (std::size_t j = 0; j < n; ++j) + Q(i, j) = Q_full(j, i); + + Matrix R(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + R(i, j) = work(i, j); + + return QRColPivResult{std::move(Q), std::move(R), std::move(perm), rank}; +} + } // namespace linalgebra diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp index 64fd25f..f8a3be5 100644 --- a/src/qr_iteration.cpp +++ b/src/qr_iteration.cpp @@ -65,6 +65,13 @@ void hessenberg_qr_step(Matrix& H, double sigma); [[nodiscard]] QRIterationResult eigenvalues_hessenberg(const Matrix& A, QRIterationOptions opts = {}); +// Francis double-shift QR — implicit bulge chasing on upper Hessenberg form. +// Handles real matrices with complex conjugate eigenvalue pairs without +// complex arithmetic. Uses robust deflation (subdiagonal + 2×2 block). +// Reference: GVL §7.5, T&B Lecture 29. +[[nodiscard]] QRIterationResult eigenvalues_francis(const Matrix& A, + QRIterationOptions opts = {}); + } // namespace linalgebra namespace { @@ -502,4 +509,255 @@ QRIterationResult eigenvalues_hessenberg(const Matrix& A, QRIterationOptions opt return result; } +// --------------------------------------------------------------------------- +// Francis double-shift QR with implicit bulge chasing +// --------------------------------------------------------------------------- + +QRIterationResult eigenvalues_francis(const Matrix& A, QRIterationOptions opts) { + require_square(A, "eigenvalues_francis"); + const std::size_t n = A.rows(); + + QRIterationResult result; + result.eigenvalues_real = Vector(n, 0.0); + result.eigenvalues_imag = Vector(n, 0.0); + + if (opts.track_convergence) + result.convergence_history.reserve( + static_cast(opts.max_iterations)); + + if (n == 1) { + result.eigenvalues_real[0] = A(0, 0); + return result; + } + + if (n == 2) { + const double a = A(0, 0), b = A(0, 1), c = A(1, 0), d = A(1, 1); + const double tr = a + d; + const double disc = (a - d) * (a - d) + 4.0 * b * c; + if (disc >= 0.0) { + const double sq = std::sqrt(disc); + result.eigenvalues_real[0] = 0.5 * (tr + sq); + result.eigenvalues_real[1] = 0.5 * (tr - sq); + } else { + result.eigenvalues_real[0] = 0.5 * tr; + result.eigenvalues_imag[0] = 0.5 * std::sqrt(-disc); + result.eigenvalues_real[1] = 0.5 * tr; + result.eigenvalues_imag[1] = -0.5 * std::sqrt(-disc); + } + return result; + } + + HessenbergResult hr = hessenberg_reduction(A); + Matrix& H = hr.H; + + std::size_t n_found = n; + std::size_t active = n; + + auto store_real = [&](double re) { + --n_found; + result.eigenvalues_real[n_found] = re; + result.eigenvalues_imag[n_found] = 0.0; + }; + + auto store_pair = [&](double re, double im) { + --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = im; + --n_found; result.eigenvalues_real[n_found] = re; result.eigenvalues_imag[n_found] = -im; + }; + + // Robust deflation: checks both subdiagonal magnitude and 2x2 block. + auto deflation_tol = [&](std::size_t i) -> double { + const double scale = std::abs(H(i - 1, i - 1)) + std::abs(H(i, i)); + return opts.tolerance * (scale > 0.0 ? scale : 1.0); + }; + + auto close_2x2 = [&]() { + const double a = H(active - 2, active - 2); + const double b = H(active - 2, active - 1); + const double c = H(active - 1, active - 2); + const double d = H(active - 1, active - 1); + const double tr = a + d; + const double disc = (a - d) * (a - d) + 4.0 * b * c; + if (disc >= 0.0) { + const double sq = std::sqrt(disc); + store_real(0.5 * (tr + sq)); + store_real(0.5 * (tr - sq)); + } else { + store_pair(0.5 * tr, 0.5 * std::sqrt(-disc)); + } + active -= 2; + }; + + // Find the start of the active unreduced block (split from the top). + auto find_block_start = [&]() -> std::size_t { + for (std::size_t i = active - 1; i >= 1; --i) { + if (std::abs(H(i, i - 1)) < deflation_tol(i)) { + H(i, i - 1) = 0.0; + return i; + } + } + return 0; + }; + + int exceptional_shift_count = 0; + + for (int k = 0; k < opts.max_iterations; ++k) { + // Deflate converged eigenvalues from bottom. + while (active >= 2) { + if (std::abs(H(active - 1, active - 2)) < deflation_tol(active - 1)) { + H(active - 1, active - 2) = 0.0; + store_real(H(active - 1, active - 1)); + --active; + } else { + break; + } + } + + if (active == 0) break; + if (active == 1) { store_real(H(0, 0)); active = 0; break; } + if (active == 2) { close_2x2(); break; } + + // Check for 2x2 block deflation (complex pair at bottom). + if (active >= 3 && std::abs(H(active - 2, active - 3)) < deflation_tol(active - 2)) { + H(active - 2, active - 3) = 0.0; + // The bottom 2x2 has converged. + const double a = H(active - 2, active - 2); + const double b = H(active - 2, active - 1); + const double c = H(active - 1, active - 2); + const double d = H(active - 1, active - 1); + const double tr = a + d; + const double disc = (a - d) * (a - d) + 4.0 * b * c; + if (disc >= 0.0) { + const double sq = std::sqrt(disc); + store_real(0.5 * (tr + sq)); + store_real(0.5 * (tr - sq)); + } else { + store_pair(0.5 * tr, 0.5 * std::sqrt(-disc)); + } + active -= 2; + exceptional_shift_count = 0; + continue; + } + + std::size_t block_start = find_block_start(); + + // Compute Francis double shift from bottom 2x2 of active block. + const double a11 = H(active - 2, active - 2); + const double a12 = H(active - 2, active - 1); + const double a21 = H(active - 1, active - 2); + const double a22 = H(active - 1, active - 1); + double s = a11 + a22; // trace of bottom 2x2 + double t = a11 * a22 - a12 * a21; // determinant of bottom 2x2 + + // Exceptional shift (Wilkinson's ad hoc) every 10 iterations to break stalls. + if (exceptional_shift_count > 0 && exceptional_shift_count % 10 == 0) { + const double w = std::abs(H(active - 1, active - 2)) + + std::abs(H(block_start + 1, block_start)); + s = 1.5 * w; + t = w * w; + } + + // First column of M = H^2 - sH + tI (implicit). + const double h00 = H(block_start, block_start); + const double h01 = H(block_start, block_start + 1); + const double h10 = H(block_start + 1, block_start); + const double h11 = H(block_start + 1, block_start + 1); + const double h21 = (block_start + 2 < active) ? H(block_start + 2, block_start + 1) : 0.0; + + double x = h00 * h00 + h01 * h10 - s * h00 + t; + double y = h10 * (h00 + h11 - s); + double z = h10 * h21; + + // Chase the bulge through the Hessenberg matrix. + for (std::size_t i = block_start; i + 2 < active; ++i) { + // Determine Householder reflector P such that P * [x; y; z]^T = [*; 0; 0]^T. + const std::size_t p = (i + 3 <= active) ? 3 : 2; + + double norm_v = std::sqrt(x * x + y * y + (p == 3 ? z * z : 0.0)); + if (norm_v == 0.0) break; + + const double sign = (x >= 0.0) ? 1.0 : -1.0; + double v0 = x + sign * norm_v; + double v1 = y; + double v2 = (p == 3) ? z : 0.0; + + const double vdot = v0 * v0 + v1 * v1 + v2 * v2; + const double tau = 2.0 / vdot; + + // Apply P from left to H rows [i, i+p-1], columns [max(i-1,0), active-1]. + const std::size_t col_start = (i > 0) ? i - 1 : 0; + for (std::size_t j = col_start; j < active; ++j) { + double d = v0 * H(i, j) + v1 * H(i + 1, j); + if (p == 3) d += v2 * H(i + 2, j); + const double coeff = tau * d; + H(i, j) -= coeff * v0; + H(i + 1, j) -= coeff * v1; + if (p == 3) H(i + 2, j) -= coeff * v2; + } + + // Apply P from right to H rows [0, min(i+p, active-1)], columns [i, i+p-1]. + const std::size_t row_end = std::min(i + p + 1, active); + for (std::size_t j = 0; j < row_end; ++j) { + double d = v0 * H(j, i) + v1 * H(j, i + 1); + if (p == 3) d += v2 * H(j, i + 2); + const double coeff = tau * d; + H(j, i) -= coeff * v0; + H(j, i + 1) -= coeff * v1; + if (p == 3) H(j, i + 2) -= coeff * v2; + } + + // Prepare for next bulge step. + if (i + 3 < active) { + x = H(i + 1, i); + y = H(i + 2, i); + z = (i + 3 < active) ? H(i + 3, i) : 0.0; + } + } + + // Final 2x2 reflector to restore Hessenberg form at bottom. + { + const std::size_t i = active - 2; + const double xi = H(i, i - 1); + const double yi = H(i + 1, i - 1); + const double r = std::hypot(xi, yi); + if (r > 0.0) { + const double c = xi / r; + const double s_val = yi / r; + // Apply Givens from left. + for (std::size_t j = i - 1; j < active; ++j) { + const double t0 = H(i, j); + const double t1 = H(i + 1, j); + H(i, j) = c * t0 + s_val * t1; + H(i + 1, j) = -s_val * t0 + c * t1; + } + // Apply Givens from right. + for (std::size_t j = 0; j < std::min(i + 3, active); ++j) { + const double t0 = H(j, i); + const double t1 = H(j, i + 1); + H(j, i) = c * t0 + s_val * t1; + H(j, i + 1) = -s_val * t0 + c * t1; + } + } + } + + ++exceptional_shift_count; + + if (opts.track_convergence) { + double s_norm = 0.0; + for (std::size_t ii = 1; ii < active; ++ii) + s_norm += H(ii, ii - 1) * H(ii, ii - 1); + result.convergence_history.push_back(std::sqrt(s_norm)); + } + ++result.iterations; + } + + if (n_found > 0) { + std::ostringstream oss; + oss << "eigenvalues_francis: did not converge in " + << opts.max_iterations << " iterations (" + << n_found << " eigenvalue(s) not yet deflated)."; + throw NonConvergenceError(oss.str()); + } + return result; +} + } // namespace linalgebra diff --git a/tests/test_qr.cpp b/tests/test_qr.cpp index 1a4114f..d78fbb1 100644 --- a/tests/test_qr.cpp +++ b/tests/test_qr.cpp @@ -223,3 +223,64 @@ TEST_CASE("QR: linearly dependent columns throw from GS methods", "[qr]") { CHECK_THROWS_AS(linalgebra::qr_modified_gs(A), SingularMatrixError); CHECK_NOTHROW(linalgebra::qr_householder(A)); } + +// --------------------------------------------------------------------------- +// qr_colpiv tests +// --------------------------------------------------------------------------- + +TEST_CASE("QR ColPiv: full rank reconstruction", "[qr][colpiv]") { + const Matrix A = random_matrix(6, 4, 100u); + auto result = linalgebra::qr_colpiv(A); + + // A * P = Q * R => reconstruct A(:, perm) = Q * R + const std::size_t m = A.rows(); + const std::size_t n = A.cols(); + + // Build A*P (permuted columns of A). + Matrix AP(m, n); + for (std::size_t j = 0; j < n; ++j) + for (std::size_t i = 0; i < m; ++i) + AP(i, j) = A(i, result.perm[j]); + + const Matrix QR = result.Q * result.R; + double err = 0.0; + for (std::size_t i = 0; i < m; ++i) + for (std::size_t j = 0; j < n; ++j) { + double d = AP(i, j) - QR(i, j); + err += d * d; + } + CHECK(std::sqrt(err) < 1e-12); + CHECK(result.rank == n); +} + +TEST_CASE("QR ColPiv: rank deficient matrix", "[qr][colpiv]") { + // Rank 2 matrix (col 2 = col 0 + col 1). + Matrix A{ + {1.0, 2.0, 3.0}, + {4.0, 5.0, 9.0}, + {7.0, 8.0, 15.0}, + {2.0, 1.0, 3.0} + }; + auto result = linalgebra::qr_colpiv(A); + CHECK(result.rank == 2); +} + +TEST_CASE("QR ColPiv: R diagonal magnitudes are non-increasing", "[qr][colpiv]") { + const Matrix A = random_matrix(8, 5, 42u); + auto result = linalgebra::qr_colpiv(A); + + for (std::size_t i = 0; i + 1 < result.rank; ++i) { + CHECK(std::abs(result.R(i, i)) >= std::abs(result.R(i + 1, i + 1)) - 1e-14); + } +} + +TEST_CASE("QR ColPiv: identity matrix", "[qr][colpiv]") { + auto I = Matrix::identity(4); + auto result = linalgebra::qr_colpiv(I); + CHECK(result.rank == 4); +} + +TEST_CASE("QR ColPiv: fat matrix throws", "[qr][colpiv]") { + Matrix A(3, 5); + CHECK_THROWS_AS(linalgebra::qr_colpiv(A), DimensionMismatchError); +} diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp index 008f109..64a5dd6 100644 --- a/tests/test_qr_iteration.cpp +++ b/tests/test_qr_iteration.cpp @@ -422,3 +422,91 @@ TEST_CASE("Hessenberg QR: faster than naive shifted QR for large n", } std::cout << "=============================================\n"; } + +// --------------------------------------------------------------------------- +// Francis double-shift QR tests +// --------------------------------------------------------------------------- + +TEST_CASE("Francis QR: 2x2 real eigenvalues", "[qr_iteration][francis]") { + Matrix A{{3.0, 1.0}, {0.0, 2.0}}; + const QRIterationResult res = linalgebra::eigenvalues_francis(A); + const EigPairs expected = {{3.0, 0.0}, {2.0, 0.0}}; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-10)); +} + +TEST_CASE("Francis QR: 2x2 complex eigenvalues", "[qr_iteration][francis]") { + // Rotation matrix — eigenvalues are cos(theta) ± i*sin(theta). + const double theta = 1.0; + Matrix A{{std::cos(theta), -std::sin(theta)}, + {std::sin(theta), std::cos(theta)}}; + const QRIterationResult res = linalgebra::eigenvalues_francis(A); + const EigPairs expected = { + {std::cos(theta), std::sin(theta)}, + {std::cos(theta), -std::sin(theta)} + }; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-10)); +} + +TEST_CASE("Francis QR: 3x3 with complex pair", "[qr_iteration][francis]") { + // Block diagonal: 2x2 rotation (complex pair) + real eigenvalue. + Matrix A{{0.0, -1.0, 0.0}, + {1.0, 0.0, 0.0}, + {0.0, 0.0, 5.0}}; + const QRIterationResult res = linalgebra::eigenvalues_francis(A); + const EigPairs expected = {{0.0, 1.0}, {0.0, -1.0}, {5.0, 0.0}}; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-10)); +} + +TEST_CASE("Francis QR: symmetric matrix (all real)", "[qr_iteration][francis]") { + const Matrix A = random_symmetric(10, 77u); + const QRIterationResult ref = linalgebra::eigenvalues_hessenberg(A); + const QRIterationResult res = linalgebra::eigenvalues_francis(A); + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, + to_pairs(ref.eigenvalues_real, ref.eigenvalues_imag), 1e-8)); +} + +TEST_CASE("Francis QR: non-symmetric with complex pairs", "[qr_iteration][francis]") { + // Random non-symmetric matrix — older single-shift methods struggle here, + // but Francis double-shift handles it natively. + std::mt19937 rng(99u); + std::uniform_real_distribution dist(-2.0, 2.0); + const std::size_t n = 8; + Matrix A(n, n); + for (std::size_t i = 0; i < n; ++i) + for (std::size_t j = 0; j < n; ++j) + A(i, j) = dist(rng); + + const QRIterationResult res = linalgebra::eigenvalues_francis(A); + + // Verify: for each eigenvalue λ, check that the characteristic polynomial + // product of (λ_i - λ_j) is consistent — i.e., sum of eigenvalues = trace. + double trace_A = 0.0; + for (std::size_t i = 0; i < n; ++i) trace_A += A(i, i); + + double trace_eigs = 0.0; + for (std::size_t i = 0; i < n; ++i) trace_eigs += res.eigenvalues_real[i]; + CHECK(trace_eigs == Catch::Approx(trace_A).margin(1e-6)); + + // All imaginary parts should come in conjugate pairs. + double imag_sum = 0.0; + for (std::size_t i = 0; i < n; ++i) imag_sum += res.eigenvalues_imag[i]; + CHECK(imag_sum == Catch::Approx(0.0).margin(1e-8)); +} + +TEST_CASE("Francis QR: companion matrix", "[qr_iteration][francis]") { + // Companion matrix for x^4 - 10x^3 + 35x^2 - 50x + 24 = (x-1)(x-2)(x-3)(x-4). + Matrix A{{0.0, 0.0, 0.0, -24.0}, + {1.0, 0.0, 0.0, 50.0}, + {0.0, 1.0, 0.0, -35.0}, + {0.0, 0.0, 1.0, 10.0}}; + const QRIterationResult res = linalgebra::eigenvalues_francis(A); + const EigPairs expected = {{1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0}, {4.0, 0.0}}; + CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8)); +} + +TEST_CASE("Francis QR: 1x1 matrix", "[qr_iteration][francis]") { + Matrix A{{7.0}}; + const QRIterationResult res = linalgebra::eigenvalues_francis(A); + CHECK(res.eigenvalues_real[0] == Catch::Approx(7.0)); + CHECK(res.eigenvalues_imag[0] == Catch::Approx(0.0).margin(1e-15)); +} -- cgit v1.2.3