Add QR eigen library: implicit Wilkinson-shifted QR for symmetric matrices

- src/QR.hpp / src/QR.cpp: fully templated QR::EigenQR (N >= 2), no heap
  allocation (3*N^2 float working buffers on stack). Givens tridiagonalization
  (bottom-up) + implicit Wilkinson-shifted QR with bulge chasing, relative
  deflation, exact-zero peeling, closed-form 2x2 termination.
- Matrix::EigenQR now delegates to QR::EigenQR (old unshifted body removed);
  eigenvalues sorted descending, eigenvectors in columns of the output.
- unit-tests/qr-build-blocks-tests.cpp: 8 building-block test cases
  (215 assertions) with scipy/numpy references.
- unit-tests/qr-reference-values.py: numpy/scipy reference generator mirroring
  every building block and the full pipeline (eigh, n=3..8).
- CMake: new 'qr' static library; Matrix links against it;
  qr-build-blocks-tests target enabled.

src/CMakeLists.txt

@@ -41,23 +41,6 @@ target_link_libraries(vector-3d
     PRIVATE
 )
 
-# Matrix
-add_library(matrix 
-    STATIC
-    Matrix.cpp
-)
-
-target_link_libraries(matrix
-    PUBLIC
-    vector-3d-intf
-    PRIVATE
-)
-
-set_target_properties(matrix
-    PROPERTIES
-    LINKER_LANGUAGE CXX
-)
-
 # SVD
 add_library(svd 
     STATIC
@@ -73,4 +56,40 @@ target_link_libraries(svd
 set_target_properties(svd
     PROPERTIES
     LINKER_LANGUAGE CXX
+)
+
+# QR (eigenvalues/eigenvectors via implicit shifted QR iteration)
+add_library(qr 
+    STATIC
+    QR.cpp
+)
+
+target_link_libraries(qr
+    PUBLIC
+    vector-3d-intf
+    PRIVATE
+)
+
+set_target_properties(qr
+    PROPERTIES
+    LINKER_LANGUAGE CXX
+)
+
+# Matrix
+add_library(matrix 
+    STATIC
+    Matrix.cpp
+)
+
+target_link_libraries(matrix
+    PUBLIC
+    vector-3d-intf
+    PRIVATE
+    svd
+    qr
+)
+
+set_target_properties(matrix
+    PROPERTIES
+    LINKER_LANGUAGE CXX
 )

src/Matrix.cpp

@@ -5,6 +5,21 @@
 #include "Matrix.hpp"
 #endif
 
+// Forward-declare QR::EigenQR so the Matrix::EigenQR implementation below can
+// call it even when Matrix.cpp is pulled in through QR.hpp's own include chain
+// (QR.cpp -> QR.hpp -> Matrix.hpp -> Matrix.cpp), where the QR namespace has
+// not been declared yet at this point. If we are not already inside that
+// chain, pull in the full QR library so its template definition is available.
+namespace QR {
+template <uint8_t N>
+void EigenQR(Matrix<N, N> &matrixToDecompose, Matrix<N, N> &eigenVectors,
+             Matrix<N, 1> &eigenValues, uint32_t maxIterations,
+             float tolerance);
+}
+#ifndef QR_H_
+#include "QR.hpp"
+#endif
+
 #ifdef MATRIX_H_ // since the .cpp file has to be included by the .hpp file this
                  // will evaluate to true
 #include "Matrix.hpp"
@@ -571,37 +586,10 @@ void Matrix<rows, columns>::EigenQR(Matrix<rows, rows> &eigenVectors,
   static_assert(rows > 1, "Matrix size must be > 1 for QR iteration");
   static_assert(rows == columns, "Matrix size must be square for QR iteration");
 
-  Matrix<rows, rows> Ak = *this; // Copy original matrix
+  // Delegate to the QR library: implicit shifted QR iteration with
-  Matrix<rows, rows> QQ{Matrix<rows, rows>::Identity()};
+  // Wilkinson shift (see src/QR.hpp for the algorithm and conventions).
-  Matrix<rows, rows> shift{0};
+  Matrix<rows, rows> A = *this; // QR::EigenQR does not modify its input
-
+  QR::EigenQR(A, eigenVectors, eigenValues, maxIterations, tolerance);
-  for (uint32_t iter = 0; iter < maxIterations; ++iter) {
-    Matrix<rows, rows> Q, R;
-
-    // // QR shift lets us "attack" the first diagonal to speed up the algorithm
-    // shift = Matrix<rows, rows>::Identity() * Ak[rows - 1][rows - 1];
-    (Ak - shift).QRDecomposition(Q, R);
-    Ak = R * Q + shift;
-    QQ = QQ * Q;
-
-    // Check convergence: off-diagonal norm
-    float offDiagSum = 0.0f;
-    for (uint32_t row = 1; row < rows; row++) {
-      for (uint32_t column = 0; column < row; column++) {
-        offDiagSum += fabs(Ak[row][column]);
-      }
-    }
-
-    if (offDiagSum < tolerance) {
-      break;
-    }
-  }
-
-  // Diagonal elements are the eigenvalues
-  for (uint8_t i = 0; i < rows; i++) {
-    eigenValues[i][0] = Ak[i][i];
-  }
-  eigenVectors = QQ;
 }
 
 #endif // MATRIX_H_

src/Matrix.hpp

@@ -5,7 +5,6 @@
 #include <string>
 #include <type_traits>
 
-// TODO: Add a function to calculate eigenvalues/vectors
 // TODO: Add a function to compute RREF
 // TODO: Add a function for SVD decomposition
 // TODO: Add a function for LQ decomposition
@@ -234,12 +233,19 @@ public:
                        Matrix<columns, columns> &R) const;
 
   /**
-   * @brief Uses QR decomposition to efficiently calculate the eigenvectors
+   * @brief Calculates the eigenvectors and values of this matrix using the
-   * and values of this matrix
+   * implicit shifted QR iteration (Wilkinson shift, Givens bulge chasing);
-   * @param eigenVectors a buffer that will contain the eigenvectors fo this
+   * see src/QR.hpp in the QR library for the full algorithm.
-   * matrix
+   * @note For a matrix larger than 2x2 the matrix MUST be symmetric.
-   * @param eigenValues a buffer that will contain the eigenValues fo this
+   *       A general (nonsymmetric) 2x2 is handled via the closed-form
-   * matrix
+   *       solution.
+   * @note The eigenvalues come out sorted DESCENDING (largest first); the
+   *       eigenvector columns are swapped to match. Eigenvector signs are
+   *       arbitrary.
+   * @param eigenVectors a buffer that will contain the eigenvectors of this
+   * matrix in its columns (column i pairs with eigenValues[i])
+   * @param eigenValues a buffer that will contain the eigenvalues of this
+   * matrix, sorted descending
    * @param maxIterations the number of iterations to perform before giving
    * up on reaching the given tolerance
    * @param tolerance the level of accuracy to obtain before stopping.

src/QR.cpp

@@ -0,0 +1,416 @@
+// This #ifndef section makes clangd happy so that it can properly do type hints
+// in this file
+#ifndef QR_H_
+#define QR_H_
+#include "QR.hpp"
+#endif
+
+#ifdef QR_H_ // since the .cpp file has to be included by the .hpp file this
+             // will evaluate to true
+#include "QR.hpp"
+#include <cmath>
+#include <cstdint>
+
+namespace QR {
+
+// ============================================================================
+// QR Building Block Implementations (fully templated, heap-free)
+// ============================================================================
+
+/**
+ * GivensRotation: R * (a, b)^T = (r, 0)^T with R = [[c, s], [-s, c]],
+ * r = +hypot(a, b), c = a/r, s = b/r.
+ */
+static void GivensRotation(float a, float b, float &c, float &s) {
+  float r = sqrtf(a * a + b * b);
+  if (r == 0.0f) {
+    c = 1.0f;
+    s = 0.0f;
+    return;
+  }
+  c = a / r;
+  s = b / r;
+}
+
+/**
+ * ApplyRotationBothSides: A <- G A G^T (similarity transform) with
+ * G = [[c, s], [-s, c]] on the (i, i+1) block, i.e. G is the ZEROING
+ * rotation G*(x, y)^T = (r, 0)^T (the orientation used by the implicit QR
+ * chase: A = Q R with Q = G^T gives the next iterate R Q = G A G^T).
+ * With (c, s) = GivensRotation(A[i][i], A[i+1][i]) this zeroes
+ * A[i+1][i] after the LEFT multiplication; the right multiplication then
+ * chases the bulge along the superdiagonal (tridiagonal chase).
+ *
+ * A must be symmetric on entry; the result stays symmetric, so both
+ * triangles are written.
+ *
+ * Block updates (with a00 = A[i][i], a01 = A[i][i+1], a11 = A[i+1][i+1]):
+ *   A[i][i]     = c^2 a00 + 2 c s a01 + s^2 a11
+ *   A[i][i+1]   = (c^2 - s^2) a01 + c s (a11 - a00)
+ *   A[i+1][i+1] = s^2 a00 - 2 c s a01 + c^2 a11
+ * Off-block updates (uniform for both sides, since the left factor G and
+ * the right factor G^T mix each side with the pattern (a, b) -> (c a + s b,
+ * -s a + c b) after transposition):
+ *   for j not in {i, i+1}:
+ *     A[i][j]   = A[j][i]   = c A[i][j]   + s A[i+1][j]
+ *     A[i+1][j] = A[j][i+1] = -s A[i][j]  + c A[i+1][j]
+ */
+template <uint8_t N>
+static void ApplyRotationBothSides(Matrix<N, N> &A, uint8_t i, float c,
+                                   float s) {
+  float a00 = A.Get(i, i);
+  float a01 = A.Get(i, i + 1);
+  float a11 = A.Get(i + 1, i + 1);
+  float c2 = c * c;
+  float s2 = s * s;
+  float cs = c * s;
+
+  A[i][i] = c2 * a00 + 2.0f * cs * a01 + s2 * a11;
+  A[i][i + 1] = (c2 - s2) * a01 + cs * (a11 - a00);
+  A[i + 1][i + 1] = s2 * a00 - 2.0f * cs * a01 + c2 * a11;
+  A[i + 1][i] = A[i][i + 1]; // keep both triangles in sync
+
+  for (uint8_t j = 0; j < N; ++j) {
+    if (j == i || j == i + 1)
+      continue;
+    float x = A.Get(i, j);
+    float y = A.Get(i + 1, j);
+    A[i][j] = c * x + s * y;
+    A[j][i] = A[i][j];
+    A[i + 1][j] = -s * x + c * y;
+    A[j][i + 1] = A[i + 1][j];
+  }
+}
+
+/**
+ * ApplyRotationToVectors: V <- V G^T with G = [[c, s], [-s, c]] on columns
+ * (i, i+1), applied to every row. G^T = [[c, -s], [s, c]], so
+ *   V[r][i]     <- c V[r][i]   + s V[r][i+1]
+ *   V[r][i+1]   <- -s V[r][i]  + c V[r][i+1]
+ *
+ * Convention pairing: if A evolves as A <- G A G^T (ApplyRotationBothSides
+ * with the SAME c, s), then V accumulates V <- V G^T. With V0 = I the
+ * invariant A0 = V A V^T is preserved at every step, so at convergence
+ * A0 = V D V^T and the columns of V are the eigenvectors. (Rationale:
+ * each chase step is A <- R Q with R = G A the upper-triangular factor and
+ * Q = G^T the orthogonal factor of A = Q R, so A = G^T A' G and the
+ * orthogonal factors multiply as G1^T G2^T ... in application order.)
+ */
+template <uint8_t N>
+static void ApplyRotationToVectors(Matrix<N, N> &V, uint8_t i, float c,
+                                   float s) {
+  for (uint8_t r = 0; r < N; ++r) {
+    float x = V.Get(r, i);
+    float y = V.Get(r, i + 1);
+    V[r][i] = c * x + s * y;
+    V[r][i + 1] = -s * x + c * y;
+  }
+}
+
+/**
+ * WilkinsonShift: eigenvalue of [[a, b], [b, d]] closest to d.
+ * mu = (a+d)/2 - sign(a-d) * sqrt(((a-d)/2)^2 + b^2), sign(0) = +1.
+ */
+static float WilkinsonShift(float a, float b, float d) {
+  float delta = 0.5f * (a - d);
+  float spread = sqrtf(delta * delta + b * b);
+  return 0.5f * (a + d) - (delta >= 0.0f ? spread : -spread);
+}
+
+/**
+ * Solve2x2Eigen: closed-form eigen-decomposition of the 2x2 block at
+ * (lo, lo+1). Works for symmetric blocks and for general 2x2 blocks with
+ * real eigenvalues (used by the N == 2 entry point).
+ *
+ * lambdaHi/lambdaLo come from the characteristic polynomial
+ *   lambda^2 - trace*lambda + det = 0.
+ * The eigenvector for lambdaHi is v = (b, lambdaHi - a) (from the first
+ * row of (A - lambda*I)v = 0), normalized to unit length. If b == 0 the
+ * block is triangular and the eigenvectors are coordinate vectors:
+ * e1 for the larger of {a, d}, e2 for the other.
+ */
+template <uint8_t N>
+static void Solve2x2Eigen(const Matrix<N, N> &A, uint8_t lo, float &lambdaHi,
+                          float &lambdaLo, float &c, float &s) {
+  float a = A.Get(lo, lo);
+  float b = A.Get(lo, lo + 1);
+  float e = A.Get(lo + 1, lo);
+  float d = A.Get(lo + 1, lo + 1);
+
+  float trace = a + d;
+  float det = a * d - b * e;
+  float disc = trace * trace - 4.0f * det;
+  if (disc < 0.0f)
+    disc = 0.0f; // round-off clamp: real 2x2 blocks have disc >= 0
+  float sqrtDisc = sqrtf(disc);
+  lambdaHi = 0.5f * (trace + sqrtDisc);
+  lambdaLo = 0.5f * (trace - sqrtDisc);
+
+  if (b != 0.0f) {
+    float v1 = lambdaHi - a;
+    float n = sqrtf(b * b + v1 * v1);
+    c = b / n;
+    s = v1 / n;
+  } else if (a >= d) {
+    c = 1.0f; // e1 is the eigenvector of a = lambdaHi
+    s = 0.0f;
+  } else {
+    c = 0.0f; // e2 is the eigenvector of d = lambdaHi
+    s = 1.0f;
+  }
+}
+
+/**
+ * Deflate: zero subdiagonal entries i in [lo, hi) whose magnitude is at or
+ * below tolerance * (|A[i][i]| + |A[i+1][i+1]|).
+ */
+template <uint8_t N>
+static void Deflate(Matrix<N, N> &A, uint8_t lo, uint8_t hi, float tolerance) {
+  for (uint8_t i = lo; i < hi; ++i) {
+    float t = A.Get(i + 1, i);
+    float scale = fabsf(A.Get(i, i)) + fabsf(A.Get(i + 1, i + 1));
+    if (fabsf(t) <= tolerance * scale) {
+      A[i + 1][i] = 0.0f;
+      A[i][i + 1] = 0.0f;
+    }
+  }
+}
+
+// ============================================================================
+// QR::EigenQR driver (implicit Wilkinson-shifted QR, bulge chasing)
+// ============================================================================
+
+/**
+ * Tridiagonalize: Givens tridiagonalization (Golub & Van Loan 8.3.1).
+ *
+ * For column k = 0..N-3 the entries A[k+2..N-1, k] are eliminated by
+ * rotations on (i, i+1) applied BOTTOM-UP, i = N-2 down to k+1, each
+ * formed from the CURRENT (already-updated) pair (A[i][k], A[i+1][k]).
+ * Bottom-up is essential: a top-down pass zeros A[i+1][k] with a rotation
+ * that would later be undone when the next rotation (i+1, i+2) is formed
+ * from an entry below, reviving A[i][k]. Each bottom-up rotation zeros the
+ * bottom of the remaining nonzero pair and the entries below stay zero
+ * (they are not mixed again, only rows i-1/i are mixed next).
+ *
+ * Already-tridiagonalized leading columns j < k are untouched: the mixed
+ * rows are both >= k+1 > j+1, so A[i][j] and A[i+1][j] are both zero there.
+ * The rotation on (i, i+1) also keeps column k+1..k+2 structure intact and
+ * does not destroy earlier columns, so after column k is done the leading
+ * (k+1)x(k+1) block is tridiagonal forever.
+ *
+ * On return: A is symmetric tridiagonal and A_orig = U A U^T (U = product
+ * of every rotation applied, in application order, as U <- U G^T).
+ */
+template <uint8_t N>
+static void Tridiagonalize(Matrix<N, N> &A, Matrix<N, N> &U) {
+  U = Matrix<N, N>{0};
+  for (uint8_t i = 0; i < N; ++i) {
+    U[i][i] = 1.0f;
+  }
+  float c = 0.0f, s = 0.0f;
+  for (uint8_t k = 0; k + 2 < N; ++k) {
+    for (int i = (int)N - 2; i >= (int)k + 1; --i) {
+      GivensRotation(A.Get(i, k), A.Get(i + 1, k), c, s);
+      ApplyRotationBothSides(A, (uint8_t)i, c, s);
+      ApplyRotationToVectors(U, (uint8_t)i, c, s);
+    }
+  }
+}
+
+/**
+ * See QR.hpp for the full contract. Implementation sketch:
+ *
+ *   Phase 0 (N >= 3): Tridiagonalize(A, U)  // A_orig = U A U^T
+ *   V = I.
+ *   while (hi > 0):
+ *     Deflate(A, 0, hi, tol); peel exact-zero trailing subdiagonals (hi--)
+ *     lo = top of the trailing unreduced block (scan down, stop at first
+ *          exact zero subdiagonal)
+ *     if lo == hi - 1:  closed-form 2x2 eigen-solve; fold Vblock into V
+ *     else:             one implicit Wilkinson-shifted QR step:
+ *                        mu = WilkinsonShift(A[hi-1][hi-1], A[hi][hi-1], A[hi][hi])
+ *                        A[lo..hi diagonal] -= mu          // whole block!
+ *                        G1 = Givens(A[lo][lo], A[lo+1][lo])
+ *                        for i = lo..hi-1:
+ *                          (i > lo: Gi = Givens(A[i][i], A[i+1][i]))
+ *                          ApplyRotationBothSides(A, i, Gi)   // A <- Gi A Gi^T
+ *                          ApplyRotationToVectors(V, i, Gi)   // V <- V Gi^T
+ *                        A[lo..hi diagonal] += mu
+ *   eigenvalues = diag(A), sorted descending with matching V column swaps.
+ *   eigenvectors = U * V.
+ *
+ * Invariant maintained for N >= 3 (symmetric input): A is symmetric
+ * tridiagonal (up to deflated zeros and ~1e-7 float roundoff in the
+ * off-tridiagonal corners) at the top of every loop iteration, and
+ * A_orig = U A U^T = (U V) A (U V)^T throughout (V = product of every
+ * rotation applied so far, in application order, as V <- V Gi^T). At
+ * convergence A = V D V^T and therefore A_orig = (U V) D (U V)^T.
+ *
+ * Orientation note: each chase rotation Gi is the ZEROING rotation
+ * (Gi * (x, y)^T = (r, 0)^T). The step A <- Gi A Gi^T equals R Q with
+ * R = Gi A upper-triangular (on the block) and Q = Gi^T -- i.e. it IS the
+ * standard QR update Q(A - mu I)Q^T with Q the orthogonal QR factor. The
+ * eigenvector accumulator therefore collects the Q factors: V <- V Gi^T.
+ */
+template <uint8_t N>
+void EigenQR(Matrix<N, N> &matrixToDecompose, Matrix<N, N> &eigenVectors,
+             Matrix<N, 1> &eigenValues, uint32_t maxIterations, float tolerance) {
+  static_assert(N >= 2, "QR::EigenQR requires N >= 2 (N = 1 is trivial)");
+
+  Matrix<N, N> A = matrixToDecompose; // input is not modified
+  Matrix<N, N> V{0};
+  // NB: Matrix::Identity() is a static factory that returns by value; a
+  // bare call would be a no-op. Set the diagonal explicitly.
+  for (uint8_t i = 0; i < N; ++i) {
+    V[i][i] = 1.0f;
+  }
+
+  // ------------------------------------------------------------------
+  // N == 2: closed-form solution (works for nonsymmetric input too)
+  // ------------------------------------------------------------------
+  if (N == 2) {
+    float l1 = 0.0f, l2 = 0.0f, c = 0.0f, s = 0.0f;
+    Solve2x2Eigen(A, 0, l1, l2, c, s);
+    // V = I * Vblock = [[c, -s], [s, c]]
+    V[0][0] = c;
+    V[0][1] = -s;
+    V[1][0] = s;
+    V[1][1] = c;
+    eigenValues[0][0] = l1;
+    eigenValues[1][0] = l2;
+    for (uint8_t r = 0; r < N; ++r)
+      for (uint8_t col = 0; col < N; ++col)
+        eigenVectors[r][col] = V.Get(r, col);
+    return;
+  }
+
+  // ------------------------------------------------------------------
+  // N >= 3: implicit shifted QR iteration (symmetric input required)
+  // ------------------------------------------------------------------
+  // Phase 0: general symmetric -> symmetric tridiagonal. The implicit
+  // QR bulge chase only preserves a tridiagonal structure, so the input
+  // must be reduced first: A_orig = U A U^T with A tridiagonal.
+  Matrix<N, N> U{};
+  Tridiagonalize(A, U);
+
+  uint32_t iter = 0;
+  uint8_t hi = N - 1;
+  while (hi > 0) {
+    Deflate(A, 0, hi, tolerance);
+    // Peel trailing rows whose subdiagonal is exactly zero (deflated or
+    // already solved). Must be re-done every iteration: a peel is only
+    // meaningful once the subdiagonal beneath it has converged.
+    while (hi > 0 && A.Get(hi, hi - 1) == 0.0f) {
+      --hi;
+    }
+    if (hi == 0) {
+      break; // fully diagonal (within tolerance)
+    }
+
+    // Find the top of the trailing unreduced block: scan down from hi-1
+    // and stop at the first exact zero subdiagonal. A[hi][hi-1] != 0 here
+    // (just peeled), so lo < hi.
+    uint8_t lo = hi;
+    for (int i = (int)hi - 1; i >= 0; --i) {
+      if (A.Get(i + 1, i) == 0.0f) {
+        break;
+      }
+      lo = (uint8_t)i;
+    }
+
+    if (lo + 1 == hi) {
+      // Trailing unreduced block is 2x2: solve in closed form.
+      float l1 = 0.0f, l2 = 0.0f, c = 0.0f, s = 0.0f;
+      Solve2x2Eigen(A, lo, l1, l2, c, s);
+      A[lo][lo] = l1;
+      A[lo + 1][lo + 1] = l2;
+      A[lo][lo + 1] = 0.0f;
+      A[lo + 1][lo] = 0.0f;
+      // Fold Vblock = [[c, -s], [s, c]] into V: V <- V * Vblock on
+      // columns (lo, lo+1). NOTE the sign convention differs from
+      // ApplyRotationToVectors (which applies [[c, s], [-s, c]]):
+      // here column 0 of Vblock is (c, s)^T, column 1 is (-s, c)^T.
+      for (uint8_t r = 0; r < N; ++r) {
+        float x = V.Get(r, lo);
+        float y = V.Get(r, lo + 1);
+        V[r][lo] = c * x + s * y;
+        V[r][lo + 1] = -s * x + c * y;
+      }
+      if (lo == 0) {
+        break; // block reached the top: matrix is fully solved
+      }
+      hi = (uint8_t)(lo - 1);
+      continue;
+    }
+
+    // One implicit Wilkinson-shifted QR step on block [lo, hi].
+    float mu = WilkinsonShift(A.Get(hi - 1, hi - 1), A.Get(hi, hi - 1),
+                              A.Get(hi, hi));
+
+    // The shift applies to the ENTIRE active block: bulge chasing
+    // triangularizes (A - mu*I), and the first Givens rotation is formed
+    // from (A[lo][lo] - mu, A[lo+1][lo]).
+    for (uint8_t i = lo; i <= hi; ++i) {
+      A[i][i] -= mu;
+    }
+
+    float c = 0.0f, s = 0.0f;
+    for (uint8_t i = lo; i < hi; ++i) {
+      if (i == lo) {
+        GivensRotation(A.Get(lo, lo), A.Get(lo + 1, lo), c, s);
+      } else {
+        GivensRotation(A.Get(i, i), A.Get(i + 1, i), c, s);
+      }
+      ApplyRotationBothSides(A, i, c, s);
+      ApplyRotationToVectors(V, i, c, s);
+    }
+
+    for (uint8_t i = lo; i <= hi; ++i) {
+      A[i][i] += mu;
+    }
+
+    if (++iter >= maxIterations) {
+      // Best-effort: fall through with the partially diagonalized A.
+      break;
+    }
+  }
+
+  // ------------------------------------------------------------------
+  // Collect eigenvalues and sort DESCENDING (swap eigenvectors to match)
+  // ------------------------------------------------------------------
+  for (uint8_t i = 0; i < N; ++i) {
+    eigenValues[i][0] = A.Get(i, i);
+  }
+  for (uint8_t i = 0; i < N - 1; ++i) {
+    uint8_t k = i;
+    for (uint8_t j = i + 1; j < N; ++j) {
+      if (eigenValues.Get(j, 0) > eigenValues.Get(k, 0)) {
+        k = j;
+      }
+    }
+    if (k != i) {
+      float t = eigenValues[i][0];
+      eigenValues[i][0] = eigenValues[k][0];
+      eigenValues[k][0] = t;
+      for (uint8_t r = 0; r < N; ++r) {
+        float x = V.Get(r, i);
+        V[r][i] = V.Get(r, k);
+        V[r][k] = x;
+      }
+    }
+  }
+
+  // True eigenvectors of the original matrix: U * V. Reuse the A buffer
+  // (its diagonal has already been collected into eigenValues).
+  U.Mult(V, A);
+
+  for (uint8_t r = 0; r < N; ++r) {
+    for (uint8_t col = 0; col < N; ++col) {
+      eigenVectors[r][col] = A.Get(r, col);
+    }
+  }
+}
+
+} // namespace QR
+
+#endif // QR_H_

src/QR.hpp

@@ -0,0 +1,224 @@
+#pragma once
+#include "Matrix.hpp"
+
+/**
+ * @brief Library that uses Matrix.hpp and computes the eigenvalues and
+ * eigenvectors of a square matrix with the implicit shifted QR iteration
+ * (Wilkinson shift, Givens bulge chasing).
+ *
+ * @note Fully templated: QR::EigenQR works for ANY Matrix<N,N> with N in
+ *       2..255 (the uint8_t range of Matrix). There is no 5x5 limit.
+ *
+ * @note N >= 3: the input matrix MUST be symmetric (A[i][j] == A[j][i]).
+ *       The implicit QR bulge chase maintains a symmetric tridiagonal
+ *       structure, which only exists for symmetric input. N = 2 handles
+ *       a general (nonsymmetric) 2x2 via the closed-form solution, so
+ *       nonsymmetric 2x2 inputs also work.
+ *
+ * @note The input matrix is NOT modified (the iteration runs on a local
+ *       copy), mirroring the SVD::SVD convention.
+ *
+ * @note EMBEDDED CONSTRAINT -- no heap. All working storage is stack
+ *       allocated as templated Matrix<N,N> buffers. Peak stack usage per
+ *       call is 3 * N^2 floats (A working copy + U and V accumulators) =
+ *       12 * N^2 bytes:
+ *         N = 5   -> ~0.3 KB
+ *         N = 10  -> ~1.2 KB
+ *         N = 20  -> ~4.8 KB
+ *         N = 50  -> ~30 KB
+ *         N = 100 -> ~120 KB
+ *         N = 255 -> ~783 KB
+ *       Instantiate only the sizes that fit your call-stack budget.
+ *
+ * @note Conventions:
+ *       - Eigenvalues come out sorted DESCENDING (largest first); the
+ *         eigenvector columns are swapped to match.
+ *       - Eigenvector signs are arbitrary (v and -v are both valid);
+ *         tests must be sign-invariant.
+ *       - Wilkinson shift: the eigenvalue of the trailing 2x2 block
+ *         closest to the bottom-right corner (Trefethen & Bau 13.4.1).
+ *
+ * @note Algorithm (Trefethen & Bau 13.4, Golub & Van Loan 8.4.3):
+ *       Phase 0 (N >= 3): Givens tridiagonalization. A general symmetric
+ *         matrix is NOT suitable for implicit QR (the bulge chase only
+ *         preserves the tridiagonal structure), so first reduce A with
+ *         adjacent Givens similarities A <- G A G^T (rotations applied
+ *         BOTTOM-UP, i = N-2 down to k+1, per column k), accumulating
+ *         U <- U G^T, until A is symmetric tridiagonal and
+ *         A_orig = U A U^T. (N = 2 needs no reduction.)
+ *       Phase 1: iterate until A is diagonal:
+ *         1. Deflate: zero out subdiagonal entries at/under the tolerance
+ *            (scaled by the adjacent diagonal magnitudes).
+ *         2. Scan for the trailing unreduced block [lo, hi].
+ *            - block of size 1: A[hi][hi] is a converged eigenvalue, done.
+ *            - block of size 2: solve the 2x2 eigenproblem in closed form
+ *              and fold its eigenvector matrix into V.
+ *            - block larger: one implicit Wilkinson-shifted QR step
+ *              (bulge chasing with Givens rotations; the shift is applied
+ *              to the ENTIRE active block [lo, hi], not just the trailing
+ *              2x2 -- the first Givens rotation must be formed from
+ *              (A[lo][lo] - mu, A[lo+1][lo])). Every rotation is folded
+ *              into V.
+ *       Phase 2: eigenvalues = diag(A), sorted DESCENDING (eigenvector
+ *         columns swapped to match), and the true eigenvectors of the
+ *         ORIGINAL matrix are U * V.
+ *
+ * @note If maxIterations is exhausted before convergence the best-effort
+ *       (partially diagonalized) values on the diagonal are returned.
+ */
+namespace QR {
+
+/**
+ * @brief Compute the eigenvalues and eigenvectors of a square matrix
+ *
+ * @param matrixToDecompose The matrix to take eigenvalues of (not
+ *        modified). MUST be symmetric for N >= 3.
+ * @param eigenVectors a buffer that will contain the eigenvectors in its
+ *        COLUMNS, sorted by descending eigenvalue (column i is the
+ *        eigenvector for eigenValues[i]).
+ * @param eigenValues a buffer that will contain the eigenvalues sorted
+ *        DESCENDING (largest first).
+ * @param maxIterations the number of QR steps to perform before giving up
+ *        on reaching the given tolerance
+ * @param tolerance the level of accuracy to obtain before stopping; a
+ *        subdiagonal entry is deflated when |A[i+1][i]| <= tolerance *
+ *        (|A[i][i]| + |A[i+1][i+1]|). For float32 arithmetic, values
+ *        around 1e-6 are a sensible choice (single-precision epsilon is
+ *        ~1.2e-7).
+ */
+template <uint8_t N>
+void EigenQR(Matrix<N, N> &matrixToDecompose, Matrix<N, N> &eigenVectors,
+             Matrix<N, 1> &eigenValues, uint32_t maxIterations, float tolerance);
+
+/**
+ * @brief Compute a Givens rotation that zeros the bottom entry of (a, b)
+ *
+ * Given the column vector (a, b), produces (c, s) defining the 2x2
+ * rotation
+ *        R = [ c  s ]
+ *            [ -s c ]
+ * such that R * (a, b)^T = (r, 0)^T with r = +hypot(a, b) >= 0, i.e.
+ * c = a / r, s = b / r.
+ *
+ * If (a, b) == (0, 0) the identity rotation (c = 1, s = 0) is returned.
+ */
+static void GivensRotation(float a, float b, float &c, float &s);
+
+/**
+ * @brief Apply the similarity transform A <- G A G^T on rows/cols (i, i+1)
+ *
+ * G = [ c  s ] on the (i, i+1) block, identity elsewhere, where G is the
+ *     [ -s c ]
+ * ZEROING rotation (G * (x, y)^T = (r, 0)^T) -- the orientation used by
+ * the implicit QR chase: A = Q R with Q = G^T gives the next iterate
+ * R Q = G A G^T. With (c, s) = GivensRotation(A[i][i], A[i+1][i]) the
+ * (i+1, i) entry is zeroed by the left multiplication and the bulge is
+ * chased along the superdiagonal by the right one. The matrix must be
+ * symmetric on entry (guaranteed by construction in the QR iteration:
+ * symmetric input stays symmetric under similarity by an orthogonal
+ * matrix). Updates the full matrix, not just the tridiagonal structure.
+ */
+template <uint8_t N>
+static void ApplyRotationBothSides(Matrix<N, N> &A, uint8_t i, float c,
+                                   float s);
+
+/**
+ * @brief Accumulate eigenvectors: V <- V G^T on columns (i, i+1)
+ *
+ * G^T = [ c  -s ] on columns (i, i+1), identity elsewhere, where G =
+ *       [ s   c ]
+ * [ c, s ] / [ -s, c ] is the zeroing rotation paired with
+ * ApplyRotationBothSides. Applied to all rows:
+ *   V[r][i]     -> c V[r][i]   + s V[r][i+1]
+ *   V[r][i+1]   -> -s V[r][i]  + c V[r][i+1]
+ *
+ * Every QR step's rotation is folded into V this way so that, together
+ * with A <- G A G^T, the invariant A_orig = V A V^T is preserved at every
+ * step (each step is A <- R Q with Q = G^T the orthogonal factor, and
+ * the orthogonal factors multiply as G1^T G2^T ... in application order).
+ * At convergence A_orig = V D V^T and the columns of V are the
+ * eigenvectors.
+ */
+template <uint8_t N>
+static void ApplyRotationToVectors(Matrix<N, N> &V, uint8_t i, float c,
+                                   float s);
+
+/**
+ * @brief Wilkinson shift for a symmetric tridiagonal
+ *
+ * Given the trailing 2x2 block
+ *        [ a  b ]
+ *        [ b  d ]
+ * returns the eigenvalue of that block that is closest to d. This is the
+ * empirically best shift for the QR iteration (Trefethen & Bau 13.4.1).
+ *
+ * mu = (a+d)/2 - sign(a-d) * sqrt(((a-d)/2)^2 + b^2)
+ * (with sign(0) taken as +1).
+ */
+static float WilkinsonShift(float a, float b, float d);
+
+/**
+ * @brief Solve the 2x2 eigenproblem of block rows/cols (lo, lo+1)
+ *
+ * Solves the (possibly nonsymmetric) 2x2 block
+ *        [ A[lo][lo]     A[lo][lo+1]   ]
+ *        [ A[lo+1][lo]   A[lo+1][lo+1] ]
+ * in closed form (characteristic polynomial + eigenvector back-substitution).
+ *
+ * @param A the matrix containing the block (not modified)
+ * @param lo the row/col index of the top-left corner of the block
+ * @param lambdaHi (out) the LARGER eigenvalue
+ * @param lambdaLo (out) the smaller eigenvalue
+ * @param c (out), s (out) eigenvector pair as an orthogonal matrix
+ *        Vblock = [ c  -s ] whose columns are the eigenvectors: column 0
+ *                 [ s   c ]
+ *        (c, s) is the unit eigenvector for lambdaHi, column 1 (-s, c) is
+ *        the unit eigenvector for lambdaLo.
+ *
+ * Note: the caller applies Vblock to its eigenvector accumulator with
+ * V <- V * Vblock (i.e. V[r][lo]   = c*x + s*y,
+ *                              V[r][lo+1] = -s*x + c*y). Vblock has the
+ * SAME [ c -s; s c ] form as the G^T factor used by
+ * ApplyRotationToVectors, so both folding operations follow one uniform
+ * convention.
+ */
+template <uint8_t N>
+static void Solve2x2Eigen(const Matrix<N, N> &A, uint8_t lo, float &lambdaHi,
+                          float &lambdaLo, float &c, float &s);
+
+/**
+ * @brief Deflate (zero out) subdiagonal entries that are at/under tolerance
+ *
+ * For each i in [lo, hi): if |A[i+1][i]| <= tolerance *
+ * (|A[i][i]| + |A[i+1][i+1]|), sets A[i+1][i] = A[i][i+1] = 0, splitting
+ * the matrix into smaller independent blocks.
+ */
+template <uint8_t N>
+static void Deflate(Matrix<N, N> &A, uint8_t lo, uint8_t hi, float tolerance);
+
+/**
+ * @brief Reduce a symmetric matrix to symmetric tridiagonal form
+ *
+ * Chases each column's entries below the subdiagonal to zero with
+ * adjacent Givens similarities (Golub & Van Loan 8.3.1, Givens variant):
+ * for column k = 0..N-3, rotations on (N-2, N-1), (N-3, N-2), ...
+ * (k+1, k+2) -- BOTTOM-UP, each formed from the current (A[i][k],
+ * A[i+1][k]) -- zero A[k+2..N-1, k] one by one. A top-down pass would not
+ * work: the rotation that zeros A[i+1][k] would be undone by the later
+ * rotation on (i+1, i+2) forming a new nonzero at A[i][k]. Each rotation
+ * is applied to A as a similarity (A <- G A G^T) and accumulated into U
+ * (U <- U G^T), so on return:
+ *   - A is symmetric tridiagonal (off-tridiagonal entries EXACTLY zero),
+ *   - A_orig = U A U^T (i.e. U^T A_orig U = A).
+ *
+ * U is initialized to the identity internally (its input contents are
+ * ignored).
+ */
+template <uint8_t N>
+static void Tridiagonalize(Matrix<N, N> &A, Matrix<N, N> &U);
+
+} // namespace QR
+
+#ifndef QR_H_
+#include "QR.cpp"
+#endif

unit-tests/CMakeLists.txt

@@ -13,6 +13,7 @@ add_executable(matrix-tests matrix-tests.cpp)
 target_link_libraries(matrix-tests 
     PRIVATE 
     matrix
+    qr
     Catch2::Catch2WithMain
 )
 
@@ -52,4 +53,14 @@ target_link_libraries(svd-integration-test
     matrix
     svd
     Catch2::Catch2WithMain
+)
+
+# QR building block tests
+add_executable(qr-build-blocks-tests qr-build-blocks-tests.cpp)
+
+target_link_libraries(qr-build-blocks-tests
+    PRIVATE
+    matrix
+    qr
+    Catch2::Catch2WithMain
 )

unit-tests/matrix-tests.cpp

@@ -4,6 +4,7 @@
 
 // include the module you're going to test next
 #include "Matrix.hpp"
+#include "QR.hpp"
 #include "SVD.hpp"
 
 // any other libraries
@@ -601,8 +602,78 @@ TEST_CASE("QR Decompositions", "Matrix") {
   }
 }
 
+// ============================================================================
+// Eigen QR Helpers (scipy references; eigenvector checks are sign-invariant)
+// ============================================================================
+
+/**
+ * @brief Normalized eigenpair residual ||A v - lambda v|| / (||A||_F + |lambda|)
+ */
+template <uint8_t N>
+static float eigenResidual(const Matrix<N, N> &A, float lambda,
+                           const Matrix<N, 1> &v) {
+  Matrix<N, 1> Av{};
+  A.Mult(v, Av);
+  float sum = 0.0f;
+  float frob = 0.0f;
+  for (uint8_t i = 0; i < N; i++) {
+    float d = Av.Get(i, 0) - lambda * v.Get(i, 0);
+    sum += d * d;
+    for (uint8_t j = 0; j < N; j++) {
+      float a = A.Get(i, j);
+      frob += a * a;
+    }
+  }
+  float scale = sqrtf(frob) + fabsf(lambda);
+  return sqrtf(sum) / scale;
+}
+
+/**
+ * @brief Column of the eigenvector matrix; used for the residual check.
+ */
+template <uint8_t N>
+static Matrix<N, 1> eigenColumn(const Matrix<N, N> &V, uint8_t col) {
+  Matrix<N, 1> v{};
+  for (uint8_t i = 0; i < N; i++) {
+    v[i][0] = V.Get(i, col);
+  }
+  return v;
+}
+
+/**
+ * @brief Check V^T V ~ I (eigenvectors orthonormal).
+ */
+template <uint8_t N>
+static bool isOrthogonal(const Matrix<N, N> &V, float tol = 1e-4f) {
+  Matrix<N, N> Vt = V.Transpose();
+  Matrix<N, N> VtV{};
+  Vt.Mult(V, VtV);
+  for (uint8_t i = 0; i < N; i++) {
+    for (uint8_t j = 0; j < N; j++) {
+      float expected = (i == j) ? 1.0f : 0.0f;
+      if (fabsf(VtV.Get(i, j) - expected) > tol) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+/**
+ * @brief Sign-invariant component check: |actual| within max(1e-4, 1e-3*|ref|)
+ *        of ref (ref is the ABSOLUTE value from the scipy reference).
+ */
+static bool componentMatches(float actual, float refAbs) {
+  float a = fabsf(actual);
+  float tol = 1e-4f;
+  if (refAbs * 1e-3f > tol) {
+    tol = refAbs * 1e-3f;
+  }
+  return fabsf(a - refAbs) <= tol;
+}
+
 TEST_CASE("Eigenvalues and Vectors", "Matrix") {
-  SECTION("2x2 Eigen") {
+  SECTION("2x2 Eigen (nonsymmetric, closed form)") {
     Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
     Matrix<2, 2> vectors{};
     Matrix<2, 1> values{};
@@ -615,28 +686,206 @@ TEST_CASE("Eigenvalues and Vectors", "Matrix") {
     REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(-0.372281f, 1e-4f));
   }
 
-  SECTION("3x3 Rank Defficient Eigen") {
+  // Reference values: numpy.linalg.eigh on float32 matrices.
-    SKIP("Skipping this because QR decomposition isn't ready for it");
+  // Eigenvector component references are ABSOLUTE values (signs arbitrary).
-    // this symmetrix tridiagonal matrix is well behaved for testing
-    Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
 
+  SECTION("3x3 Symmetric Eigen") {
+    Matrix<3, 3> A{1, 2, 3, 2, 5, 8, 3, 8, 9};
     Matrix<3, 3> vectors{};
     Matrix<3, 1> values{};
-    A.EigenQR(vectors, values, 1000000, 1e-8f);
+    A.EigenQR(vectors, values, 10000, 1e-6f);
 
-    std::string strBuf1 = "";
+    // eigenvalues (descending)
-    vectors.ToString(strBuf1);
+    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(16.102417f, 1e-4f));
-    std::cout << "Vectors:\n" << strBuf1 << std::endl;
+    REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(0.191920f, 1e-4f));
-    strBuf1 = "";
+    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(-1.2943381f, 1e-4f));
-    values.ToString(strBuf1);
-    std::cout << "Values:\n" << strBuf1 << std::endl;
 
-    REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.23197f, 1e-4f));
+    // eigenvector |components| (sign-invariant)
-    REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.525322f, 1e-4f));
+    REQUIRE(componentMatches(vectors[0][0], 0.231657207f));
-    REQUIRE_THAT(vectors[2][0], Catch::Matchers::WithinRel(0.81867f, 1e-4f));
+    REQUIRE(componentMatches(vectors[1][0], 0.59582746f));
-    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(-1.11684f, 1e-4f));
+    REQUIRE(componentMatches(vectors[2][0], 0.768976331f));
+    REQUIRE(componentMatches(vectors[0][1], 0.956842422f));
+    REQUIRE(componentMatches(vectors[1][1], 0.282139271f));
+    REQUIRE(componentMatches(vectors[2][1], 0.0696421042f));
+    REQUIRE(componentMatches(vectors[0][2], 0.175463736f));
+    REQUIRE(componentMatches(vectors[1][2], 0.75192225f));
+    REQUIRE(componentMatches(vectors[2][2], 0.635472536f));
+
+    // eigenvectors orthonormal; eigenpair residuals small
+    REQUIRE(isOrthogonal(vectors));
+    for (uint8_t col = 0; col < 3; col++) {
+      REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) <
+              1e-4f);
+    }
+  }
+
+  SECTION("3x3 Rank Deficient Eigen") {
+    // A = v v^T with v = [1, 2, 3]: eigenvalues {14, 0, 0}
+    Matrix<3, 3> A{1, 2, 3, 2, 4, 6, 3, 6, 9};
+    Matrix<3, 3> vectors{};
+    Matrix<3, 1> values{};
+    A.EigenQR(vectors, values, 10000, 1e-6f);
+
+    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(14.0f, 1e-4f));
     REQUIRE_THAT(values[1][0], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
-    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(16.1168f, 1e-4f));
+    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+
+    // dominant eigenvector is v/|v| (sign-invariant); the two null-space
+    // eigenvectors may be ANY orthonormal basis of the null plane, so only
+    // orthogonality + residuals are checked for the full matrix.
+    REQUIRE(componentMatches(vectors[0][0], 0.267261237f));
+    REQUIRE(componentMatches(vectors[1][0], 0.534522474f));
+    REQUIRE(componentMatches(vectors[2][0], 0.801783741f));
+    REQUIRE(isOrthogonal(vectors));
+    for (uint8_t col = 0; col < 3; col++) {
+      REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) <
+              1e-4f);
+    }
+  }
+
+  SECTION("4x4 Symmetric Eigen") {
+    Matrix<4, 4> A{2, 1, 0, 1, 1, 3, 1, 0, 0, 1, 4, 1, 1, 0, 1, 5};
+    Matrix<4, 4> vectors{};
+    Matrix<4, 1> values{};
+    A.EigenQR(vectors, values, 10000, 1e-6f);
+
+    // eigenvalues are exactly {6, 4, 3, 1}
+    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(6.0f, 1e-4f));
+    REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(4.0f, 1e-4f));
+    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(3.0f, 1e-4f));
+    REQUIRE_THAT(values[3][0], Catch::Matchers::WithinRel(1.0f, 1e-4f));
+
+    // eigenvector |components| (sign-invariant)
+    REQUIRE(componentMatches(vectors[0][0], 0.258198887f));
+    REQUIRE(componentMatches(vectors[1][0], 0.258198887f));
+    REQUIRE(componentMatches(vectors[2][0], 0.516397774f));
+    REQUIRE(componentMatches(vectors[3][0], 0.774596691f));
+    REQUIRE(componentMatches(vectors[0][1], 0.0f));
+    REQUIRE(componentMatches(vectors[1][1], 0.577350259f));
+    REQUIRE(componentMatches(vectors[2][1], 0.577350259f));
+    REQUIRE(componentMatches(vectors[3][1], 0.577350259f));
+    REQUIRE(componentMatches(vectors[0][2], 0.577350259f));
+    REQUIRE(componentMatches(vectors[1][2], 0.577350259f));
+    REQUIRE(componentMatches(vectors[2][2], 0.577350259f));
+    REQUIRE(componentMatches(vectors[3][2], 0.0f));
+    REQUIRE(componentMatches(vectors[0][3], 0.774596691f));
+    REQUIRE(componentMatches(vectors[1][3], 0.516397774f));
+    REQUIRE(componentMatches(vectors[2][3], 0.258198887f));
+    REQUIRE(componentMatches(vectors[3][3], 0.258198887f));
+
+    REQUIRE(isOrthogonal(vectors));
+    for (uint8_t col = 0; col < 4; col++) {
+      REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) <
+              1e-4f);
+    }
+  }
+
+  SECTION("5x5 Symmetric Eigen") {
+    Matrix<5, 5> A{3, 1, 0, 0, 1, 1, 4, 1, 0, 0, 0, 1, 5, 1, 0, 0, 0, 1, 6, 1,
+                   1, 0, 0, 1, 7};
+    Matrix<5, 5> vectors{};
+    Matrix<5, 1> values{};
+    A.EigenQR(vectors, values, 10000, 1e-6f);
+
+    // eigenvalues (descending)
+    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(7.90154457f, 1e-4f));
+    REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(6.20044184f, 1e-4f));
+    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(5.14503145f, 1e-4f));
+    REQUIRE_THAT(values[3][0], Catch::Matchers::WithinRel(3.61823463f, 1e-4f));
+    REQUIRE_THAT(values[4][0], Catch::Matchers::WithinRel(2.13474774f, 1e-4f));
+
+    // eigenvector |components| (sign-invariant)
+    REQUIRE(componentMatches(vectors[0][0], 0.182430908f));
+    REQUIRE(componentMatches(vectors[1][0], 0.102749094f));
+    REQUIRE(componentMatches(vectors[2][0], 0.21844925f));
+    REQUIRE(componentMatches(vectors[3][0], 0.531091094f));
+    REQUIRE(componentMatches(vectors[4][0], 0.791444063f));
+    REQUIRE(componentMatches(vectors[0][1], 0.0877681747f));
+    REQUIRE(componentMatches(vectors[1][1], 0.245861098f));
+    REQUIRE(componentMatches(vectors[2][1], 0.628771126f));
+    REQUIRE(componentMatches(vectors[3][1], 0.508941948f));
+    REQUIRE(componentMatches(vectors[4][1], 0.526758015f));
+    REQUIRE(componentMatches(vectors[0][2], 0.349721253f));
+    REQUIRE(componentMatches(vectors[1][2], 0.628706098f));
+    REQUIRE(componentMatches(vectors[2][2], 0.370167077f));
+    REQUIRE(componentMatches(vectors[3][2], 0.575020194f));
+    REQUIRE(componentMatches(vectors[4][2], 0.121457018f));
+    REQUIRE(componentMatches(vectors[0][3], 0.429638386f));
+    REQUIRE(componentMatches(vectors[1][3], 0.498553723f));
+    REQUIRE(componentMatches(vectors[2][3], 0.619968951f));
+    REQUIRE(componentMatches(vectors[3][3], 0.35809797f));
+    REQUIRE(componentMatches(vectors[4][3], 0.232936427f));
+    REQUIRE(componentMatches(vectors[0][4], 0.807540476f));
+    REQUIRE(componentMatches(vectors[1][4], 0.534011006f));
+    REQUIRE(componentMatches(vectors[2][4], 0.188524753f));
+    REQUIRE(componentMatches(vectors[3][4], 0.00615991838f));
+    REQUIRE(componentMatches(vectors[4][4], 0.164715111f));
+
+    REQUIRE(isOrthogonal(vectors));
+    for (uint8_t col = 0; col < 5; col++) {
+      REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) <
+              1e-4f);
+    }
+  }
+
+  SECTION("6x6 Symmetric Eigen") {
+    Matrix<6, 6> A{4, 1, 0, 0, 0, 1, 1, 5, 1, 0, 0, 0, 0, 1, 6, 1, 0, 0, 0, 0,
+                   1, 7, 1, 0, 0, 0, 0, 1, 8, 1, 1, 0, 0, 0, 1, 3};
+    Matrix<6, 6> vectors{};
+    Matrix<6, 1> values{};
+    A.EigenQR(vectors, values, 10000, 1e-6f);
+
+    // eigenvalues (descending)
+    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(8.86080551f, 1e-4f));
+    REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(7.25410175f, 1e-4f));
+    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(6.11490774f, 1e-4f));
+    REQUIRE_THAT(values[3][0], Catch::Matchers::WithinRel(4.88509226f, 1e-4f));
+    REQUIRE_THAT(values[4][0], Catch::Matchers::WithinRel(3.74589825f, 1e-4f));
+    REQUIRE_THAT(values[5][0], Catch::Matchers::WithinRel(2.13919425f, 1e-4f));
+
+    // eigenvector |components| (sign-invariant)
+    REQUIRE(componentMatches(vectors[0][0], 0.0430923924f));
+    REQUIRE(componentMatches(vectors[1][0], 0.0662503168f));
+    REQUIRE(componentMatches(vectors[2][0], 0.212687209f));
+    REQUIRE(componentMatches(vectors[3][0], 0.542206466f));
+    REQUIRE(componentMatches(vectors[4][0], 0.7962538f));
+    REQUIRE(componentMatches(vectors[5][0], 0.143213451f));
+    REQUIRE(componentMatches(vectors[0][1], 0.0623276457f));
+    REQUIRE(componentMatches(vectors[1][1], 0.307613879f));
+    REQUIRE(componentMatches(vectors[2][1], 0.631065309f));
+    REQUIRE(componentMatches(vectors[3][1], 0.483806193f));
+    REQUIRE(componentMatches(vectors[4][1], 0.508129358f));
+    REQUIRE(componentMatches(vectors[5][1], 0.104793385f));
+    REQUIRE(componentMatches(vectors[0][2], 0.374228716f));
+    REQUIRE(componentMatches(vectors[1][2], 0.605694294f));
+    REQUIRE(componentMatches(vectors[2][2], 0.301064402f));
+    REQUIRE(componentMatches(vectors[3][2], 0.571099699f));
+    REQUIRE(componentMatches(vectors[4][2], 0.204411641f));
+    REQUIRE(componentMatches(vectors[5][2], 0.185764849f));
+    REQUIRE(componentMatches(vectors[0][3], 0.571099699f));
+    REQUIRE(componentMatches(vectors[1][3], 0.301064402f));
+    REQUIRE(componentMatches(vectors[2][3], 0.605694294f));
+    REQUIRE(componentMatches(vectors[3][3], 0.374228716f));
+    REQUIRE(componentMatches(vectors[4][3], 0.185764849f));
+    REQUIRE(componentMatches(vectors[5][3], 0.204411641f));
+    REQUIRE(componentMatches(vectors[0][4], 0.483806193f));
+    REQUIRE(componentMatches(vectors[1][4], 0.631065309f));
+    REQUIRE(componentMatches(vectors[2][4], 0.307613879f));
+    REQUIRE(componentMatches(vectors[3][4], 0.0623276457f));
+    REQUIRE(componentMatches(vectors[4][4], 0.104793385f));
+    REQUIRE(componentMatches(vectors[5][4], 0.508129358f));
+    REQUIRE(componentMatches(vectors[0][5], 0.542206466f));
+    REQUIRE(componentMatches(vectors[1][5], 0.212687209f));
+    REQUIRE(componentMatches(vectors[2][5], 0.0662503168f));
+    REQUIRE(componentMatches(vectors[3][5], 0.0430923924f));
+    REQUIRE(componentMatches(vectors[4][5], 0.143213451f));
+    REQUIRE(componentMatches(vectors[5][5], 0.7962538f));
+
+    REQUIRE(isOrthogonal(vectors));
+    for (uint8_t col = 0; col < 6; col++) {
+      REQUIRE(eigenResidual(A, values[col][0], eigenColumn(vectors, col)) <
+              1e-4f);
+    }
   }
 }
 

unit-tests/qr-build-blocks-tests.cpp

@@ -0,0 +1,581 @@
+// include the unit test framework first
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+// include the module you're going to test next
+#include "Matrix.hpp"
+#include "QR.hpp"
+
+// any other libraries
+#include <array>
+#include <cmath>
+#include <iostream>
+
+// ============================================================================
+// Helpers
+// ============================================================================
+
+/**
+ * @brief Frobenius norm of an N x N matrix.
+ */
+template <uint8_t N>
+static float frob(const Matrix<N, N> &M) {
+  float sum = 0.0f;
+  for (uint8_t i = 0; i < N; i++)
+    for (uint8_t j = 0; j < N; j++) {
+      float v = M.Get(i, j);
+      sum += v * v;
+    }
+  return sqrtf(sum);
+}
+
+/**
+ * @brief Check M is orthogonal (M^T M ~ I).
+ */
+template <uint8_t N>
+static bool isOrthogonal(const Matrix<N, N> &M, float tol = 1e-5f) {
+  Matrix<N, N> Mt = M.Transpose();
+  Matrix<N, N> MtM{};
+  Mt.Mult(M, MtM);
+  for (uint8_t i = 0; i < N; i++)
+    for (uint8_t j = 0; j < N; j++) {
+      float expected = (i == j) ? 1.0f : 0.0f;
+      if (fabsf(MtM.Get(i, j) - expected) > tol)
+        return false;
+    }
+  return true;
+}
+
+/**
+ * @brief 3x3 trace.
+ */
+static float trace3(const Matrix<3, 3> &A) {
+  return A.Get(0, 0) + A.Get(1, 1) + A.Get(2, 2);
+}
+
+/**
+ * @brief 3x3 sum of principal 2x2 minors (2nd elementary invariant).
+ */
+static float e2_3x3(const Matrix<3, 3> &A) {
+  return A.Get(0, 0) * A.Get(1, 1) - A.Get(0, 1) * A.Get(0, 1) +
+         A.Get(0, 0) * A.Get(2, 2) - A.Get(0, 2) * A.Get(0, 2) +
+         A.Get(1, 1) * A.Get(2, 2) - A.Get(1, 2) * A.Get(1, 2);
+}
+
+/**
+ * @brief 3x3 determinant.
+ */
+static float det3(const Matrix<3, 3> &A) {
+  return A.Get(0, 0) *
+             (A.Get(1, 1) * A.Get(2, 2) - A.Get(1, 2) * A.Get(2, 1)) -
+         A.Get(0, 1) *
+             (A.Get(1, 0) * A.Get(2, 2) - A.Get(1, 2) * A.Get(2, 0)) +
+         A.Get(0, 2) *
+             (A.Get(1, 0) * A.Get(2, 1) - A.Get(1, 1) * A.Get(2, 0));
+}
+
+/**
+ * @brief Sign-invariant comparison of |actual| against refAbs.
+ */
+static bool matchesAbs(float actual, float refAbs, float relTol = 1e-5f,
+                       float absTol = 1e-6f) {
+  float a = fabsf(actual);
+  if (refAbs < 1e-3f)
+    return a < absTol + relTol;
+  return fabsf(a - refAbs) <= relTol * refAbs;
+}
+
+// ============================================================================
+// TEST 1: GivensRotation
+// ============================================================================
+TEST_CASE("QR Building Block: GivensRotation", "[Matrix][QR]") {
+  // R = [[c, s], [-s, c]] must satisfy R * (a, b)^T = (r, 0)^T.
+
+  {
+    // Reference: hypot(2, 1) = sqrt(5) = 2.236067977
+    float c = 0, s = 0;
+    QR::GivensRotation(2.0f, 1.0f, c, s);
+    REQUIRE_THAT(c, Catch::Matchers::WithinRel(0.894427191f, 1e-6f));
+    REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.447213595f, 1e-6f));
+    REQUIRE_THAT(c * 2.0f + s * 1.0f,
+                 Catch::Matchers::WithinRel(2.236067977f, 1e-6f));
+    REQUIRE_THAT(-s * 2.0f + c * 1.0f, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+  }
+
+  {
+    // Reference: hypot(3, 4) = 5 exactly
+    float c = 0, s = 0;
+    QR::GivensRotation(3.0f, 4.0f, c, s);
+    REQUIRE_THAT(c, Catch::Matchers::WithinRel(0.6f, 1e-6f));
+    REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.8f, 1e-6f));
+    REQUIRE_THAT(c * 3.0f + s * 4.0f, Catch::Matchers::WithinRel(5.0f, 1e-6f));
+    REQUIRE_THAT(-s * 3.0f + c * 4.0f, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+  }
+
+  {
+    // Pure second component: c = 0, s = 1
+    float c = 1, s = 1;
+    QR::GivensRotation(0.0f, 5.0f, c, s);
+    REQUIRE_THAT(c, Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+    REQUIRE_THAT(s, Catch::Matchers::WithinRel(1.0f, 1e-6f));
+  }
+
+  {
+    // Zero vector: identity rotation
+    float c = 0, s = 0;
+    QR::GivensRotation(0.0f, 0.0f, c, s);
+    REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-7f));
+    REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+  }
+
+  {
+    // Negative first component preserves the sign of c
+    float c = 0, s = 0;
+    QR::GivensRotation(-2.0f, 1.0f, c, s);
+    REQUIRE_THAT(c, Catch::Matchers::WithinRel(-0.894427191f, 1e-6f));
+    REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.447213595f, 1e-6f));
+    REQUIRE_THAT(-s * -2.0f + c * 1.0f, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+  }
+}
+
+// ============================================================================
+// TEST 2: ApplyRotationBothSides (similarity A <- G A G^T)
+// ============================================================================
+TEST_CASE("QR Building Block: ApplyRotationBothSides", "[Matrix][QR]") {
+  // Reference (numpy, float64): A = [[2,1,0],[1,3,1],[0,1,4]], i = 0,
+  // Givens(2,1) -> G A G^T =
+  // [[ 3.0,         1.0,           0.447213595],
+  //  [ 1.0,         2.0,           0.894427191],
+  //  [ 0.447213595, 0.894427191,   4.0]]
+  // (Note: G A G^T with G zeroing (2,1) sends the A[0][1] coupling into the
+  // (0,2) corner, NOT into the subdiagonal -- the subdiagonal-zeroing happens
+  // in the QR chase context where the bulge column has the right shape.)
+  {
+    Matrix<3, 3> A{2, 1, 0, 1, 3, 1, 0, 1, 4};
+    float c = 0.894427191f, s = 0.447213595f;
+
+    QR::ApplyRotationBothSides(A, 0, c, s);
+
+    REQUIRE_THAT(A.Get(0, 0), Catch::Matchers::WithinRel(3.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(0, 1), Catch::Matchers::WithinRel(1.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(0, 2),
+                 Catch::Matchers::WithinRel(0.447213595f, 1e-5f));
+    REQUIRE_THAT(A.Get(1, 1), Catch::Matchers::WithinRel(2.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(1, 2),
+                 Catch::Matchers::WithinRel(0.894427191f, 1e-5f));
+    REQUIRE_THAT(A.Get(2, 2), Catch::Matchers::WithinRel(4.0f, 1e-5f));
+
+    // Symmetry must be preserved exactly in both triangles
+    for (uint8_t i = 0; i < 3; i++)
+      for (uint8_t j = 0; j < 3; j++)
+        REQUIRE(A.Get(i, j) == A.Get(j, i));
+  }
+
+  // Same check at i = 1.
+  // Reference (numpy, float64): B = [[5,0,1],[0,6,2],[1,2,7]], i = 1,
+  // Givens(6,2) -> G B G^T =
+  // [[ 5.0,          0.316227766,  0.948683298],
+  //  [ 0.316227766,   7.3,          1.9],
+  //  [ 0.948683298,   1.9,          5.7]]
+  {
+    Matrix<3, 3> B{5, 0, 1, 0, 6, 2, 1, 2, 7};
+    float c = 0.948683298f, s = 0.316227766f;
+
+    QR::ApplyRotationBothSides(B, 1, c, s);
+
+    REQUIRE_THAT(B.Get(0, 0), Catch::Matchers::WithinRel(5.0f, 1e-5f));
+    REQUIRE_THAT(B.Get(0, 1),
+                 Catch::Matchers::WithinRel(0.316227766f, 1e-5f));
+    REQUIRE_THAT(B.Get(0, 2),
+                 Catch::Matchers::WithinRel(0.948683298f, 1e-5f));
+    REQUIRE_THAT(B.Get(1, 1), Catch::Matchers::WithinRel(7.3f, 1e-5f));
+    REQUIRE_THAT(B.Get(1, 2), Catch::Matchers::WithinRel(1.9f, 1e-5f));
+    REQUIRE_THAT(B.Get(2, 2), Catch::Matchers::WithinRel(5.7f, 1e-5f));
+
+    for (uint8_t i = 0; i < 3; i++)
+      for (uint8_t j = 0; j < 3; j++)
+        REQUIRE(B.Get(i, j) == B.Get(j, i));
+  }
+
+  // Identity rotation leaves the matrix unchanged
+  {
+    Matrix<3, 3> C{1, 2, 3, 2, 4, 5, 3, 5, 6};
+    QR::ApplyRotationBothSides(C, 1, 1.0f, 0.0f);
+    REQUIRE(C.Get(0, 0) == 1.0f);
+    REQUIRE(C.Get(0, 1) == 2.0f);
+    REQUIRE(C.Get(0, 2) == 3.0f);
+    REQUIRE(C.Get(1, 1) == 4.0f);
+    REQUIRE(C.Get(1, 2) == 5.0f);
+    REQUIRE(C.Get(2, 2) == 6.0f);
+  }
+
+  // Spectrum invariants (trace, Frobenius norm) are preserved. (c, s)
+  // must be a unit vector for G A G^T to be a similarity transform.
+  {
+    Matrix<3, 3> D{1, 2, 3, 2, 5, 8, 3, 8, 9};
+    float tr = trace3(D);
+    float fn = frob(D);
+    float c = 0.6f, s = 0.8f;
+    QR::ApplyRotationBothSides(D, 0, c, s);
+    REQUIRE_THAT(trace3(D), Catch::Matchers::WithinRel(tr, 1e-5f));
+    REQUIRE_THAT(frob(D), Catch::Matchers::WithinRel(fn, 1e-5f));
+  }
+}
+// ============================================================================
+// TEST 3: ApplyRotationToVectors (V <- V G^T)
+// ============================================================================
+TEST_CASE("QR Building Block: ApplyRotationToVectors", "[Matrix][QR]") {
+  // V = I, i = 0, Givens(2,1): V <- I * G^T with G^T = [[c, -s], [s, c]] =
+  // [[ c,  -s,  0],
+  //  [ s,   c,  0],
+  //  [ 0,   0,  1]]
+  {
+    Matrix<3, 3> V{0};
+    V[0][0] = 1;
+    V[1][1] = 1;
+    V[2][2] = 1;
+    float c = 0.894427191f, s = 0.447213595f;
+
+    QR::ApplyRotationToVectors(V, 0, c, s);
+
+    REQUIRE_THAT(V.Get(0, 0), Catch::Matchers::WithinRel(0.894427191f, 1e-6f));
+    REQUIRE_THAT(V.Get(0, 1), Catch::Matchers::WithinRel(-0.447213595f, 1e-6f));
+    REQUIRE_THAT(V.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+    REQUIRE_THAT(V.Get(1, 0), Catch::Matchers::WithinRel(0.447213595f, 1e-6f));
+    REQUIRE_THAT(V.Get(1, 1), Catch::Matchers::WithinRel(0.894427191f, 1e-6f));
+    REQUIRE_THAT(V.Get(1, 2), Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+    REQUIRE_THAT(V.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+    REQUIRE_THAT(V.Get(2, 1), Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+    REQUIRE_THAT(V.Get(2, 2), Catch::Matchers::WithinRel(1.0f, 1e-7f));
+
+    // Product of rotations must stay orthogonal
+    REQUIRE(isOrthogonal(V));
+  }
+
+  // Two successive rotations accumulate (V <- V G1^T G2^T)
+  // Reference (numpy, float64):
+  // [[ 0.894427191, -0.424264069,  0.141421356],
+  //  [ 0.447213595,  0.848528137, -0.282842712],
+  //  [ 0.0,          0.316227766,  0.948683298]]
+  {
+    Matrix<3, 3> V{0};
+    V[0][0] = 1;
+    V[1][1] = 1;
+    V[2][2] = 1;
+    QR::ApplyRotationToVectors(V, 0, 0.894427191f, 0.447213595f);
+    QR::ApplyRotationToVectors(V, 1, 0.948683298f, 0.316227766f);
+    REQUIRE(isOrthogonal(V));
+    // Column 0 was only touched by the first rotation
+    REQUIRE_THAT(V.Get(0, 0), Catch::Matchers::WithinRel(0.894427191f, 1e-5f));
+    REQUIRE_THAT(V.Get(1, 0), Catch::Matchers::WithinRel(0.447213595f, 1e-5f));
+    REQUIRE_THAT(V.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+    REQUIRE_THAT(V.Get(0, 1), Catch::Matchers::WithinRel(-0.424264069f, 1e-5f));
+    REQUIRE_THAT(V.Get(0, 2), Catch::Matchers::WithinRel(0.141421356f, 1e-5f));
+    REQUIRE_THAT(V.Get(1, 2), Catch::Matchers::WithinRel(-0.282842712f, 1e-5f));
+    REQUIRE_THAT(V.Get(2, 1), Catch::Matchers::WithinRel(0.316227766f, 1e-5f));
+    REQUIRE_THAT(V.Get(2, 2), Catch::Matchers::WithinRel(0.948683298f, 1e-5f));
+  }
+}
+
+// ============================================================================
+// TEST 4: WilkinsonShift
+// ============================================================================
+TEST_CASE("QR Building Block: WilkinsonShift", "[Matrix][QR]") {
+  // mu = (a+d)/2 - sign(a-d) * sqrt(((a-d)/2)^2 + b^2)
+  // Reference: eigenvalues of [[2,1],[1,4]] are 1.5858, 4.4142; closest
+  // to d = 4 is 4.414213562.
+  REQUIRE_THAT(QR::WilkinsonShift(2.0f, 1.0f, 4.0f),
+               Catch::Matchers::WithinRel(4.414213562f, 1e-6f));
+
+  // [[5,2],[2,1]]: eigenvalues 0.1716, 5.8284; closest to d = 1 is 0.171572875
+  REQUIRE_THAT(QR::WilkinsonShift(5.0f, 2.0f, 1.0f),
+               Catch::Matchers::WithinRel(0.171572875f, 1e-5f));
+
+  // Zero off-diagonal: returns d itself (sign(0) = +1 picks d, not a)
+  REQUIRE_THAT(QR::WilkinsonShift(3.0f, 0.0f, 7.0f),
+               Catch::Matchers::WithinRel(7.0f, 1e-7f));
+  REQUIRE_THAT(QR::WilkinsonShift(7.0f, 0.0f, 3.0f),
+               Catch::Matchers::WithinRel(3.0f, 1e-7f));
+
+  // a == d: shift is the larger-magnitude off-diagonal combination
+  // [[1,3],[3,1]]: eigenvalues -2, 4; closest to d = 1 is -2
+  REQUIRE_THAT(QR::WilkinsonShift(1.0f, 3.0f, 1.0f),
+               Catch::Matchers::WithinRel(-2.0f, 1e-6f));
+}
+
+// ============================================================================
+// TEST 5: Solve2x2Eigen
+// ============================================================================
+TEST_CASE("QR Building Block: Solve2x2Eigen", "[Matrix][QR]") {
+  // Symmetric block [[2,1],[1,3]]:
+  // eigenvalues 1.381966011, 3.618033989;
+  // eigenvector of 3.618033989 is +/- (0.525731112, 0.850650808)
+  {
+    Matrix<2, 2> A{2, 1, 1, 3};
+    float lHi = 0, lLo = 0, c = 0, s = 0;
+    QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s);
+
+    REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(3.618033989f, 1e-6f));
+    REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(1.381966011f, 1e-6f));
+    REQUIRE(matchesAbs(c, 0.525731112f));
+    REQUIRE(matchesAbs(s, 0.850650808f));
+
+    // Residual: A * vHi = lHi * vHi with vHi = (c, s)
+    REQUIRE_THAT(c * 2.0f + s * 1.0f,
+                 Catch::Matchers::WithinRel(lHi * c, 1e-5f));
+    REQUIRE_THAT(c * 1.0f + s * 3.0f,
+                 Catch::Matchers::WithinRel(lHi * s, 1e-5f));
+    // Second eigenvector vLo = (-s, c)
+    REQUIRE_THAT(-s * 2.0f + c * 1.0f,
+                 Catch::Matchers::WithinRel(lLo * -s, 1e-5f));
+    REQUIRE_THAT(-s * 1.0f + c * 3.0f,
+                 Catch::Matchers::WithinRel(lLo * c, 1e-5f));
+  }
+
+  // Nonsymmetric block [[1,2],[3,4]] (used by the N == 2 entry point):
+  // eigenvalues 5.372281323, -0.372281323;
+  // eigenvector of 5.372281323 is +/- (0.415973558, 0.909376709)
+  {
+    Matrix<2, 2> A{1, 2, 3, 4};
+    float lHi = 0, lLo = 0, c = 0, s = 0;
+    QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s);
+
+    REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(5.372281323f, 1e-6f));
+    REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(-0.372281323f, 1e-6f));
+    REQUIRE(matchesAbs(c, 0.415973558f));
+    REQUIRE(matchesAbs(s, 0.909376709f));
+
+    // Both-row residual with vHi = (c, s): A v = l v
+    REQUIRE_THAT(c * 1.0f + s * 2.0f,
+                 Catch::Matchers::WithinRel(lHi * c, 1e-5f));
+    REQUIRE_THAT(c * 3.0f + s * 4.0f,
+                 Catch::Matchers::WithinRel(lHi * s, 1e-5f));
+  }
+
+  // Diagonal blocks: eigenvectors are coordinate vectors
+  {
+    Matrix<2, 2> A{5, 0, 0, 2};
+    float lHi = 0, lLo = 0, c = 0, s = 0;
+    QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s);
+    REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(5.0f, 1e-7f));
+    REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(2.0f, 1e-7f));
+    REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-7f));
+    REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+
+    A = Matrix<2, 2>{2, 0, 0, 5};
+    QR::Solve2x2Eigen(A, 0, lHi, lLo, c, s);
+    REQUIRE_THAT(lHi, Catch::Matchers::WithinRel(5.0f, 1e-7f));
+    REQUIRE_THAT(lLo, Catch::Matchers::WithinRel(2.0f, 1e-7f));
+    REQUIRE_THAT(c, Catch::Matchers::WithinAbs(0.0f, 1e-7f));
+    REQUIRE_THAT(s, Catch::Matchers::WithinRel(1.0f, 1e-7f));
+  }
+}
+// ============================================================================
+// TEST 6: Deflate
+// ============================================================================
+TEST_CASE("QR Building Block: Deflate", "[Matrix][QR]") {
+  // subdiag[0] = 1e-9 <= 1e-6 * (|2| + |3|) = 5e-6  -> deflated
+  // subdiag[1] = 0.5  > 1e-6 * (|3| + |4|) = 7e-6  -> kept
+  {
+    Matrix<3, 3> A{2, 1e-9f, 0, 1e-9f, 3, 0.5f, 0, 0.5f, 4};
+    QR::Deflate(A, 0, 2, 1e-6f);
+
+    REQUIRE(A.Get(1, 0) == 0.0f);
+    REQUIRE(A.Get(0, 1) == 0.0f);
+    REQUIRE_THAT(A.Get(2, 1), Catch::Matchers::WithinRel(0.5f, 1e-7f));
+    REQUIRE_THAT(A.Get(1, 2), Catch::Matchers::WithinRel(0.5f, 1e-7f));
+    // Diagonals untouched
+    REQUIRE_THAT(A.Get(0, 0), Catch::Matchers::WithinRel(2.0f, 1e-7f));
+    REQUIRE_THAT(A.Get(1, 1), Catch::Matchers::WithinRel(3.0f, 1e-7f));
+    REQUIRE_THAT(A.Get(2, 2), Catch::Matchers::WithinRel(4.0f, 1e-7f));
+  }
+
+  // Nothing deflated when all subdiagonals are well above tolerance
+  {
+    Matrix<3, 3> A{2, 0.1f, 0, 0.1f, 3, 0.2f, 0, 0.2f, 4};
+    QR::Deflate(A, 0, 2, 1e-6f);
+    REQUIRE_THAT(A.Get(1, 0), Catch::Matchers::WithinRel(0.1f, 1e-7f));
+    REQUIRE_THAT(A.Get(2, 1), Catch::Matchers::WithinRel(0.2f, 1e-7f));
+  }
+}
+
+// ============================================================================
+// TEST 7: Tridiagonalize
+// ============================================================================
+TEST_CASE("QR Building Block: Tridiagonalize", "[Matrix][QR]") {
+  // 4x4 symmetric with a full (0,3) corner coupling
+  {
+    Matrix<4, 4> A{2, 1, 0, 1, 1, 3, 1, 0, 0, 1, 4, 1, 1, 0, 1, 5};
+    Matrix<4, 4> Aorig = A;
+    Matrix<4, 4> U{0};
+
+    QR::Tridiagonalize(A, U);
+
+    // Off-tridiagonal entries must be zero up to float32 roundoff (the
+    // Givens zeroing cancels only in exact arithmetic; residuals are
+    // ~1e-7 for O(1) entries).
+    REQUIRE_THAT(A.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(0, 3), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(1, 3), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+    REQUIRE_THAT(A.Get(3, 1), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+
+    // Symmetry preserved exactly
+    for (uint8_t i = 0; i < 4; i++)
+      for (uint8_t j = 0; j < 4; j++)
+        REQUIRE(A.Get(i, j) == A.Get(j, i));
+
+    // U must be orthogonal
+    REQUIRE(isOrthogonal(U));
+
+    // Reconstruction: U * A_tri * U^T == Aorig (absolute check for
+    // originally-zero entries: WithinRel has no absolute fallback there)
+    Matrix<4, 4> UAt{};
+    U.Mult(A, UAt);
+    Matrix<4, 4> UAtU{};
+    UAt.Mult(U.Transpose(), UAtU);
+    for (uint8_t i = 0; i < 4; i++)
+      for (uint8_t j = 0; j < 4; j++) {
+        float actual = UAtU.Get(i, j);
+        float expected = Aorig.Get(i, j);
+        if (fabsf(expected) < 1e-3f)
+          REQUIRE_THAT(actual, Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+        else
+          REQUIRE_THAT(actual,
+                       Catch::Matchers::WithinRel(expected, 1e-5f));
+      }
+
+    // Spectrum invariants match the original
+    {
+      float tr0 = Aorig.Get(0, 0) + Aorig.Get(1, 1) + Aorig.Get(2, 2) +
+                  Aorig.Get(3, 3);
+      float tr1 = A.Get(0, 0) + A.Get(1, 1) + A.Get(2, 2) + A.Get(3, 3);
+      REQUIRE_THAT(tr1, Catch::Matchers::WithinRel(tr0, 1e-6f));
+      REQUIRE_THAT(frob(A), Catch::Matchers::WithinRel(frob(Aorig), 1e-6f));
+    }
+
+    // Eigenvalues of the tridiagonal match the original (scipy reference):
+    // 6.0, 4.0, 3.0, 1.0
+    {
+      Matrix<4, 1> vals{};
+      Matrix<4, 4> vecs{};
+      QR::EigenQR(A, vecs, vals, 10000, 1e-6f);
+      REQUIRE_THAT(vals[0][0], Catch::Matchers::WithinRel(6.0f, 1e-4f));
+      REQUIRE_THAT(vals[1][0], Catch::Matchers::WithinRel(4.0f, 1e-4f));
+      REQUIRE_THAT(vals[2][0], Catch::Matchers::WithinRel(3.0f, 1e-4f));
+      REQUIRE_THAT(vals[3][0], Catch::Matchers::WithinRel(1.0f, 1e-4f));
+    }
+  }
+
+  // 5x5 symmetric
+  {
+    Matrix<5, 5> A{3, 1, 0, 0, 1, 1, 4, 1, 0, 0, 0, 1, 5, 1, 0, 0, 0, 1, 6, 1,
+                   1, 0, 0, 1, 7};
+    Matrix<5, 5> Aorig = A;
+    Matrix<5, 5> U{0};
+
+    QR::Tridiagonalize(A, U);
+
+    // All |i - j| >= 2 entries zero up to float32 roundoff
+    for (uint8_t i = 0; i < 5; i++)
+      for (uint8_t j = 0; j < 5; j++)
+        if (i > j + 1 || j > i + 1)
+          REQUIRE_THAT(A.Get(i, j), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+
+    REQUIRE(isOrthogonal(U));
+
+    Matrix<5, 5> UAt{};
+    U.Mult(A, UAt);
+    Matrix<5, 5> UAtU{};
+    UAt.Mult(U.Transpose(), UAtU);
+    for (uint8_t i = 0; i < 5; i++)
+      for (uint8_t j = 0; j < 5; j++) {
+        float actual = UAtU.Get(i, j);
+        float expected = Aorig.Get(i, j);
+        if (fabsf(expected) < 1e-3f)
+          REQUIRE_THAT(actual, Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+        else
+          REQUIRE_THAT(actual,
+                       Catch::Matchers::WithinRel(expected, 1e-5f));
+      }
+  }
+
+  // Already tridiagonal: U must come out as the identity
+  {
+    Matrix<3, 3> A{1, 2, 0, 2, 5, 2, 0, 2, 9};
+    Matrix<3, 3> U{0};
+    QR::Tridiagonalize(A, U);
+    for (uint8_t i = 0; i < 3; i++)
+      for (uint8_t j = 0; j < 3; j++) {
+        float expected = (i == j) ? 1.0f : 0.0f;
+        REQUIRE_THAT(U.Get(i, j), Catch::Matchers::WithinAbs(expected, 1e-7f));
+      }
+  }
+}
+// ============================================================================
+// TEST 8: One full shifted QR step (integration of the blocks)
+// ============================================================================
+TEST_CASE("QR Building Block: Full Shifted QR Step", "[Matrix][QR]") {
+  // One Wilkinson-shifted QR step on the whole 3x3 block is a similarity
+  // transform, so all spectrum invariants (trace, sum of principal 2x2
+  // minors, determinant) must be preserved.
+  //
+  // A = [[1,2,3],[2,5,8],[3,8,9]]:  tr = 15, e2 = -18, det = -4
+  {
+    Matrix<3, 3> A{1, 2, 3, 2, 5, 8, 3, 8, 9};
+    float tr0 = trace3(A);   // 15
+    float e20 = e2_3x3(A);   // -18
+    float det0 = det3(A);    // -4
+
+    // mu from the trailing 2x2 [[5,8],[8,9]]: eigenvalues
+    // -1.246211251, 15.246211251; closest to d = 9 is 15.246211251 (Wilkinson)
+    float mu = QR::WilkinsonShift(A.Get(1, 1), A.Get(2, 1), A.Get(2, 2));
+    REQUIRE_THAT(mu, Catch::Matchers::WithinRel(15.246211251f, 1e-5f));
+
+    for (uint8_t i = 0; i < 3; i++)
+      A[i][i] -= mu;
+
+    // Bulge chase: rotations on (0,1) then (1,2)
+    float c = 0, s = 0;
+    QR::GivensRotation(A.Get(0, 0), A.Get(1, 0), c, s);
+    QR::ApplyRotationBothSides(A, 0, c, s);
+    QR::GivensRotation(A.Get(1, 1), A.Get(2, 1), c, s);
+    QR::ApplyRotationBothSides(A, 1, c, s);
+
+    for (uint8_t i = 0; i < 3; i++)
+      A[i][i] += mu;
+
+    // Symmetry preserved
+    for (uint8_t i = 0; i < 3; i++)
+      for (uint8_t j = 0; j < 3; j++)
+        REQUIRE(A.Get(i, j) == A.Get(j, i));
+
+    // Spectrum invariants preserved
+    REQUIRE_THAT(trace3(A), Catch::Matchers::WithinRel(tr0, 1e-5f));
+    REQUIRE_THAT(e2_3x3(A), Catch::Matchers::WithinRel(e20, 1e-5f));
+    REQUIRE_THAT(det3(A), Catch::Matchers::WithinRel(det0, 1e-5f));
+  }
+
+  // For TRIDIAGONAL input a single step keeps the tridiagonal structure
+  {
+    Matrix<3, 3> T{1, 2, 0, 2, 5, 2, 0, 2, 9};
+    float mu = QR::WilkinsonShift(T.Get(1, 1), T.Get(2, 1), T.Get(2, 2));
+    for (uint8_t i = 0; i < 3; i++)
+      T[i][i] -= mu;
+    float c = 0, s = 0;
+    QR::GivensRotation(T.Get(0, 0), T.Get(1, 0), c, s);
+    QR::ApplyRotationBothSides(T, 0, c, s);
+    QR::GivensRotation(T.Get(1, 1), T.Get(2, 1), c, s);
+    QR::ApplyRotationBothSides(T, 1, c, s);
+    for (uint8_t i = 0; i < 3; i++)
+      T[i][i] += mu;
+
+    // Corners must vanish up to float32 roundoff: tridiagonal form
+    // maintained. The cancellation is exact in exact arithmetic (the
+    // corner is s1*a - c1*b times a factor, and Givens gives s1*a = c1*b),
+    // so the residual is pure rounding, ~1e-6 for O(1) entries.
+    REQUIRE_THAT(T.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+    REQUIRE_THAT(T.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
+  }
+}

unit-tests/qr-reference-values.py

@@ -0,0 +1,246 @@
+#!/usr/bin/env python3
+"""
+Reference values for the QR eigen-decomposition building block tests
+(unit-tests/qr-build-blocks-tests.cpp). Run this to verify/implement the
+C++ implementation in src/QR.hpp / src/QR.cpp against numpy/scipy.
+
+Conventions (match the C++ exactly):
+  * Givens zeroing rotation: G = [[c, s], [-s, c]], c = x/r, s = y/r,
+    r = hypot(x, y).  G * (x, y)^T = (r, 0)^T.
+  * Similarity transform: A <- G A G^T (ApplyRotationBothSides).
+  * Eigenvector accumulation: V <- V G^T (ApplyRotationToVectors).
+    Vblock in the 2x2 closed form is [[c, -s], [s, c]] (same shape as G^T).
+  * Tridiagonalization: bottom-up Givens (i = N-2 down to k+1 per column k).
+  * Shifted QR loop: Wilkinson shift mu from the trailing 2x2, chase on the
+    trailing unreduced block [lo, hi], deflate by relative tolerance, peel
+    exact-zero subdiagonals, 2x2 closed-form termination.
+  * Pipeline: M0 = U * Mtri * U^T and Mtri = V * D * V^T  =>
+    eigenvectors of M0 = U * V (columns), eigenvalues = diag(D).
+
+Usage: python3 qr-reference-values.py
+"""
+
+import numpy as np
+import scipy.linalg as sla
+
+np.set_printoptions(precision=9, linewidth=120)
+
+
+def givens(x, y):
+    """c = x/r, s = y/r with r = hypot(x, y)."""
+    r = np.hypot(x, y)
+    if r == 0.0:
+        return 1.0, 0.0
+    return x / r, y / r
+
+
+def rot(n, i, c, s):
+    """G = I with [[c, s], [-s, c]] embedded at (i, i+1)."""
+    G = np.eye(n)
+    G[i:i + 2, i:i + 2] = np.array([[c, s], [-s, c]])
+    return G
+
+
+def tridiagonalize(M0):
+    """Bottom-up Givens tridiagonalization. Returns (Mtri, U) with
+    M0 = U Mtri U^T."""
+    n = len(M0)
+    M = M0.copy()
+    U = np.eye(n)
+    for k in range(n - 2):
+        for i in range(n - 2, k, -1):
+            c, s = givens(M[i, k], M[i + 1, k])
+            G = rot(n, i, c, s)
+            M = G @ M @ G.T
+            U = U @ G.T
+    return M, U
+
+
+def wilkinson(a, b, d):
+    """Eigenvalue of [[a, b], [b, d]] closest to d."""
+    delta = 0.5 * (a - d)
+    spread = np.sqrt(delta * delta + b * b)
+    return 0.5 * (a + d) - (spread if delta >= 0 else -spread)
+
+
+def solve2x2(A, lo):
+    """Closed form for the block at (lo, lo+1): (lHi, lLo, c, s) with
+    vHi = (c, s), vLo = (-s, c)."""
+    a = A[lo, lo]
+    b = A[lo, lo + 1]
+    e = A[lo + 1, lo]
+    d = A[lo + 1, lo + 1]
+    tr = a + d
+    det = a * d - b * e
+    disc = max(0.0, tr * tr - 4 * det)
+    lhi = 0.5 * (tr + np.sqrt(disc))
+    llo = 0.5 * (tr - np.sqrt(disc))
+    if b != 0.0:
+        v1 = lhi - a
+        nn = np.hypot(b, v1)
+        c, s = b / nn, v1 / nn
+    elif a >= d:
+        c, s = 1.0, 0.0
+    else:
+        c, s = 0.0, 1.0
+    return lhi, llo, c, s
+
+
+def eigenqr(M0, tol=1e-12, max_iter=100000):
+    """Full pipeline mirroring QR::EigenQR. Returns (eigs, W) where W has
+    the eigenvectors of M0 as columns."""
+    n = len(M0)
+    if n == 2:
+        l1, l2, c, s = solve2x2(M0, 0)
+        return np.array([l1, l2]), np.array([[c, -s], [s, c]])
+    M, U = tridiagonalize(M0)
+    V = np.eye(n)
+    hi = n - 1
+    for _ in range(max_iter):
+        # deflate: zero tiny subdiagonals (relative test)
+        for i in range(hi):
+            t = M[i + 1, i]
+            scale = abs(M[i, i]) + abs(M[i + 1, i + 1])
+            if abs(t) <= tol * scale:
+                M[i + 1, i] = M[i, i + 1] = 0.0
+        # peel exact-zero trailing subdiagonals
+        while hi > 0 and M[hi, hi - 1] == 0.0:
+            hi -= 1
+        if hi == 0:
+            break
+        # find start of trailing unreduced block
+        lo = hi
+        for i in range(hi - 1, -1, -1):
+            if M[i + 1, i] == 0.0:
+                break
+            lo = i
+        if lo + 1 == hi:
+            # closed-form 2x2 termination: set diagonal, fold Vblock in
+            l1, l2, c, s = solve2x2(M, lo)
+            Vb = np.eye(n)
+            Vb[lo:lo + 2, lo:lo + 2] = np.array([[c, -s], [s, c]])
+            V = V @ Vb
+            M[lo, lo] = l1
+            M[lo + 1, lo + 1] = l2
+            M[lo + 1, lo] = M[lo, lo + 1] = 0.0
+            if lo == 0:
+                break
+            hi = lo - 1
+            continue
+        # full shifted step on [lo, hi] (shift applies to the active block)
+        mu = wilkinson(M[hi - 1, hi - 1], M[hi, hi - 1], M[hi, hi])
+        diag = M.diagonal().copy()
+        diag[lo:hi + 1] -= mu
+        np.fill_diagonal(M, diag)
+        c, s = givens(M[lo, lo], M[lo + 1, lo])
+        G = rot(n, lo, c, s)
+        M = G @ M @ G.T
+        V = V @ G.T
+        for i in range(lo + 1, hi):
+            c, s = givens(M[i, i], M[i + 1, i])
+            G = rot(n, i, c, s)
+            M = G @ M @ G.T
+            V = V @ G.T
+        diag = M.diagonal().copy()
+        diag[lo:hi + 1] += mu
+        np.fill_diagonal(M, diag)
+
+    eigs = np.diag(M).astype(float)
+    order = np.argsort(eigs)[::-1]  # descending, like the C++ test harness
+    eigs = eigs[order]
+    W = U @ V
+    W = W[:, order]
+    return eigs, W
+
+
+def report(name, val, ref=None, tol=1e-6):
+    ok = "OK " if ref is None or np.allclose(val, ref, rtol=tol, atol=tol) else "FAIL"
+    print(f"[{ok}] {name} = {val}")
+    if ref is not None:
+        print(f"        scipy/numpy ref = {ref}")
+
+
+def main():
+    print("=== TEST 1: GivensRotation ===")
+    c, s = givens(2.0, 1.0)
+    print(f"  c = {c}  s = {s}")
+    # G * (x, y)^T = (r, 0)^T:  G = [[c, s], [-s, c]]
+    assert abs(c * 2 + s * 1 - np.sqrt(5)) < 1e-15
+    assert abs(-s * 2 + c * 1) < 1e-15
+
+    print("\n=== TEST 2: ApplyRotationBothSides A <- G A G^T ===")
+    A = np.array([[3.0, 4.0, 5.0], [6.0, 7.0, 8.0], [9.0, 10.0, 11.0]])
+    G = rot(3, 0, 0.6, 0.8)
+    B = G @ A @ G.T
+    print(B)
+
+    A = np.array([[5.0, 0.0, 1.0], [0.0, 6.0, 2.0], [1.0, 2.0, 7.0]])
+    c, s = givens(6.0, 2.0)
+    G = rot(3, 1, c, s)
+    B = G @ A @ G.T
+    print(B)
+
+    print("\n=== TEST 3: V accumulation V <- V G^T ===")
+    V = np.eye(3)
+    G = rot(3, 0, 0.894427191, 0.447213595)
+    V = V @ G.T
+    print(V)
+    V2 = V @ rot(3, 1, 0.848874681, 0.528748047).T
+    print(V2)
+
+    print("\n=== TEST 4: Solve2x2Eigen ===")
+    for A in (np.array([[5.0, 8.0], [8.0, 9.0]]), np.array([[1.0, 2.0], [3.0, 4.0]])):
+        l1, l2, c, s = solve2x2(A, 0)
+        ref = np.linalg.eigvalsh(A) if np.allclose(A, A.T) else np.linalg.eigvals(A)
+        print(f"  A={A.ravel()}  lHi={l1} lLo={l2}  c={c} s={s}  ref={np.sort(ref)[::-1]}")
+
+    print("\n=== TEST 8: WilkinsonShift ===")
+    print(f"  W(5, 8, 9) = {wilkinson(5, 8, 9)}")
+    print(f"  W(4, 2, 7) = {wilkinson(4, 2, 7)}")
+    print(f"  W(9, 2, 5) = {wilkinson(9, 2, 5)}")
+
+    print("\n=== TEST 8b: one full shifted chase step on tridiagonal 3x3 ===")
+    T = np.array([[1.0, 2.0, 0.0], [2.0, 5.0, 2.0], [0.0, 2.0, 9.0]])
+    mu = wilkinson(5, 2, 9)
+    M = T - mu * np.eye(3)
+    c, s = givens(M[0, 0], M[1, 0])
+    M = rot(3, 0, c, s) @ M @ rot(3, 0, c, s).T
+    c, s = givens(M[1, 1], M[2, 1])
+    M = rot(3, 1, c, s) @ M @ rot(3, 1, c, s).T
+    M = M + mu * np.eye(3)
+    print(f"  mu = {mu}")
+    print(M)
+    print(f"  corners: {M[0, 2]}, {M[2, 0]} (exact-arithmetic zeros)")
+    print(f"  trace {M.trace():.15f} (was {T.trace()})")
+
+    print("\n=== TEST 7: Tridiagonalize ===")
+    M4 = np.array([[2.0, 1, 0, 1], [1, 3, 1, 0], [0, 1, 4, 1], [1, 0, 1, 5]])
+    M, U = tridiagonalize(M4)
+    print("  M4 tridiagonalized:\n", M)
+    print(f"  reconstruction U M U^T == M4: {np.allclose(U @ M @ U.T, M4, atol=1e-9)}")
+    M5 = np.array([[3.0, 1, 0, 0, 1], [1, 4, 1, 0, 0], [0, 1, 5, 1, 0],
+                   [0, 0, 1, 6, 1], [1, 0, 0, 1, 7]])
+    M, U = tridiagonalize(M5)
+    print("  M5 tridiagonalized:\n", M)
+    print(f"  reconstruction: {np.allclose(U @ M @ U.T, M5, atol=1e-9)}")
+
+    print("\n=== End-to-end: random symmetric vs scipy.linalg.eigh ===")
+    rng = np.random.default_rng(12345)
+    worst = 0.0
+    for n in range(3, 9):
+        M0 = rng.normal(size=(n, n))
+        M0 = (M0 + M0.T) / 2
+        eigs, W = eigenqr(M0.astype(float))
+        ref = sla.eigh(M0)
+        e_err = np.max(np.abs(np.sort(eigs) - ref[0]))
+        resid = np.linalg.norm(W @ np.diag(eigs) @ W.T - M0)
+        ortho = np.linalg.norm(W.T @ W - np.eye(n))
+        print(f"  n={n}: eigs_err={e_err:.2e}  resid={resid:.2e}  ortho={ortho:.2e}")
+        worst = max(worst, e_err, resid, ortho)
+    print(f"\nworst over all n: {worst:.2e}")
+    assert worst < 1e-10, "end-to-end reference FAILED"
+    print("ALL REFERENCES OK")
+
+
+if __name__ == "__main__":
+    main()