Add some QR algos

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<std::size_t> 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<std::size_t> perm(n);
+    std::iota(perm.begin(), perm.end(), std::size_t{0});
+
+    // Precompute column norms squared.
+    std::vector<double> 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<double> 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<std::size_t>(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<double> 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));
+}