Add QR iteration and cleanup

10 files changed, 742 insertions, 110 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1d2374c..1f40b0c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,7 +8,6 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 option(LINEAR_ALGEBRA_BUILD_TESTS "Build unit tests" ON)
-option(LINEAR_ALGEBRA_BUILD_EXAMPLES "Build example programs" ON)
 set(LINEAR_ALGEBRA_SIMD "AUTO" CACHE STRING "SIMD backend for matmul: AUTO, NONE, AVX, AVX2, AVX512")
 set_property(CACHE LINEAR_ALGEBRA_SIMD PROPERTY STRINGS AUTO NONE AVX AVX2 AVX512)
 
@@ -19,6 +18,7 @@ add_library(linear_algebra
     src/triangular_solve.cpp
     src/lu.cpp
     src/qr.cpp
+    src/qr_iteration.cpp
 )
 
 add_library(linear_algebra::core ALIAS linear_algebra)
@@ -54,14 +54,6 @@ elseif(MSVC)
     endif()
 endif()
 
-if(LINEAR_ALGEBRA_BUILD_EXAMPLES)
-    add_executable(linear_system examples/linear_system.cpp)
-    target_link_libraries(linear_system PRIVATE linear_algebra::core)
-
-    add_executable(matmul examples/matmul.cpp)
-    target_link_libraries(matmul PRIVATE linear_algebra::core)
-endif()
-
 option(LINEAR_ALGEBRA_BUILD_EXPERIMENTS "Build experiment programs" ON)
 if(LINEAR_ALGEBRA_BUILD_EXPERIMENTS)
     add_executable(pivoting_vs_no_pivoting experiments/pivoting_vs_no_pivoting.cpp)
@@ -69,6 +61,9 @@ if(LINEAR_ALGEBRA_BUILD_EXPERIMENTS)
 
     add_executable(hilbert_qr experiments/hilbert_qr.cpp)
     target_link_libraries(hilbert_qr PRIVATE linear_algebra::core)
+
+    add_executable(matmul experiments/matmul.cpp)
+    target_link_libraries(matmul PRIVATE linear_algebra::core)
 endif()
 
 if(LINEAR_ALGEBRA_BUILD_TESTS)
@@ -93,6 +88,7 @@ if(LINEAR_ALGEBRA_BUILD_TESTS)
         tests/test_triangular_solve.cpp
         tests/test_lu.cpp
         tests/test_qr.cpp
+        tests/test_qr_iteration.cpp
     )
 
     target_link_libraries(linear_algebra_tests
diff --git a/examples/linear_system.cpp b/examples/linear_system.cpp
deleted file mode 100644
index 184f576..0000000
--- a/examples/linear_system.cpp
+++ /dev/null
@@ -1,36 +0,0 @@
-#include "matrix.hpp"
-#include "norms.hpp"
-#include "triangular_solve.hpp"
-#include "vector.hpp"
-
-#include <iomanip>
-#include <iostream>
-
-int main() {
-    const linalg::Matrix basis = linalg::Matrix::identity(3);
-    const linalg::Matrix upper{
-        {4.0, -2.0, 1.0},
-        {0.0, 3.0, 5.0},
-        {0.0, 0.0, -2.0}
-    };
-    const linalg::Vector expected{2.0, -1.0, 3.0};
-    const linalg::Vector rhs = upper * expected;
-    const linalg::Vector x = linalg::backward_substitution(upper, rhs);
-
-    std::cout << "Week 3 triangular solve demo\n";
-    std::cout << "Identity matrix diagonal: ";
-    for (std::size_t i = 0; i < basis.rows(); ++i) {
-        std::cout << basis(i, i) << (i + 1 == basis.rows() ? '\n' : ' ');
-    }
-
-    std::cout << "Recovered solution x: ";
-    for (std::size_t i = 0; i < x.size(); ++i) {
-        std::cout << std::fixed << std::setprecision(2) << x[i]
-                  << (i + 1 == x.size() ? '\n' : ' ');
-    }
-
-    const linalg::Vector residual = (upper * x) - rhs;
-    std::cout << "Residual 2-norm = " << linalg::norm2(residual) << '\n';
-
-    return 0;
-}
diff --git a/examples/matmul.cpp b/examples/matmul.cpp
deleted file mode 100644
index c6245f3..0000000
--- a/examples/matmul.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#include "matrix.hpp"
-
-#include <iomanip>
-#include <iostream>
-
-int main() {
-    const linalg::Matrix a{
-        {1.0, 2.0, 3.0},
-        {4.0, 5.0, 6.0}
-    };
-    const linalg::Matrix b{
-        {7.0, 8.0},
-        {9.0, 10.0},
-        {11.0, 12.0}
-    };
-
-    const linalg::Matrix c = a * b;
-
-    std::cout << "Matrix multiplication\n";
-    std::cout << "A is " << a.rows() << " x " << a.cols() << '\n';
-    std::cout << "B is " << b.rows() << " x " << b.cols() << '\n';
-    std::cout << "C = A * B:\n";
-
-    for (std::size_t i = 0; i < c.rows(); ++i) {
-        for (std::size_t j = 0; j < c.cols(); ++j) {
-            std::cout << std::fixed << std::setprecision(2) << c(i, j)
-                      << (j + 1 == c.cols() ? '\n' : ' ');
-        }
-    }
-
-    return 0;
-}
diff --git a/experiments/hilbert_qr.cpp b/experiments/hilbert_qr.cpp
index 9209a30..218b35c 100644
--- a/experiments/hilbert_qr.cpp
+++ b/experiments/hilbert_qr.cpp
@@ -1,9 +1,3 @@
-// Experiment: QR methods on Hilbert matrices
-//
-// The Hilbert matrix H[i][j] = 1/(i+j+1) is the canonical ill-conditioned
-// dense matrix.  Its condition number grows roughly as (3.5 * e)^n / sqrt(n),
-// reaching ~10^13 at n=10 and ~10^18 at n=14.
-//
 // We compare classical GS, modified GS, and Householder QR on:
 //   - reconstruction error   ||A - QR||_F
 //   - orthogonality error    ||Q^T Q - I||_F
@@ -68,9 +62,6 @@ double orthogonality_error(const QRResult& qr) {
     return std::sqrt(err);
 }
 
-// ---------------------------------------------------------------------------
-// Timing
-// ---------------------------------------------------------------------------
 
 using Clock = std::chrono::high_resolution_clock;
 using Seconds = std::chrono::duration<double>;
@@ -88,10 +79,6 @@ double min_time(Fn fn, int trials = 5) {
     return best;
 }
 
-// ---------------------------------------------------------------------------
-// Run one method on one size, return {recon, ortho, time} or nullopt on failure
-// ---------------------------------------------------------------------------
-
 using QRFn = std::function<QRResult(const Matrix&)>;
 
 struct Result { double recon, ortho, time_s; };
@@ -110,10 +97,6 @@ std::optional<Result> measure(const Matrix& A, QRFn fn) {
     }
 }
 
-// ---------------------------------------------------------------------------
-// Pretty printing
-// ---------------------------------------------------------------------------
-
 void print_row(const std::string& method, std::optional<Result> r) {
     std::cout << std::left << std::setw(16) << method;
     if (!r) {
diff --git a/experiments/matmul.cpp b/experiments/matmul.cpp
new file mode 100644
index 0000000..8ef0ac5
--- /dev/null
+++ b/experiments/matmul.cpp
@@ -0,0 +1,149 @@
+#include "linalg_error.hpp"
+#include "matrix.hpp"
+
+#include <chrono>
+#include <cstddef>
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+#if !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX512F__)
+#   define MATMUL_BACKEND "AVX512"
+#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX2__)
+#   define MATMUL_BACKEND "AVX2"
+#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && defined(__AVX__)
+#   define MATMUL_BACKEND "AVX"
+#elif !defined(LINEAR_ALGEBRA_FORCE_SCALAR_MATMUL) && \
+      defined(__ARM_NEON) && defined(__aarch64__) && \
+      defined(__ARM_FEATURE_FP64_VECTOR_ARITHMETIC)
+#   define MATMUL_BACKEND "NEON"
+#else
+#   define MATMUL_BACKEND "scalar"
+#endif
+
+using linalg::Matrix;
+using Clock   = std::chrono::high_resolution_clock;
+using Seconds = std::chrono::duration<double>;
+
+
+Matrix naive_matmul(const Matrix& lhs, const Matrix& rhs) {
+    if (lhs.cols() != rhs.rows()) {
+        throw linalg::DimensionMismatchError(
+            "naive_matmul: lhs.cols() != rhs.rows()");
+    }
+    const std::size_t m = lhs.rows();
+    const std::size_t n = rhs.cols();
+    const std::size_t k = lhs.cols();
+
+    Matrix result(m, n, 0.0);
+
+    const double* A = lhs.data();
+    const double* B = rhs.data();
+    double*       C = result.data();
+
+    for (std::size_t i = 0; i < m; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            for (std::size_t p = 0; p < k; ++p)
+                C[i * n + j] += A[i * k + p] * B[p * n + j];
+
+    return result;
+}
+
+// ---------------------------------------------------------------------------
+// Helpers
+// ---------------------------------------------------------------------------
+
+// Fill an n×n matrix with a deterministic pattern so the compiler cannot
+// optimise multiplications away.
+Matrix make_matrix(std::size_t n) {
+    Matrix M(n, n);
+    const double inv = 1.0 / static_cast<double>(n + 1);
+    for (std::size_t i = 0; i < n; ++i)
+        for (std::size_t j = 0; j < n; ++j)
+            M(i, j) = static_cast<double>(i + j + 1) * inv;
+    return M;
+}
+
+template <typename Fn>
+double min_time_s(Fn fn, int trials) {
+    double best = 1e30;
+    for (int t = 0; t < trials; ++t) {
+        const auto t0 = Clock::now();
+        fn();
+        const auto t1 = Clock::now();
+        const double elapsed = Seconds(t1 - t0).count();
+        if (elapsed < best) best = elapsed;
+    }
+    return best;
+}
+
+double flops(std::size_t n) {
+    const double nd = static_cast<double>(n);
+    return 2.0 * nd * nd * nd;
+}
+
+
+int main() {
+    const std::vector<std::size_t> sizes = {
+        8, 16, 32, 64, 128, 256, 512
+    };
+
+    constexpr std::size_t small_threshold = 128;
+    constexpr int         trials_small    = 9;
+    constexpr int         trials_large    = 3;
+
+    std::cout << std::string(72, '*') << "\n";
+    std::cout << "  Matmul benchmark: naive (ijk) vs SIMD (" MATMUL_BACKEND
+                 ") + transpose\n";
+    std::cout << "  C = A * B,  A and B both n×n\n";
+    std::cout << std::string(72, '*') << "\n\n";
+
+    // Column header.
+    std::cout << std::left
+              << std::setw(6)  << "n"
+              << std::setw(14) << "naive ms"
+              << std::setw(14) << "naive GFLOP/s"
+              << std::setw(14) << "SIMD ms"
+              << std::setw(14) << "SIMD GFLOP/s"
+              << std::setw(10) << "speedup"
+              << "\n";
+    std::cout << std::string(72, '-') << "\n";
+
+    for (const std::size_t n : sizes) {
+        const Matrix A = make_matrix(n);
+        const Matrix B = make_matrix(n);
+
+        const int trials = (n <= small_threshold) ? trials_small : trials_large;
+
+        volatile double sink_naive = naive_matmul(A, B)(0, 0);
+        volatile double sink_simd  = (A * B)(0, 0);
+        (void)sink_naive;
+        (void)sink_simd;
+
+        const double t_naive = min_time_s([&] { (void)naive_matmul(A, B); }, trials);
+        const double t_simd  = min_time_s([&] { (void)(A * B); },           trials);
+
+        const double fp       = flops(n);
+        const double gf_naive = fp / t_naive / 1e9;
+        const double gf_simd  = fp / t_simd  / 1e9;
+        const double speedup  = t_naive / t_simd;
+
+        std::cout << std::left  << std::setw(6) << n
+                  << std::fixed << std::setprecision(3)
+                  << std::setw(14) << t_naive * 1e3
+                  << std::setprecision(2)
+                  << std::setw(14) << gf_naive
+                  << std::setprecision(3)
+                  << std::setw(14) << t_simd * 1e3
+                  << std::setprecision(2)
+                  << std::setw(14) << gf_simd
+                  << std::setprecision(2) << std::setw(10) << speedup
+                  << "x\n";
+    }
+
+    std::cout << "\n(each cell = minimum over "
+              << trials_small << " trials for n<=" << small_threshold
+              << ", " << trials_large << " trials for larger n)\n";
+
+    return 0;
+}
diff --git a/experiments/pivoting_vs_no_pivoting.cpp b/experiments/pivoting_vs_no_pivoting.cpp
index 0e71dc1..57281d8 100644
--- a/experiments/pivoting_vs_no_pivoting.cpp
+++ b/experiments/pivoting_vs_no_pivoting.cpp
@@ -1,8 +1,3 @@
-// Experiment: partial pivoting vs no-pivot LU
-//
-// Demonstrates why partial pivoting is essential for numerical stability.
-// Run the binary and inspect the residuals printed to stdout.
-
 #include "lu.hpp"
 #include "matrix.hpp"
 #include "norms.hpp"
@@ -23,7 +18,6 @@ using linalg::Vector;
 
 // ---------------------------------------------------------------------------
 // Local no-pivot LU for comparison only.
-// This is intentionally naive — it is here to show what breaks without pivoting.
 // ---------------------------------------------------------------------------
 
 struct NoPivotLU {
@@ -55,8 +49,6 @@ NoPivotLU lu_no_pivot(const Matrix& A, double tol = 1e-14) {
     return NoPivotLU{std::move(L), std::move(U), false, 0};
 }
 
-// Solve using a no-pivot LU (L unit lower triangular, U upper triangular).
-// If the factorization failed or U is numerically singular, returns nullopt.
 std::optional<Vector> solve_no_pivot(const NoPivotLU& f, const Vector& b) {
     if (f.failed) return std::nullopt;
     try {
@@ -178,8 +170,6 @@ void exp_random(std::size_t n = 8) {
 }
 
 // 2. Badly row-scaled matrix
-//    Rows differ in magnitude by ~10^14.  Without pivoting, tiny early pivots
-//    amplify round-off; with pivoting, the large-row is selected first.
 void exp_badly_scaled() {
     const Matrix A{
         {1e-14, 1.0,   2.0  },
@@ -191,21 +181,15 @@ void exp_badly_scaled() {
 }
 
 // 3. Classic pathological example for no-pivot LU.
-//    With epsilon = 1e-15, no-pivot computes a huge multiplier (1/epsilon),
-//    which causes catastrophic cancellation in the updated rows.
-//    With pivoting, we swap first and the multiplier is bounded by 1.
 void exp_epsilon_pathology() {
     constexpr double eps = 1e-15;
     const Matrix A{{eps, 1.0}, {1.0, 2.0}};
-    //  True solution of [eps 1; 1 2] * x = [1+eps; 3] is x = [1; 1].
     const Vector b{1.0 + eps, 3.0};
     run_case("Epsilon pathology [[1e-15,1],[1,2]] (classic)", A, b);
     std::cout << "  Note: exact solution is x = [1, 1]\n";
 }
 
 // 4. Matrix where no-pivot LU diverges visibly on a 4x4 example.
-//    The first pivot is small (0.001) but rows below have entries ~1000.
-//    No pivot causes multipliers of magnitude 10^6, annihilating subdiagonal info.
 void exp_amplified_multiplier() {
     const Matrix A{
         {0.001, 1.0,   0.0,   0.0  },
diff --git a/include/linalg_error.hpp b/include/linalg_error.hpp
index 25c48ad..ac9b456 100644
--- a/include/linalg_error.hpp
+++ b/include/linalg_error.hpp
@@ -22,4 +22,10 @@ public:
         : LinAlgError(message) {}
 };
 
+class NonConvergenceError : public LinAlgError {
+public:
+    explicit NonConvergenceError(const std::string& message)
+        : LinAlgError(message) {}
+};
+
 }  // namespace linalg
diff --git a/include/qr_iteration.hpp b/include/qr_iteration.hpp
new file mode 100644
index 0000000..08e17fa
--- /dev/null
+++ b/include/qr_iteration.hpp
@@ -0,0 +1,82 @@
+#pragma once
+
+// QR iteration for eigenvalue computation
+//
+// Refs:
+//   Trefethen & Bau, "Numerical Linear Algebra" (T&B)
+//     Lecture 25 — Eigenvalue algorithms
+//     Lecture 26 — Schur factorisation
+//     Lecture 28 — The QR algorithm (unshifted)
+//     Lecture 29 — The QR algorithm with shifts
+//   Golub & Van Loan, "Matrix Computations" 4th ed. (GVL)
+//     §7.3  — The Unshifted QR Algorithm
+//     §7.4  — The Shifted QR Algorithm
+//     §7.4.2 — Wilkinson shift
+
+#include <vector>
+
+#include "linalg_error.hpp"
+#include "matrix.hpp"
+#include "vector.hpp"
+
+namespace linalg {
+
+// ---------------------------------------------------------------------------
+// Options
+// ---------------------------------------------------------------------------
+
+// All defaults are consistent with the recommendations in T&B Lecture 28.
+struct QRIterationOptions {
+    // Convergence threshold.  Iteration halts once the Frobenius norm of the
+    // strict lower triangle of A_k falls below this value.
+    // Ref: T&B §28; GVL §7.3.
+    double tolerance = 1e-10;
+
+    int max_iterations = 1000;
+
+    // When true, the Frobenius norm of the strict lower triangle is recorded
+    // after every QR step and returned in QRIterationResult::convergence_history.
+    bool track_convergence = false;
+};
+
+
+struct QRIterationResult {
+    // Real and imaginary parts of the n eigenvalues.
+    // For symmetric inputs all imaginary parts are zero.
+    // Complex-conjugate pairs from 2×2 Schur blocks appear as ±imag entries.
+    // Both vectors always have length n (the matrix dimension).
+    Vector eigenvalues_real;
+    Vector eigenvalues_imag;
+
+    // Total number of QR steps performed before convergence or max_iterations.
+    int iterations = 0;
+
+    // Populated only when QRIterationOptions::track_convergence is true.
+    // Entry k is ||lower(A_k)||_F after the k-th QR step.
+    // std::vector is used here because linalg::Vector has no push_back;
+    // convergence_history is a plain time-series container, not a math object.
+    std::vector<double> convergence_history;
+};
+
+// Algorithm (T&B Algorithm 28.1):
+//
+//   A_0 = A
+//   for k = 1, 2, ...:
+//       factor  A_{k-1} = Q_k R_k   (Householder QR)
+//       set     A_k     = R_k Q_k   (orthogonal similarity: preserves eigenvalues)
+//
+// The iterates A_k converge to the real Schur form of A: a quasi-upper-
+// triangular matrix with 1×1 blocks (real eigenvalue) and 2×2 blocks
+// (complex-conjugate pair) on the diagonal.
+//
+// Convergence rate: linear.  Per-step factor ≈ |lambda_{j+1} / lambda_j|
+// for the off-diagonal entries linking eigenvalue clusters j and j+1.
+// (T&B Lecture 28, Theorem 28.2)
+//
+// Throws DimensionMismatchError  if A is not square.
+// Throws NonConvergenceError     if convergence is not achieved within
+//                                opts.max_iterations steps.
+[[nodiscard]] QRIterationResult eigenvalues_unshifted(const Matrix& A,
+                                                      QRIterationOptions opts = {});
+
+}  // namespace linalg
diff --git a/src/qr_iteration.cpp b/src/qr_iteration.cpp
new file mode 100644
index 0000000..f488a5f
--- /dev/null
+++ b/src/qr_iteration.cpp
@@ -0,0 +1,207 @@
+#include "qr_iteration.hpp"
+
+#include <cassert>
+#include <cmath>
+#include <sstream>
+
+#include "linalg_error.hpp"
+#include "matrix.hpp"
+#include "qr.hpp"
+#include "vector.hpp"
+
+// References used throughout this file:
+//   T&B  — Trefethen & Bau, "Numerical Linear Algebra"
+//   GVL  — Golub & Van Loan, "Matrix Computations" 4th ed.
+
+namespace linalg {
+
+namespace {
+
+// ---------------------------------------------------------------------------
+// Internal helpers
+// ---------------------------------------------------------------------------
+
+// Frobenius norm of the strict lower triangle of an n×n matrix.
+// This is the standard convergence diagnostic for QR iteration: as A_k
+// approaches the real Schur form, all entries below the main diagonal
+// (excluding 2×2 block sub-diagonals) tend to zero.
+//
+// ||lower(A)||_F = sqrt( sum_{i > j} A(i,j)^2 )
+//
+// Ref: T&B §28; used as the convergence criterion in Algorithm 28.1.
+double lower_triangle_norm(const Matrix& A) {
+    const std::size_t n = A.rows();
+    double s = 0.0;
+    for (std::size_t i = 1; i < n; ++i)      // row 1 .. n-1
+        for (std::size_t j = 0; j < i; ++j)  // col 0 .. i-1  (strict lower)
+            s += A(i, j) * A(i, j);
+    return std::sqrt(s);
+}
+
+// Extract eigenvalues from a quasi-upper-triangular matrix (real Schur form).
+//
+// Scans the diagonal from top-left to bottom-right.  At each position i:
+//   — |A(i+1, i)| < tol  → 1×1 block: real eigenvalue A(i,i), imag = 0.
+//   — otherwise           → 2×2 block [A(i..i+1, i..i+1)]: eigenvalues via
+//                           quadratic formula.  When the discriminant is
+//                           negative the result is a complex-conjugate pair,
+//                           stored as (re, +im) and (re, -im) in the real
+//                           and imaginary part Vectors.
+//
+// Fills positions 0..n-1 of `real_out` and `imag_out` (pre-sized to n).
+//
+// Ref: T&B Lecture 28; GVL §7.4.1.
+void extract_eigenvalues(const Matrix& T, double tol,
+                         Vector& real_out, Vector& imag_out) {
+    const std::size_t n = T.rows();
+    std::size_t out = 0;  // next write position in real_out / imag_out
+    std::size_t i   = 0;  // current scan position in T
+
+    while (i < n) {
+        const bool is_last   = (i + 1 == n);
+        const bool sub_small = is_last || (std::abs(T(i + 1, i)) < tol);
+
+        if (sub_small) {
+            // 1×1 block: real eigenvalue.
+            real_out[out] = T(i, i);
+            imag_out[out] = 0.0;
+            ++out;
+            ++i;
+        } else {
+            // 2×2 block:
+            //   | a  b |
+            //   | c  d |
+            // Characteristic polynomial: lambda^2 - (a+d)*lambda + (ad - bc) = 0.
+            // Discriminant: (a-d)^2 + 4*b*c.
+            // Ref: GVL §7.4.1.
+            const double a    = T(i,     i);
+            const double b    = T(i,     i + 1);
+            const double c    = T(i + 1, i);
+            const double d    = T(i + 1, i + 1);
+            const double tr   = a + d;
+            const double disc = (a - d) * (a - d) + 4.0 * b * c;
+
+            if (disc >= 0.0) {
+                // Real eigenvalues — unusual in converged real Schur form, but
+                // handled robustly in case the block didn't fully split.
+                const double sq   = std::sqrt(disc);
+                real_out[out]     = 0.5 * (tr + sq);
+                imag_out[out]     = 0.0;
+                real_out[out + 1] = 0.5 * (tr - sq);
+                imag_out[out + 1] = 0.0;
+            } else {
+                // Complex-conjugate pair: real part ± imaginary part.
+                const double re       = 0.5 * tr;
+                const double im       = 0.5 * std::sqrt(-disc);
+                real_out[out]         = re;
+                imag_out[out]         =  im;
+                real_out[out + 1]     = re;
+                imag_out[out + 1]     = -im;
+            }
+            out += 2;
+            i   += 2;
+        }
+    }
+
+    assert(out == n);
+}
+
+// Verify that A is square; throw DimensionMismatchError otherwise.
+void require_square(const Matrix& A, const char* fname) {
+    if (A.rows() != A.cols()) {
+        std::ostringstream oss;
+        oss << fname << ": requires a square matrix, got "
+            << A.rows() << "x" << A.cols();
+        throw DimensionMismatchError(oss.str());
+    }
+}
+
+}  // namespace
+
+// ---------------------------------------------------------------------------
+// Stage 1: Unshifted QR iteration
+// ---------------------------------------------------------------------------
+//
+// Each step performs an orthogonal similarity transformation:
+//   A_{k-1} = Q_k R_k   (Householder QR; backward-stable)
+//   A_k     = R_k Q_k   = Q_k^T A_{k-1} Q_k
+//
+// Similarity preserves eigenvalues (GVL §7.3.1, Theorem 7.3.1).
+// The iterates converge to the real Schur form: a quasi-upper-triangular
+// matrix whose 1×1 blocks give real eigenvalues and 2×2 blocks give
+// complex-conjugate pairs.
+//
+// Convergence rate: linear, with per-step reduction factor
+//   |lambda_{j+1} / lambda_j|  for the (j, j+1) coupling.
+// (T&B Lecture 28, Theorem 28.2; GVL §7.3.2)
+//
+// Each iteration costs O(n^3) due to full Householder QR; Hessenberg
+// reduction (Stage 3) reduces subsequent steps to O(n^2).
+
+QRIterationResult eigenvalues_unshifted(const Matrix& A,
+                                        QRIterationOptions opts) {
+    require_square(A, "eigenvalues_unshifted");
+    const std::size_t n = A.rows();
+
+    // Threshold for classifying a sub-diagonal entry as "zero" when reading
+    // eigenvalues out of the converged Schur form.  Using the same value as
+    // the convergence tolerance is appropriate; we only reach extraction once
+    // ||lower(A_k)||_F < opts.tolerance.  Ref: GVL §7.4.1.
+    const double extract_tol = opts.tolerance;
+
+    QRIterationResult result;
+    // Pre-size eigenvalue Vectors; they are always length n.
+    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));
+    }
+
+    // Handle the trivial 1×1 case immediately.
+    if (n == 1) {
+        result.eigenvalues_real[0] = A(0, 0);
+        return result;
+    }
+
+    // Working copy; becomes the quasi-upper-triangular Schur form A_k.
+    Matrix Ak = A;
+
+    for (int k = 0; k < opts.max_iterations; ++k) {
+        // --- QR step ---
+        // Factor A_{k-1} = Q R using backward-stable Householder reflections.
+        const QRResult qr = qr_householder(Ak);
+
+        // A_k = R Q  (orthogonal similarity: Q^T A_{k-1} Q)
+        Ak = qr.R * qr.Q;
+
+        // --- Convergence check ---
+        const double lower_norm = lower_triangle_norm(Ak);
+
+        if (opts.track_convergence) {
+            result.convergence_history.push_back(lower_norm);
+        }
+        ++result.iterations;
+
+        if (lower_norm < opts.tolerance) {
+            extract_eigenvalues(Ak, extract_tol,
+                                result.eigenvalues_real,
+                                result.eigenvalues_imag);
+            return result;
+        }
+    }
+
+    // Maximum iterations reached without convergence — fail loudly.
+    // Possible causes: eigenvalues too close in magnitude, or complex
+    // eigenvalue pairs that require a double shift (see Stage 2).
+    std::ostringstream oss;
+    oss << "eigenvalues_unshifted: did not converge in "
+        << opts.max_iterations << " iterations "
+        << "(final ||lower(A_k)||_F = " << lower_triangle_norm(Ak)
+        << ", tolerance = " << opts.tolerance << "). "
+        << "Try eigenvalues_shifted (Stage 2) or increase max_iterations.";
+    throw NonConvergenceError(oss.str());
+}
+
+}  // namespace linalg
diff --git a/tests/test_qr_iteration.cpp b/tests/test_qr_iteration.cpp
new file mode 100644
index 0000000..c994c95
--- /dev/null
+++ b/tests/test_qr_iteration.cpp
@@ -0,0 +1,293 @@
+// Tests for qr_iteration.hpp / qr_iteration.cpp
+//
+// Stage 1: Unshifted QR iteration.
+//
+// All Stage 1 tests use symmetric matrices (only real eigenvalues) because
+// the unshifted algorithm converges to upper-triangular form — not merely
+// quasi-upper-triangular — only when all eigenvalues are real.  A matrix
+// with a complex-conjugate pair would stall: its 2×2 Schur block keeps a
+// non-negligible subdiagonal entry indefinitely, so ||lower(A_k)||_F never
+// falls below the tolerance.  Proper handling of complex pairs requires the
+// double-shift strategy introduced in Stage 2.
+//
+// Refs: T&B Lecture 28; GVL §7.3–7.4.
+
+#include "linalg_error.hpp"
+#include "matrix.hpp"
+#include "qr_iteration.hpp"
+#include "vector.hpp"
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <iostream>
+#include <utility>
+#include <vector>
+
+using linalg::Matrix;
+using linalg::NonConvergenceError;
+using linalg::QRIterationOptions;
+using linalg::QRIterationResult;
+using linalg::Vector;
+
+// ---------------------------------------------------------------------------
+// Test helpers
+// ---------------------------------------------------------------------------
+
+namespace {
+
+// Sort (real, imag) eigenvalue pairs by real part (ascending), then by imag.
+// Returns a std::vector<std::pair<double,double>> — a plain container of
+// pairs, not a math vector.
+using EigPairs = std::vector<std::pair<double, double>>;
+
+EigPairs to_pairs(const Vector& real_v, const Vector& imag_v) {
+    EigPairs out;
+    out.reserve(real_v.size());
+    for (std::size_t i = 0; i < real_v.size(); ++i)
+        out.emplace_back(real_v[i], imag_v[i]);
+    std::sort(out.begin(), out.end(),
+              [](const std::pair<double, double>& a,
+                 const std::pair<double, double>& b) {
+                  return a.first != b.first ? a.first < b.first
+                                            : a.second < b.second;
+              });
+    return out;
+}
+
+// Return true when every computed eigenvalue is within `tol` of the
+// corresponding expected eigenvalue (after sorting both sets).
+// `expected` is a plain std::vector of (real, imag) pairs used as test data.
+bool eigs_match(const Vector& computed_real, const Vector& computed_imag,
+                const EigPairs& expected, double tol) {
+    if (computed_real.size() != expected.size()) return false;
+    const EigPairs computed = to_pairs(computed_real, computed_imag);
+    EigPairs exp_sorted      = expected;
+    std::sort(exp_sorted.begin(), exp_sorted.end(),
+              [](const std::pair<double, double>& a,
+                 const std::pair<double, double>& b) {
+                  return a.first != b.first ? a.first < b.first
+                                            : a.second < b.second;
+              });
+    for (std::size_t i = 0; i < computed.size(); ++i) {
+        const double dr = computed[i].first  - exp_sorted[i].first;
+        const double di = computed[i].second - exp_sorted[i].second;
+        if (std::sqrt(dr * dr + di * di) > tol) return false;
+    }
+    return true;
+}
+
+}  // namespace
+
+// ---------------------------------------------------------------------------
+// Test 1: 2×2 symmetric matrix with known eigenvalues
+// ---------------------------------------------------------------------------
+//
+// A = | 2  1 |   is symmetric positive definite.
+//     | 1  2 |
+//
+// Characteristic polynomial: (2-λ)^2 - 1 = 0  →  λ = 1, 3.
+// Eigenvectors: [1,-1]/√2 (λ=1) and [1,1]/√2 (λ=3).
+//
+// The unshifted iteration converges at rate |λ_1/λ_2| = 1/3 per step,
+// so only a handful of iterations are needed.
+// Ref: T&B Theorem 28.2.
+
+TEST_CASE("QR iteration (unshifted): 2x2 symmetric known eigenvalues",
+          "[qr_iteration][stage1]") {
+    const Matrix A{
+        {2.0, 1.0},
+        {1.0, 2.0}
+    };
+
+    const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+
+    REQUIRE(res.eigenvalues_real.size() == 2);
+    REQUIRE(res.eigenvalues_imag.size() == 2);
+    REQUIRE(res.iterations > 0);
+
+    // All eigenvalues of a symmetric matrix must be real.
+    CHECK(std::abs(res.eigenvalues_imag[0]) < 1e-8);
+    CHECK(std::abs(res.eigenvalues_imag[1]) < 1e-8);
+
+    const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}};
+    CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
+}
+
+// ---------------------------------------------------------------------------
+// Test 2: 4×4 symmetric tridiagonal matrix — reference eigenvalues
+// ---------------------------------------------------------------------------
+//
+// The n×n symmetric tridiagonal matrix with 2 on the diagonal and -1 on the
+// first super- and sub-diagonals has known eigenvalues (discrete Laplacian):
+//
+//   λ_k = 2 - 2 cos(k π / (n+1)),   k = 1, …, n
+//
+// Ref: Golub & Van Loan §4.4.2 (discrete sine transform).
+//
+// For n = 4:
+//   λ_1 = 2 - 2 cos(π/5)   ≈ 0.3820
+//   λ_2 = 2 - 2 cos(2π/5)  ≈ 1.3820
+//   λ_3 = 2 - 2 cos(3π/5)  ≈ 2.6180
+//   λ_4 = 2 - 2 cos(4π/5)  ≈ 3.6180
+
+TEST_CASE("QR iteration (unshifted): 4x4 symmetric tridiagonal",
+          "[qr_iteration][stage1]") {
+    const Matrix A{
+        { 2.0, -1.0,  0.0,  0.0},
+        {-1.0,  2.0, -1.0,  0.0},
+        { 0.0, -1.0,  2.0, -1.0},
+        { 0.0,  0.0, -1.0,  2.0}
+    };
+
+    const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+
+    REQUIRE(res.eigenvalues_real.size() == 4);
+    REQUIRE(res.eigenvalues_imag.size() == 4);
+
+    // All eigenvalues of a symmetric matrix must be real.
+    for (std::size_t k = 0; k < 4; ++k)
+        CHECK(std::abs(res.eigenvalues_imag[k]) < 1e-8);
+
+    // Compare against the closed-form reference.
+    constexpr double pi = 3.14159265358979323846;
+    const EigPairs expected = {
+        {2.0 - 2.0 * std::cos(      pi / 5.0), 0.0},
+        {2.0 - 2.0 * std::cos(2.0 * pi / 5.0), 0.0},
+        {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0},
+        {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0}
+    };
+    CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
+}
+
+// ---------------------------------------------------------------------------
+// Test 3: 5×5 symmetric tridiagonal — convergence history
+// ---------------------------------------------------------------------------
+//
+// Uses a 5×5 symmetric tridiagonal (discrete Laplacian) to guarantee all
+// real eigenvalues and predictable linear convergence.  The Frobenius norm
+// of the strict lower triangle is printed at every step so the convergence
+// rate can be observed directly.
+//
+// Expected behaviour: ||lower(A_k)||_F decreases geometrically each step
+// (linear convergence), with ratio ≈ max_j |λ_{j+1}/λ_j|.
+// Ref: T&B Theorem 28.2.
+//
+// 5×5 tridiagonal eigenvalues: λ_k = 2 - 2cos(kπ/6), k = 1..5.
+
+TEST_CASE("QR iteration (unshifted): 5x5 convergence history",
+          "[qr_iteration][stage1]") {
+    const Matrix A{
+        { 2.0, -1.0,  0.0,  0.0,  0.0},
+        {-1.0,  2.0, -1.0,  0.0,  0.0},
+        { 0.0, -1.0,  2.0, -1.0,  0.0},
+        { 0.0,  0.0, -1.0,  2.0, -1.0},
+        { 0.0,  0.0,  0.0, -1.0,  2.0}
+    };
+
+    QRIterationOptions opts;
+    opts.track_convergence = true;
+
+    const QRIterationResult res = linalg::eigenvalues_unshifted(A, opts);
+
+    REQUIRE_FALSE(res.convergence_history.empty());
+    REQUIRE(res.eigenvalues_real.size() == 5);
+
+    // Print convergence history so the linear rate is visible.
+    std::cout << "\n=== Stage 1: Unshifted QR — 5x5 convergence history ===\n";
+    std::cout << "  Converged in " << res.iterations << " iteration(s)\n";
+    for (std::size_t k = 0; k < res.convergence_history.size(); ++k) {
+        std::cout << "  iter " << (k + 1)
+                  << ": ||lower(A_k)||_F = "
+                  << res.convergence_history[k] << "\n";
+    }
+    std::cout << "=======================================================\n";
+
+    // The final recorded norm must be below the default tolerance.
+    CHECK(res.convergence_history.back() < opts.tolerance);
+}
+
+// ---------------------------------------------------------------------------
+// Test 4: Eigenvalue residuals below 1e-8
+// ---------------------------------------------------------------------------
+//
+// For several symmetric matrices with analytically known eigenvalues, verify
+// that every computed eigenvalue is within 1e-8 of its expected value.
+//
+// Residual means the absolute error |λ_computed - λ_exact| (eigenvalue
+// accuracy), not a matrix residual ||A x - λ x||, which would require
+// eigenvectors unavailable in Stage 1.
+
+TEST_CASE("QR iteration (unshifted): residuals below 1e-8",
+          "[qr_iteration][stage1]") {
+
+    SECTION("2x2: eigenvalues 1 and 3") {
+        const Matrix A{{2.0, 1.0}, {1.0, 2.0}};
+        const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+        const EigPairs expected = {{1.0, 0.0}, {3.0, 0.0}};
+        CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
+    }
+
+    SECTION("3x3 diagonal: eigenvalues 1, 4, 9") {
+        // Diagonal matrix — already in Schur form; converges in one step.
+        const Matrix D{
+            {1.0, 0.0, 0.0},
+            {0.0, 4.0, 0.0},
+            {0.0, 0.0, 9.0}
+        };
+        const QRIterationResult res = linalg::eigenvalues_unshifted(D);
+        const EigPairs expected = {{1.0, 0.0}, {4.0, 0.0}, {9.0, 0.0}};
+        CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
+    }
+
+    SECTION("4x4 tridiagonal: closed-form eigenvalues") {
+        const Matrix A{
+            { 2.0, -1.0,  0.0,  0.0},
+            {-1.0,  2.0, -1.0,  0.0},
+            { 0.0, -1.0,  2.0, -1.0},
+            { 0.0,  0.0, -1.0,  2.0}
+        };
+        constexpr double pi = 3.14159265358979323846;
+        const EigPairs expected = {
+            {2.0 - 2.0 * std::cos(      pi / 5.0), 0.0},
+            {2.0 - 2.0 * std::cos(2.0 * pi / 5.0), 0.0},
+            {2.0 - 2.0 * std::cos(3.0 * pi / 5.0), 0.0},
+            {2.0 - 2.0 * std::cos(4.0 * pi / 5.0), 0.0}
+        };
+        const QRIterationResult res = linalg::eigenvalues_unshifted(A);
+        CHECK(eigs_match(res.eigenvalues_real, res.eigenvalues_imag, expected, 1e-8));
+    }
+
+    SECTION("5x5 identity: all eigenvalues == 1") {
+        const Matrix I = Matrix::identity(5);
+        const QRIterationResult res = linalg::eigenvalues_unshifted(I);
+        REQUIRE(res.eigenvalues_real.size() == 5);
+        for (std::size_t k = 0; k < 5; ++k) {
+            CHECK(std::abs(res.eigenvalues_real[k] - 1.0) < 1e-8);
+            CHECK(std::abs(res.eigenvalues_imag[k])       < 1e-8);
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Failure cases
+// ---------------------------------------------------------------------------
+
+TEST_CASE("QR iteration (unshifted): non-square matrix throws",
+          "[qr_iteration][stage1]") {
+    const Matrix A(3, 4);  // non-square
+    CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A),
+                    linalg::DimensionMismatchError);
+}
+
+TEST_CASE("QR iteration (unshifted): max_iterations exceeded throws",
+          "[qr_iteration][stage1]") {
+    // Cap at zero iterations — any non-trivial matrix fails immediately.
+    const Matrix A{{2.0, 1.0}, {1.0, 2.0}};
+    QRIterationOptions opts;
+    opts.max_iterations = 0;
+    CHECK_THROWS_AS(linalg::eigenvalues_unshifted(A, opts), NonConvergenceError);
+}