Add QR iteration and cleanup

1 files changed, 82 insertions, 0 deletions

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