Fixes for control systems AND SVD and QR decomposition (#9)

Reviewed-on: #9
Co-authored-by: Cynopolis <quinn.henthorne@gmail.com>

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
@@ -74,3 +57,39 @@ 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,35 @@
 #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
+
+// Forward-declare SVD::SVD so the Matrix::SVD implementation below can call
+// it even when Matrix.cpp is pulled in through SVD.hpp's own include chain
+// (SVD.hpp -> Matrix.hpp -> Matrix.cpp), where the SVD namespace has not
+// been declared yet at this point. If we are not already inside that chain,
+// pull in the full SVD library so its template definition is available.
+namespace SVD {
+template <uint8_t rows, uint8_t columns>
+void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
+         Matrix<columns, 1> &sigma, Matrix<columns, columns> &Vt);
+}
+#ifndef SVD_H_
+#include "SVD.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 +600,24 @@ 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) {
+template <uint8_t rows, uint8_t columns>
-    Matrix<rows, rows> Q, R;
+void Matrix<rows, columns>::SVD(Matrix<rows, columns> &U,
-
+                                Matrix<columns, 1> &sigma,
-    // // QR shift lets us "attack" the first diagonal to speed up the algorithm
+                                Matrix<columns, columns> &Vt) const {
-    // shift = Matrix<rows, rows>::Identity() * Ak[rows - 1][rows - 1];
+  // Delegate to the SVD library (see src/SVD.hpp for the algorithm and
-    (Ak - shift).QRDecomposition(Q, R);
+  // conventions). NB: the fully-qualified ::SVD is required here — inside
-    Ak = R * Q + shift;
+  // this member the unqualified name SVD refers to this method, which
-    QQ = QQ * Q;
+  // would shadow the namespace in a qualified lookup. SVD::SVD takes its
-
+  // input by non-const reference but does not modify it; pass a copy so
-    // Check convergence: off-diagonal norm
+  // the const-ness of *this is preserved.
-    float offDiagSum = 0.0f;
+  Matrix<rows, columns> A = *this;
-    for (uint32_t row = 1; row < rows; row++) {
+  ::SVD::SVD<rows, columns>(A, U, sigma, Vt);
-      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,9 +5,7 @@
 #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
 
 template <uint8_t rows, uint8_t columns> class Matrix {
@@ -234,12 +232,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.
@@ -247,6 +252,24 @@ public:
   void EigenQR(Matrix<rows, rows> &eigenVectors, Matrix<rows, 1> &eigenValues,
                uint32_t maxIterations = 1000, float tolerance = 1e-6f) const;
 
+  /**
+   * @brief Compute the Singular Value Decomposition (SVD) of this matrix.
+   *
+   * Wrapper around SVD::SVD (see SVD.hpp for the full algorithm
+   * description, output storage conventions, and stack-usage notes).
+   * Decomposes A = U · Σ · Vᵀ where U is rows×columns, Σ is the vector
+   * of singular values (columns×1, sorted descending), and Vᵀ is
+   * columns×columns. Works for any shape (wide matrices are handled
+   * internally by computing SVD(Aᵀ) and swapping the factors back).
+   * This matrix is not modified.
+   *
+   * @param U Output: left singular vectors (rows×columns)
+   * @param sigma Output: singular values in descending order (columns×1)
+   * @param Vt Output: right singular vectors, transposed (columns×columns)
+   */
+  void SVD(Matrix<rows, columns> &U, Matrix<columns, 1> &sigma,
+           Matrix<columns, columns> &Vt) const;
+
 protected:
   std::array<float, rows * columns> matrix;
 

src/QR.cpp

@@ -0,0 +1,422 @@
+// 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.
+ */
+// [[maybe_unused]]: this helper is only referenced from template
+// (EigenQR/Tridiagonalize), so in translation units that include this file
+// but never instantiate those templates, the definition is legitimately
+// unused. The attribute silences -Wunused-function there without hiding
+// real dead code in TUs that do use the algorithm.
+[[maybe_unused]] 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.
+ */
+[[maybe_unused]] 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,196 @@
+#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 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 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

src/SVD.hpp

@@ -3,6 +3,21 @@
 
 /**
  * @brief library that uses Matrix.hpp and performs SVD on a matrix
+ *
+ * @note Fully templated: SVD works for ANY Matrix<R, C> with R, C in
+ *       1..255 (the uint8_t range of Matrix). There is no 5×5 limit.
+ *
+ * @note EMBEDDED CONSTRAINT — no heap. All working storage is stack
+ *       allocated as templated Matrix<N,N> buffers where
+ *       N = max(R, C). Peak stack usage per SVD call is
+ *       ≈ 11·N² floats (≈ 44·N² bytes):
+ *         N = 5   → ~1.1 KB
+ *         N = 10  → ~4.4 KB
+ *         N = 20  → ~18 KB
+ *         N = 50  → ~110 KB
+ *         N = 100 → ~440 KB
+ *         N = 255 → ~2.9 MB
+ *       Instantiate only the sizes that fit your call-stack budget.
  */
 namespace SVD {
 /**
@@ -18,14 +33,28 @@ namespace SVD {
  * The decomposition satisfies: A ≈ U × diag(Σ) × Vᵀ
  * Singular values are returned in descending order.
  *
- * @param U Output: left singular vectors (m×k orthogonal matrix)
+ * Output storage conventions:
- * @param sigma Output: singular values as k×1 column vector, sorted descending
+ *   - U: Matrix<rows, columns> — first k columns are meaningful
- * @param Vt Output: right singular vectors transposed (k×n matrix)
+ *     (rows k..columns−1 are zero in the wide case)
+ *   - sigma: Matrix<columns, 1> — first k entries are the singular
+ *     values; entries beyond k (wide matrices only) are zero
+ *   - Vt: Matrix<columns, columns> — first k rows are meaningful
+ *     (zero-padded in the tall case)
  *
- * @note This implementation uses the Golub-Kahan-Reinsch algorithm:
+ * For wide matrices (rows < columns) the SVD is computed on Aᵀ and the
- *       1. Householder bidiagonalization of A
+ * factors are swapped back.
- *       2. Implicit QR iteration on the bidiagonal matrix
+ *
- *       3. Accumulation of U and V factors throughout
+ * @tparam rows Number of rows in A (1..255)
+ * @tparam columns Number of columns in A (1..255)
+ * @param matrixToDecompose Input: the matrix A
+ * @param U Output: left singular vectors (rows×columns matrix)
+ * @param sigma Output: singular values (columns×1 vector, sorted descending)
+ * @param Vt Output: right singular vectors transposed (columns×columns)
+ *
+ * @note This implementation uses Householder bidiagonalization followed
+ *       by block reduction: 2×2 blocks via closed form, larger blocks
+ *       via cyclic Jacobi eigen-decomposition of BᵀB with residual
+ *       singular values σᵢ = ‖B·vᵢ‖ (see docs/svd-refactor.md).
  */
 template <uint8_t rows, uint8_t columns>
 void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
@@ -33,7 +62,14 @@ void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
 
 // ========================================================================
 // SVD Building Block Functions (for unit testing)
-// These operate on internal 5×5 working arrays for maximum flexibility.
+//
+// Templated on the working-buffer size N. All block operations work on
+// N×N matrices with runtime bounds (m, n, p, blockSize, ...) — the
+// regions beyond the bounds are zero-padded working space.
+//
+// N is deduced from the Matrix arguments at the call site, e.g.
+//   Matrix<8, 8> W, QL, QR;
+//   SVD::Bidiagonalize(W, 6, 8, 6, QL, QR);   // N = 8 deduced
 // ========================================================================
 
 /**
@@ -42,9 +78,9 @@ void SVD(Matrix<rows, columns> &matrixToDecompose, Matrix<rows, columns> &U,
  * Given input vector x, computes normalized v and scalar alpha such that:
  *   (I - 2·v·vᵀ) · x = [alpha, 0, 0, ...]ᵀ
  *
- * @param x Input vector (up to 5 elements)
+ * @param x Input vector (up to len elements)
  * @param len Number of valid elements in x
- * @param v Output: normalized Householder vector (v[0] is the first element)
+ * @param v Output: normalized Householder vector (length ≥ len)
  * @param alpha Output: the resulting first element after reflection
  * @return The norm of the input vector x
  */
@@ -54,30 +90,270 @@ static float ComputeHouseholder(const float *x, uint8_t len, float *v,
 /**
  * @brief Apply a Householder reflection from the left.
  *
- * Transforms W = (I - 2·v·vᵀ) · W where v operates on rows [startRow..endRow].
+ * Transforms W = (I - 2·v·vᵀ) · W where v operates on rows [startRow..endRow]
+ * and is applied across all N columns (zero-padded columns are a no-op).
  *
- * @param W Input/output: matrix to transform (5×5 working array)
+ * @tparam N Working buffer size
+ * @param W Input/output: matrix to transform
  * @param v Householder vector (length = endRow - startRow + 1)
  * @param startRow First row index
  * @param endRow Last row index
  */
-static void ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v,
+template <uint8_t N>
+static void ApplyHouseholderLeft(Matrix<N, N> &W, const float *v,
                                  uint8_t startRow, uint8_t endRow);
 
 /**
  * @brief Apply a Householder reflection from the right.
  *
  * Transforms W = W · (I - 2·v·vᵀ) where v operates on columns
- * [startCol..endCol].
+ * [startCol..endCol] and is applied across all N rows (zero-padded rows
+ * are a no-op).
  *
- * @param W Input/output: matrix to transform (5×5 working array)
+ * @tparam N Working buffer size
+ * @param W Input/output: matrix to transform
  * @param v Householder vector (length = endCol - startCol + 1)
  * @param startCol First column index
  * @param endCol Last column index
  */
-static void ApplyHouseholderRight(Matrix<5, 5> &W, const float *v,
+template <uint8_t N>
+static void ApplyHouseholderRight(Matrix<N, N> &W, const float *v,
                                   uint8_t startCol, uint8_t endCol);
 
+/**
+ * @brief Reduce a matrix to upper bidiagonal form using Householder reflections.
+ *
+ * Applies a sequence of Householder reflections to reduce the input matrix
+ * W (m×q, where q ≥ p, stored in N×N working space) to upper bidiagonal
+ * form B (p×q), accumulating the left and right transformation matrices
+ * in QL and QR respectively.
+ *
+ * Algorithm (Golub-Kahan bidiagonalization):
+ *   For k = 0 to p-1:
+ *     1. Left HH on column k, rows k..m-1: zero out subdiagonal below B[k+1][k]
+ *     2. Right HH on row k, cols k+2..q-1: zero out superdiagonal above B[k][k+1]
+ *
+ * The accumulated transformations satisfy:
+ *   QLᵀ · W_original · QR = B (upper bidiagonal)
+ *
+ * @tparam N Working buffer size (≥ m and ≥ q)
+ * @param W Input/output: matrix to bidiagonalize (first m×q used)
+ * @param m Number of rows in the working matrix
+ * @param q Number of columns in the working matrix (q ≥ p)
+ * @param p Rank = min(m, original_columns) — number of bidiagonalization steps
+ * @param QL Input/output: left Householder accumulation (initialized to identity)
+ * @param QR Input/output: right Householder accumulation (initialized to identity)
+ */
+template <uint8_t N>
+static void Bidiagonalize(Matrix<N, N> &W, uint8_t m, uint8_t q, uint8_t p,
+                          Matrix<N, N> &QL, Matrix<N, N> &QR);
+
+/**
+ * @brief Deflate a bidiagonal matrix by zeroing negligible superdiagonals.
+ *
+ * Scans the p×p upper-bidiagonal matrix stored in W and zeros out any
+ * superdiagonal element W[i][i+1] whose magnitude is negligible relative
+ * to the local diagonal scale (|W[i][i]| + |W[i+1][i+1]|). Deflating
+ * splits the matrix into independent unreduced blocks that can each be
+ * solved separately.
+ *
+ * @tparam N Working buffer size
+ * @param W Input/output: bidiagonal matrix (first p×p used)
+ * @param p Size of the bidiagonal matrix (min(rows, columns))
+ * @param tol Relative deflation tolerance (e.g. 1e-8f)
+ */
+template <uint8_t N>
+static void DeflateBidiagonal(Matrix<N, N> &W, uint8_t p, float tol);
+
+/**
+ * @brief Check whether a bidiagonal matrix has fully reduced to diagonal.
+ *
+ * Returns true when every superdiagonal element of the p×p bidiagonal
+ * matrix in W is (numerically) zero, i.e. the diagonal entries are the
+ * (unsorted) singular values and no unreduced blocks remain.
+ *
+ * @tparam N Working buffer size
+ * @param W Input: bidiagonal matrix (first p×p used)
+ * @param p Size of the bidiagonal matrix (min(rows, columns))
+ * @param tol Numerical zero threshold multiplier
+ * @return true when all superdiagonal elements are ~0
+ */
+template <uint8_t N>
+static bool BidiagonalIsDiagonal(const Matrix<N, N> &W, uint8_t p, float tol);
+
+/**
+ * @brief Compute the full SVD of a 2×2 upper-bidiagonal block (pure).
+ *
+ * Decomposes B = [[a, b], [0, d]] as:
+ *   B = Ublock · diag(sigma[0], sigma[1]) · Vblockᵀ
+ *
+ * Guarantees:
+ *   - sigma[0] ≥ sigma[1] ≥ 0 (singular values, from eigenvalues of BᵀB)
+ *   - Ublock and Vblock are orthogonal (columns are the left/right
+ *     singular vectors respectively; Vblock = scipy's Vᵀᵀ)
+ *   - Ublock · diag(sigma) · Vblockᵀ == B (within float tolerance)
+ *
+ * Math: eigenvectors of BᵀB = [[a², ab], [ab, b²+d²]] give the right
+ * singular vectors (v1 = normalize(ab, σ1²−a²) with a safe fallback when
+ * that vector is ~0; v2 = (−v1y, v1x)); left singular vectors are
+ * uᵢ = B·vᵢ/σᵢ with a rank-deficiency guard: when σᵢ ≈ 0 (i.e. ~1e-30),
+ * that U column is filled with the signed orthogonal complement of the
+ * other U column instead of dividing by ~0.
+ *
+ * @param a B[0][0] (first diagonal element)
+ * @param b B[0][1] (superdiagonal element)
+ * @param d B[1][1] (second diagonal element)
+ * @param Ublock Output: 2×2 left singular vectors (columns)
+ * @param Vblock Output: 2×2 right singular vectors (columns)
+ * @param sigma Output: singular values, sigma[0] ≥ sigma[1] ≥ 0
+ */
+static void SolveBidiagonalBlock2x2(float a, float b, float d,
+                                    float Ublock[2][2], float Vblock[2][2],
+                                    float sigma[2]);
+
+/**
+ * @brief Cyclic Jacobi eigenvalue algorithm for a symmetric matrix (pure).
+ *
+ * Reduces symmetric n×n matrix T to (near-)diagonal form IN PLACE using
+ * cyclic Jacobi rotations, accumulating the eigenvectors in V.
+ *
+ * On return:
+ *   - T's diagonal entries are the eigenvalues (off-diagonals ~0)
+ *   - evals[i] = T[i][i], UNSORTED, SIGNED (this is a general symmetric
+ *     eigen solver, not just for PSD matrices like T = BᵀB)
+ *   - columns of V are the corresponding eigenvectors (T·V = V·Λ)
+ *
+ * Convergence: relative off-diagonal tolerance 1e-10, hard-capped at
+ * 100 sweeps.
+ *
+ * @tparam N Working buffer size (≥ n)
+ * @param T Input/output: symmetric matrix (first n×n used, destroyed in place)
+ * @param n Matrix size
+ * @param evals Output: eigenvalues, unsorted, length ≥ n
+ * @param V Output: eigenvector matrix (first n×n used), columns are eigenvectors
+ */
+template <uint8_t N>
+static void JacobiEigenSymmetric(Matrix<N, N> &T, uint8_t n, float *evals,
+                                 Matrix<N, N> &V);
+
+/**
+ * @brief Fold a block SVD's factors into the QL/QR accumulators.
+ *
+ * Given the block SVD of a bidiagonal block, B = Ublock·Σ·Vblockᵀ, the
+ * accumulated Householder matrices must absorb the block factors:
+ *   QL[:, blockStart..blockStart+blockSize−1] ← QL[:, ...] · Ublock
+ *     (rows 0..rowsQL−1)
+ *   QR[:, blockStart..blockStart+blockSize−1] ← QR[:, ...] · Vblock
+ *     (rows 0..rowsQR−1)
+ *
+ * rowsQL / rowsQR are the meaningful row extents of the accumulators
+ * (e.g. for a wide matrix W = Aᵀ, QL carries n = rows(W) meaningful
+ * rows while QR is read back over its first m rows).
+ *
+ * In-place update is done through temporary buffers (updating QL's block
+ * columns while still reading them corrupts the result).
+ *
+ * @tparam N Working buffer size
+ * @param blockStart First column/row index of the block in W
+ * @param blockSize Size of the block (2, or > 2 for the Jacobi path)
+ * @param Ublock Left singular-vector factor of the block (first blockSize×blockSize used)
+ * @param Vblock Right singular-vector factor of the block (first blockSize×blockSize used)
+ * @param rowsQL Number of meaningful rows of QL
+ * @param rowsQR Number of meaningful rows of QR
+ * @param QL Input/output: left transformation accumulator
+ * @param QR Input/output: right transformation accumulator
+ */
+template <uint8_t N>
+static void ApplyBlockFactorsToAccumulators(uint8_t blockStart,
+                                            uint8_t blockSize,
+                                            const Matrix<N, N> &Ublock,
+                                            const Matrix<N, N> &Vblock,
+                                            uint8_t rowsQL, uint8_t rowsQR,
+                                            Matrix<N, N> &QL,
+                                            Matrix<N, N> &QR);
+
+/**
+ * @brief Solve a bidiagonal block larger than 2×2 via Jacobi eigen of BᵀB.
+ *
+ * Computes the full SVD of the unreduced upper-bidiagonal block
+ * W[blockStart..blockStart+blockSize−1] via eigendecomposition of the
+ * tridiagonal T = BᵀB:
+ *   1. Snapshot the ORIGINAL block diagonal/superdiagonal from W
+ *   2. Form T = BᵀB (tridiagonal symmetric)
+ *   3. JacobiEigenSymmetric on T → eigenvalues (unsorted) + V
+ *   4. Sort eigenvalues descending, reordering V columns
+ *   5. Compute RESIDUAL singular values: σᵢ = ‖B_orig · vᵢ‖
+ *      (NOT sqrt(eigenvalue) — forming BᵀB squares the condition number,
+ *       causing float noise to swamp true tiny eigenvalues for
+ *       rank-deficient blocks)
+ *   6. Re-sort σ descending, keeping V and B·v consistent
+ *   7. Build Ublock: uᵢ = B_orig · vᵢ / σᵢ (unit norm); for σᵢ ≈ 0,
+ *      use Gram-Schmidt orthogonal completion against prior U columns
+ *   8. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators
+ *   9. Write residual norms into W's diagonal and zero the block's
+ *      superdiagonals
+ *
+ * @tparam N Working buffer size (≥ blockSize)
+ * @param W Input/output: bidiagonal matrix; the block's diagonal holds
+ *        the singular values and its superdiagonals are zeroed on return
+ * @param blockStart First column/row index of the block
+ * @param blockSize Size of the block (> 2)
+ * @param rowsQL Number of meaningful rows of QL
+ * @param rowsQR Number of meaningful rows of QR
+ * @param QL Input/output: left transformation accumulator
+ * @param QR Input/output: right transformation accumulator
+ */
+template <uint8_t N>
+static void SolveBidiagonalBlockJacobi(Matrix<N, N> &W, uint8_t blockStart,
+                                       uint8_t blockSize, uint8_t rowsQL,
+                                       uint8_t rowsQR, Matrix<N, N> &QL,
+                                       Matrix<N, N> &QR);
+
+/**
+ * @brief Extract singular values from bidiagonal matrix diagonal and sort.
+ *
+ * Extracts absolute values of diagonal elements of W as singular values,
+ * then sorts them in descending order while reordering columns of QL
+ * and QR to maintain consistency. A negative diagonal element flips the
+ * sign of the corresponding QL column to keep A = U·Σ·Vᵀ.
+ *
+ * @tparam N Working buffer size
+ * @param W Input: bidiagonal matrix (first p×p used)
+ * @param sigma Output: sorted singular values (N×1 column vector, only first p used)
+ * @param p Number of singular values (min(rows, columns))
+ * @param QL Input/output: left transformation matrix (modified during sort)
+ * @param QR Input/output: right transformation matrix (modified during sort)
+ */
+template <uint8_t N>
+static void ExtractAndSortSingularValues(Matrix<N, N> &W, Matrix<N, 1> &sigma,
+                                         uint8_t p, Matrix<N, N> &QL,
+                                         Matrix<N, N> &QR);
+
+/**
+ * @brief Assemble final U and Vt matrices from QL/QR.
+ *
+ * Computes the final left singular vectors (U) and right singular vectors
+ * transposed (Vt) from the accumulated Householder transformations.
+ *
+ * For non-transpose case: U = QL[:,0:p], Vt = QR[:,0:p]ᵀ
+ * For transpose case:    U = QR[:,0:p], Vt = full QLᵀ (all n rows)
+ *
+ * @tparam N Working buffer size (≥ m and ≥ n)
+ * @param m Number of rows in original matrix
+ * @param n Number of columns in original matrix
+ * @param p Rank = min(m, n)
+ * @param transposeNeeded True if we computed SVD(Aᵀ) instead of SVD(A)
+ * @param QL Left Householder accumulation (N×N)
+ * @param QR Right Householder accumulation (N×N)
+ * @param U Output: left singular vectors (N×N, first m×p used)
+ * @param Vt Output: right singular vectors transposed (N×N, first p×n used)
+ */
+template <uint8_t N>
+static void AssembleUAndVt(uint8_t m, uint8_t n, uint8_t p,
+                           bool transposeNeeded, const Matrix<N, N> &QL,
+                           const Matrix<N, N> &QR, Matrix<N, N> &U,
+                           Matrix<N, N> &Vt);
+
 /**
  * @brief Compute a Givens rotation that zeros out y.
  *
@@ -98,6 +374,7 @@ static void ComputeGivens(float x, float y, float &c, float &s);
  *
  * Applies [c s; -s c] to rows i, j of W (columns startCol..endCol).
  *
+ * @tparam N Working buffer size
  * @param W Input/output: matrix to transform
  * @param i First row index
  * @param j Second row index
@@ -106,7 +383,8 @@ static void ComputeGivens(float x, float y, float &c, float &s);
  * @param startCol First column to transform
  * @param endCol Last column to transform
  */
-static void ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
+template <uint8_t N>
+static void ApplyGivensLeft(Matrix<N, N> &W, uint8_t i, uint8_t j, float c,
                             float s, uint8_t startCol, uint8_t endCol);
 
 /**
@@ -114,6 +392,7 @@ static void ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
  *
  * Applies [c -s; s c]ᵀ to columns i, j of W (rows startRow..endRow).
  *
+ * @tparam N Working buffer size
  * @param W Input/output: matrix to transform
  * @param i First column index
  * @param j Second column index
@@ -122,10 +401,11 @@ static void ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
  * @param startRow First row to transform
  * @param endRow Last row to transform
  */
-static void ApplyGivensRight(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
+template <uint8_t N>
+static void ApplyGivensRight(Matrix<N, N> &W, uint8_t i, uint8_t j, float c,
                              float s, uint8_t startRow, uint8_t endRow);
 } // namespace SVD
 
 #ifndef SVD_H_
 #include "SVD.cpp"
-#endif // SVD_H_
+#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
 )
 
@@ -53,3 +54,13 @@ target_link_libraries(svd-integration-test
     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
@@ -390,7 +391,7 @@ TEST_CASE("Identity Matrix", "Matrix") {
         if (oneColumnIndex == column) {
           REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
         } else {
-          REQUIRE_THAT(value, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+          REQUIRE_THAT(value, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
         }
       }
       oneColumnIndex++;
@@ -407,7 +408,7 @@ TEST_CASE("Identity Matrix", "Matrix") {
         if (oneColumnIndex == column && row < 3) {
           REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
         } else {
-          REQUIRE_THAT(value, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+          REQUIRE_THAT(value, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
         }
       }
       oneColumnIndex++;
@@ -423,7 +424,7 @@ TEST_CASE("Identity Matrix", "Matrix") {
         if (oneColumnIndex == column) {
           REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
         } else {
-          REQUIRE_THAT(value, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+          REQUIRE_THAT(value, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
         }
       }
       oneColumnIndex++;
@@ -519,7 +520,7 @@ TEST_CASE("QR Decompositions", "Matrix") {
     // check that all R values are correct
     REQUIRE_THAT(R[0][0], Catch::Matchers::WithinRel(3.16228f, 1e-4f));
     REQUIRE_THAT(R[0][1], Catch::Matchers::WithinRel(4.42719f, 1e-4f));
-    REQUIRE_THAT(R[1][0], Catch::Matchers::WithinRel(0.0f, 1e-4f));
+    REQUIRE_THAT(R[1][0], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
     REQUIRE_THAT(R[1][1], Catch::Matchers::WithinRel(0.63246f, 1e-4f));
   }
 
@@ -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_THAT(values[1][0], Catch::Matchers::WithinRel(0.0f, 1e-4f));
+    REQUIRE(componentMatches(vectors[0][1], 0.956842422f));
-    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(16.1168f, 1e-4f));
+    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::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);
+    }
   }
 }
 
@@ -753,14 +1002,14 @@ TEST_CASE("SVD: Simple 2x2 Matrix", "Matrix") {
   REQUIRE(isSortedDescending(sigma, 2));
 
   // Verify U is orthogonal: UᵀU ≈ I
-  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 
   // Verify Vt is orthogonal: VtVᵀ ≈ I
-  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 
   // Verify reconstruction: A ≈ U Σ Vᵀ
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 }
 
 TEST_CASE("SVD: Symmetric Positive Definite 2x2", "Matrix") {
@@ -777,7 +1026,7 @@ TEST_CASE("SVD: Symmetric Positive Definite 2x2", "Matrix") {
 
   // For symmetric PD matrices, U ≈ V (up to sign)
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 }
 
 TEST_CASE("SVD: Full-Rank 3x3 Matrix", "Matrix") {
@@ -797,11 +1046,11 @@ TEST_CASE("SVD: Full-Rank 3x3 Matrix", "Matrix") {
                Catch::Matchers::WithinRel(0.1968665211f, 1e-4f));
 
   REQUIRE(isSortedDescending(sigma, 3));
-  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
-  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 }
 
 TEST_CASE("SVD: Rank-Deficient 3x3 Matrix", "Matrix") {
@@ -821,12 +1070,13 @@ TEST_CASE("SVD: Rank-Deficient 3x3 Matrix", "Matrix") {
   REQUIRE(sigma.Get(2, 0) < 1e-3f);
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 }
 
 TEST_CASE("SVD: Diagonal 3x3 Matrix", "Matrix") {
   // For a diagonal matrix, σ = diagonal entries, U = V = I
-  Matrix<3, 3> A{10.0f, 0.0f, 0.0f, 5.0f, 0.0f, 0.0f, 0.0f, 0.0f, 2.0f};
+  // Row-major init: [10,0,0, 0,5,0, 0,0,2] = diag(10,5,2)
+  Matrix<3, 3> A{10.0f, 0.0f, 0.0f, 0.0f, 5.0f, 0.0f, 0.0f, 0.0f, 2.0f};
   Matrix<3, 3> U{}, Vt{};
   Matrix<3, 1> sigma{};
 
@@ -837,7 +1087,7 @@ TEST_CASE("SVD: Diagonal 3x3 Matrix", "Matrix") {
   REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.0f, 1e-4f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 }
 
 TEST_CASE("SVD: Tall Matrix (4×3)", "Matrix") {
@@ -858,10 +1108,10 @@ TEST_CASE("SVD: Tall Matrix (4×3)", "Matrix") {
   REQUIRE(sigma.Get(2, 0) < 1e-3f);
 
   // U should be 4×3 with orthonormal columns
-  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 }
 
 TEST_CASE("SVD: Wide Matrix (3×5)", "Matrix") {
@@ -882,10 +1132,10 @@ TEST_CASE("SVD: Wide Matrix (3×5)", "Matrix") {
   REQUIRE(sigma.Get(2, 0) < 1e-3f);
 
   // Vt should be 5×5 with orthonormal rows (first k)
-  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 }
 
 TEST_CASE("SVD: 5×5 Symmetric Tridiagonal", "Matrix") {
@@ -908,11 +1158,137 @@ TEST_CASE("SVD: 5×5 Symmetric Tridiagonal", "Matrix") {
                Catch::Matchers::WithinRel(0.2679491924f, 1e-4f));
 
   REQUIRE(isSortedDescending(sigma, 5));
-  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
-  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("SVD: 5×5 Full-Rank Random", "Matrix") {
+  // Reference: scipy.linalg.svd, np.random.default_rng(7).standard_normal((5,5))
+  // cond ≈ 11.5
+  // σ = [3.04651784, 2.22681732, 1.84290662, 1.02101969, 0.264826749]
+  Matrix<5, 5> A{0.00123015f, 0.298746f, -0.274138f, -0.890592f, -0.454671f,
+                 -0.991647f, 0.0601436f, 1.34022f, -0.492207f, -0.620475f,
+                 0.489842f, 0.356887f, 0.105414f, -0.930468f, -0.0292518f,
+                 0.695303f, -1.34421f, -0.457616f, -1.90122f, -1.28954f,
+                 -1.84174f, -0.235091f, -1.26745f, 0.271264f, 0.156751f};
+  Matrix<5, 5> U{}, Vt{};
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(3.04651784f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.22681732f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.84290662f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0),
+               Catch::Matchers::WithinRel(1.02101969f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(0.264826749f, 1e-4f));
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("SVD: 5×5 Rank-Deficient (rank 3)", "Matrix") {
+  // Reference: scipy.linalg.svd of rng.standard_normal((5,3)) @
+  // rng.standard_normal((3,5)) — exactly rank 3
+  // σ = [7.29829771, 2.76487864, 1.57392325, ~1e-16, ~1e-16]
+  Matrix<5, 5> A{3.46252f, -1.52873f, -0.111526f, 1.28954f, -5.18688f,
+                 -1.34702f, 1.92936f, -0.0410797f, -0.958791f, 0.449623f,
+                 0.844592f, 0.0986352f, 0.408213f, 0.124867f, -2.45393f,
+                 1.34177f, -0.587312f, -1.39847f, 0.580032f, -0.167692f,
+                 0.915462f, -0.311165f, -1.16141f, 0.377283f, 0.055373f};
+  Matrix<5, 5> U{}, Vt{};
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(7.29829771f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.76487864f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.57392325f, 1e-4f));
+  // The two rank-deficient singular values must be at noise level
+  REQUIRE(sigma.Get(3, 0) < 1e-3f);
+  REQUIRE(sigma.Get(4, 0) < 1e-3f);
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  // Rank-3 matrix: the top-3 SVD terms must reproduce A
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 5e-3f));
+}
+
+TEST_CASE("SVD: 5×5 Wide Dynamic Range (cond ≈ 9000)", "Matrix") {
+  // Symmetric banded, diagonal decays 50 → 1e-3, off-diagonals 3 → 0.01
+  // Reference: scipy.linalg.svd
+  // σ = [50.1496395, 10.1719501, 1.02846061, 0.0207253626, 0.00557535757]
+  Matrix<5, 5> A{50.0f, 3.0f, 0.0f, 0.0f, 0.0f,
+                 -3.0f, 10.0f, 0.5f, 0.0f, 0.0f,
+                 0.0f, -0.5f, 1.0f, 0.08f, 0.0f,
+                 0.0f, 0.0f, -0.08f, 0.01f, 0.01f,
+                 0.0f, 0.0f, 0.0f, -0.01f, 0.001f};
+  Matrix<5, 5> U{}, Vt{};
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(50.1496395f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(10.1719501f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0),
+               Catch::Matchers::WithinRel(1.02846061f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0),
+               Catch::Matchers::WithinRel(0.0207253626f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(4, 0),
+               Catch::Matchers::WithinRel(0.00557535757f, 1e-3f));
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  // Relative to the largest entry (‖A‖F ≈ 50.1)
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 5e-2f));
+}
+
+TEST_CASE("SVD: 5×5 Symmetric Indefinite", "Matrix") {
+  // Symmetric with negative eigenvalues — σ must equal |eigenvalues|
+  // Reference: scipy.linalg.svd
+  // σ = [4.70141723, 3.76392521, 2.73426097, 1.15773083, 0.642665756]
+  Matrix<5, 5> A{2.0f, -1.0f, 0.0f, 0.0f, 0.5f,
+                 -1.0f, 2.0f, -1.0f, 0.0f, 0.0f,
+                 0.0f, -1.0f, 3.0f, -1.0f, 0.0f,
+                 0.0f, 0.0f, -1.0f, 2.0f, -1.0f,
+                 0.5f, 0.0f, 0.0f, -1.0f, 4.0f};
+  Matrix<5, 5> U{}, Vt{};
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(4.70141723f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(3.76392521f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0),
+               Catch::Matchers::WithinRel(2.73426097f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0),
+               Catch::Matchers::WithinRel(1.15773083f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0),
+               Catch::Matchers::WithinRel(0.642665756f, 1e-4f));
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 }
 
 TEST_CASE("SVD: Non-Square with Negative Values (2×3)", "Matrix") {
@@ -931,7 +1307,7 @@ TEST_CASE("SVD: Non-Square with Negative Values (2×3)", "Matrix") {
                Catch::Matchers::WithinRel(0.6646227432f, 1e-4f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 }
 
 TEST_CASE("SVD: Near-Singular 2×2 Matrix", "Matrix") {
@@ -948,7 +1324,7 @@ TEST_CASE("SVD: Near-Singular 2×2 Matrix", "Matrix") {
   REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1e-6f, 1e-2f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
 }
 
 TEST_CASE("SVD: Orthogonal Matrix (3×3)", "Matrix") {
@@ -968,7 +1344,7 @@ TEST_CASE("SVD: Orthogonal Matrix (3×3)", "Matrix") {
   REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 }
 
 TEST_CASE("SVD: Identity Matrix", "Matrix") {
@@ -984,7 +1360,7 @@ TEST_CASE("SVD: Identity Matrix", "Matrix") {
   REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
 }
 
 TEST_CASE("SVD: Zero Matrix", "Matrix") {
@@ -1000,7 +1376,7 @@ TEST_CASE("SVD: Zero Matrix", "Matrix") {
   REQUIRE(sigma.Get(2, 0) < 1e-6f);
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
 }
 
 TEST_CASE("SVD: 2×1 Column Vector", "Matrix") {
@@ -1016,7 +1392,7 @@ TEST_CASE("SVD: 2×1 Column Vector", "Matrix") {
   REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.0f, 1e-4f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
 }
 
 TEST_CASE("SVD: 1×2 Row Vector", "Matrix") {
@@ -1032,5 +1408,142 @@ TEST_CASE("SVD: 1×2 Row Vector", "Matrix") {
   REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.0f, 1e-4f));
 
   float reconErr = svdReconstructionError(A, U, sigma, Vt);
-  REQUIRE_THAT(reconErr, Catch::Matchers::WithinRel(0.0f, 1e-4f));
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+}
+// ============================================================================
+// SVD Tests — Large-Size Instantiations (N > 5)
+//
+// The SVD is templated on N = max(rows, cols) with stack-only buffers, so
+// these cases exercise instantiations beyond the old 5×5 hard limit:
+// 7×5 (N=7, tall), 6×6 (N=6, square), 5×8 (N=8, wide/transpose path),
+// 6×4 (N=6, tall, near rank-deficiency → deflation path).
+// Reference singular values: scipy.linalg.svd.
+// ============================================================================
+
+TEST_CASE("SVD: Tall 7×5 Matrix (N=7)", "Matrix") {
+  // Reference: scipy.linalg.svd
+  // σ = [7.9180769443, 4.6593687008, 4.2921645616, 2.6009010840, 1.9842770351]
+  Matrix<7, 5> A{-0.7528f, 2.7043f, 1.392f, 0.592f, -2.0639f,
+                 -2.064f, -2.6515f, 2.1971f, 0.6067f, 1.2484f,
+                 -2.8765f, 2.8195f, 1.9947f, -1.726f, -1.9091f,
+                 -1.8996f, -1.1745f, 0.1485f, -0.4083f, -1.2526f,
+                 0.6711f, -2.163f, -1.2471f, -0.8018f, -0.2636f,
+                 1.7111f, -1.802f, 0.0854f, 0.5545f, -2.7213f,
+                 0.6453f, -1.9769f, -2.6097f, 2.6933f, 2.7938f};
+  Matrix<7, 5> U{};
+  Matrix<5, 5> Vt{};
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(7.9180769443f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(4.6593687008f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(4.2921645616f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(2.6009010840f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(1.9842770351f, 1e-4f));
+
+  REQUIRE(isSortedDescending(sigma, 5));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("SVD: Square 6×6 Matrix (N=6)", "Matrix") {
+  // Reference: scipy.linalg.svd (float32 inputs)
+  // σ = [5.018912792, 4.244967461, 2.505512476,
+  //      1.838801861, 0.9111995101, 0.4580149353]
+  Matrix<6, 6> A{1.2336f, -0.7815f, -1.6093f, 0.7369f, -0.2394f, -1.5118f,
+                 -0.0193f, -1.8624f, 1.6373f, -0.9649f, 0.6501f, -0.7532f,
+                 0.0803f, 0.1868f, -1.2606f, 1.8783f,
+                 1.1005f, 1.758f, 1.5793f, 0.3916f, 1.6875f, -1.646f,
+                 -1.2161f, -1.8191f, -0.6987f, -0.4453f, -0.9146f, 1.315f,
+                 -0.573f, -0.8763f, 0.1708f, -1.4363f, 1.2088f, -1.7018f,
+                 1.089f, 1.9475f};
+  Matrix<6, 6> U{}, Vt{};
+  Matrix<6, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.018912792f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(4.244967461f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.505512476f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(1.838801861f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(0.9111995101f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(5, 0), Catch::Matchers::WithinRel(0.4580149353f, 1e-4f));
+
+  REQUIRE(isSortedDescending(sigma, 6));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("SVD: Wide 5×8 Matrix (N=8, transpose path)", "Matrix") {
+  // Reference: scipy.linalg.svd
+  // σ = [5.8027782929, 4.1105282764, 3.7755966048, 3.3208483982, 2.0321410547]
+  //
+  // Wide matrices take the Aᵀ transpose path; Vt must be the FULL 8×8
+  // orthogonal matrix (all 8 rows meaningful), not just the top 5.
+  Matrix<5, 8> A{-1.5064f, -2.4724f, 1.5773f, 1.0343f, 1.145f, 1.3564f, -2.1298f, -0.7077f,
+                 -1.9207f, 1.8155f, 0.6165f, -0.8455f, -2.1822f, -0.9451f, -0.8741f, 1.148f,
+                 0.6878f, 1.9361f, -0.1389f, -1.902f, 1.0662f, 1.3039f, 0.3064f, 1.3548f,
+                 -0.031f, 0.1137f, -0.3623f, -2.3729f, -1.9605f, -2.3429f, 0.6821f, -0.9282f,
+                 0.0429f, 2.0378f, -1.2535f, -0.4481f, 1.2778f, -1.356f, -2.1151f, -1.0512f};
+  Matrix<5, 8> U{};
+  Matrix<8, 8> Vt{};
+  Matrix<8, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.8027782929f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(4.1105282764f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(3.7755966048f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(3.3208483982f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0), Catch::Matchers::WithinRel(2.0321410547f, 1e-4f));
+  // Remaining singular values must be at noise level
+  REQUIRE(sigma.Get(5, 0) < 1e-3f);
+  REQUIRE(sigma.Get(6, 0) < 1e-3f);
+  REQUIRE(sigma.Get(7, 0) < 1e-3f);
+
+  REQUIRE(isSortedDescending(sigma, 8));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("SVD: Tall 6×4 Near Rank-Deficient (N=6, deflation path)", "Matrix") {
+  // Reference: scipy.linalg.svd
+  // σ = [5.9434060901, 3.2857910666, 0.3066158795, 6.48e-07]
+  //
+  // σ₄ ≈ 6.5e-7 forces the deflation logic to zero the last
+  // superdiagonal and isolate the trailing 1×1 block.
+  Matrix<6, 4> A{-0.086904f, 1.410225f, 1.308323f, 2.234762f,
+                 0.022123f, 0.896751f, 0.324176f, 0.773607f,
+                 -0.473015f, 1.555111f, 0.290059f, 1.157726f,
+                 -0.78371f, 1.398884f, -1.930606f, -1.548717f,
+                 0.201518f, -0.626835f, 0.976596f, 0.875294f,
+                 -1.24206f, 1.60595f, -3.078089f, -2.73695f};
+  Matrix<6, 4> U{};
+  Matrix<4, 4> Vt{};
+  Matrix<4, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.9434060901f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(3.2857910666f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(0.3066158795f, 1e-4f));
+  // Fourth singular value is at noise level (matrix is ~rank 3)
+  REQUIRE(sigma.Get(3, 0) < 1e-4f);
+
+  REQUIRE(isSortedDescending(sigma, 4));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+  REQUIRE_THAT(orthogonalityError(Vt), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 }

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()

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

@@ -783,3 +783,865 @@ TEST_CASE("SVD Building Block: Givens preserves orthogonality",
     REQUIRE(isOrthogonal5(M_right));
   }
 }
+
+// ============================================================================
+// TEST 11: ExtractAndSortSingularValues (Phase 3)
+// ===========================================================================
+TEST_CASE("SVD Phase 3: ExtractAndSortSingularValues", "[Matrix][SVD]") {
+  // Test case 1: Diagonal matrix with unordered singular values
+  {
+    Matrix<5, 5> W{0};
+    W[0][0] = 2.0f;
+    W[1][1] = 10.0f;
+    W[2][2] = 5.0f;
+
+    Matrix<5, 1> sigma{0};
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::ExtractAndSortSingularValues(W, sigma, 3, QL, QR);
+
+    // Singular values should be sorted descending: [10, 5, 2]
+    REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(10.0f, 1e-6f));
+    REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(5.0f, 1e-6f));
+    REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.0f, 1e-6f));
+
+    // QL and QR columns should have been swapped to match the sort order
+    // Original: col 0 → σ=2, col 1 → σ=10, col 2 → σ=5
+    // After sort: col 0 has σ=10 (was orig col 1), col 1 has σ=5 (was orig col 2),
+    //             col 2 has σ=2 (was orig col 0)
+    // Starting from identity: QL[:,0] = e₀, QL[:,1] = e₁, QL[:,2] = e₂
+    // After swaps: QL[:,0] = e₁, QL[:,1] = e₂, QL[:,2] = e₀
+    REQUIRE_THAT(QL.Get(0, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(1, 0), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(1, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(2, 1), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(0, 2), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(1, 2), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(QL.Get(2, 2), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+  }
+
+  // Test case 2: Negative diagonal elements (absolute value extraction)
+  {
+    Matrix<5, 5> W{0};
+    W[0][0] = -3.0f;
+    W[1][1] = -7.0f;
+    W[2][2] = 5.0f;
+
+    Matrix<5, 1> sigma{0};
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::ExtractAndSortSingularValues(W, sigma, 3, QL, QR);
+
+    // Should extract absolute values and sort: [7, 5, 3]
+    REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(7.0f, 1e-6f));
+    REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(5.0f, 1e-6f));
+    REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(3.0f, 1e-6f));
+  }
+
+  // Test case 3: Already sorted (no swaps needed)
+  {
+    Matrix<5, 5> W{0};
+    W[0][0] = 9.0f;
+    W[1][1] = 6.0f;
+    W[2][2] = 3.0f;
+
+    Matrix<5, 1> sigma{0};
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::ExtractAndSortSingularValues(W, sigma, 3, QL, QR);
+
+    REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(9.0f, 1e-6f));
+    REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(6.0f, 1e-6f));
+    REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(3.0f, 1e-6f));
+
+    // QL and QR should be unchanged (identity)
+    for (uint8_t i = 0; i < 5; i++) {
+      for (uint8_t j = 0; j < 5; j++) {
+        float expected = (i == j) ? 1.0f : 0.0f;
+        REQUIRE_THAT(QL.Get(i, j), Catch::Matchers::WithinRel(expected, 1e-6f));
+        REQUIRE_THAT(QR.Get(i, j), Catch::Matchers::WithinRel(expected, 1e-6f));
+      }
+    }
+  }
+
+  // Test case 4: Single singular value
+  {
+    Matrix<5, 5> W{0};
+    W[0][0] = 42.0f;
+
+    Matrix<5, 1> sigma{0};
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::ExtractAndSortSingularValues(W, sigma, 1, QL, QR);
+
+    REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(42.0f, 1e-6f));
+  }
+}
+
+// ============================================================================
+// TEST 12: AssembleUAndVt (Phase 4)
+// ===========================================================================
+TEST_CASE("SVD Phase 4: AssembleUAndVt", "[Matrix][SVD]") {
+  // Test case 1: Non-transpose case (m ≥ n) — U from QL, Vt from QRᵀ
+  {
+    uint8_t m = 3, n = 2, p = 2;
+    bool transposeNeeded = false;
+
+    Matrix<5, 5> QL{0};
+    // Make columns orthonormal
+    QL[0][0] = 3.0f / 5.0f;
+    QL[1][0] = 4.0f / 5.0f;
+    QL[2][0] = 0.0f;
+    QL[0][1] = 4.0f / 5.0f;
+    QL[1][1] = -3.0f / 5.0f;
+    QL[2][1] = 0.0f;
+
+    Matrix<5, 5> QR{0};
+    QR[0][0] = 1.0f;   // Vt[:,0]ᵀ
+    QR[1][0] = 0.0f;
+    QR[0][1] = 0.0f;   // Vt[:,1]ᵀ
+    QR[1][1] = 1.0f;
+
+    Matrix<5, 5> U{0};
+    Matrix<5, 5> Vt{0};
+
+    SVD::AssembleUAndVt(m, n, p, transposeNeeded, QL, QR, U, Vt);
+
+    // U should be QL[:,0:2]
+    REQUIRE_THAT(U.Get(0, 0), Catch::Matchers::WithinRel(3.0f / 5.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(1, 0), Catch::Matchers::WithinRel(4.0f / 5.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(0, 1), Catch::Matchers::WithinRel(4.0f / 5.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(1, 1), Catch::Matchers::WithinRel(-3.0f / 5.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(2, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+
+    // Vt should be QR[:,0:2]ᵀ
+    REQUIRE_THAT(Vt.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(1, 1), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+
+    // Verify U is orthogonal (first p columns)
+    Matrix<5, 5> Ut = U.Transpose();
+    Matrix<5, 5> UtU{0};
+    Ut.Mult(U, UtU);
+    REQUIRE_THAT(UtU.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(UtU.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(UtU.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(UtU.Get(1, 1), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+  }
+
+  // Test case 2: Transpose case (m < n) — U from QRᵀ, Vt from QLᵀ
+  {
+    uint8_t m = 2, n = 3, p = 2;
+    bool transposeNeeded = true;
+
+    Matrix<5, 5> QL{0};
+    QL[0][0] = 1.0f;   // Vt[:,0]ᵀ
+    QL[1][0] = 0.0f;
+    QL[0][1] = 0.0f;   // Vt[:,1]ᵀ
+    QL[1][1] = 1.0f;
+
+    Matrix<5, 5> QR{0};
+    QR[0][0] = 3.0f / 5.0f;   // U[:,0]
+    QR[1][0] = 4.0f / 5.0f;
+    QR[2][0] = 0.0f;
+    QR[0][1] = 4.0f / 5.0f;
+    QR[1][1] = -3.0f / 5.0f;
+    QR[2][1] = 0.0f;
+
+    Matrix<5, 5> U{0};
+    Matrix<5, 5> Vt{0};
+
+    SVD::AssembleUAndVt(m, n, p, transposeNeeded, QL, QR, U, Vt);
+
+    // U should be QR[:,0:2]ᵀ → U[i][j] = QR[j][i]
+    REQUIRE_THAT(U.Get(0, 0), Catch::Matchers::WithinRel(3.0f / 5.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(0, 1), Catch::Matchers::WithinRel(4.0f / 5.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(1, 0), Catch::Matchers::WithinRel(4.0f / 5.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(1, 1), Catch::Matchers::WithinRel(-3.0f / 5.0f, 1e-6f));
+
+    // Vt should be QL[:,0:2]ᵀ → Vt[i][j] = QL[j][i]
+    REQUIRE_THAT(Vt.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(1, 1), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+
+    // Verify U has correct dimensions (m×n = 2×3)
+    REQUIRE(U.Get(0, 2) == 0.0f);
+    REQUIRE(U.Get(1, 2) == 0.0f);
+
+    // Verify Vt is orthogonal (first p rows)
+    Matrix<5, 5> VtVtT{0};
+    Vt.Mult(Vt.Transpose(), VtVtT);
+    REQUIRE_THAT(VtVtT.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(VtVtT.Get(1, 1), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    // Row 2 of Vt is all zeros (p=2 < n=3), so VtVtT[2][2] = 0 is expected
+  }
+
+  // Test case 3: Square matrix (m = n)
+  {
+    uint8_t m = 2, n = 2, p = 2;
+    bool transposeNeeded = false;
+
+    Matrix<5, 5> QL{0};
+    QL[0][0] = 1.0f; QL[1][1] = 1.0f;
+
+    Matrix<5, 5> QR{0};
+    QR[0][0] = 0.6f; QR[0][1] = 0.8f;
+    QR[1][0] = 0.8f; QR[1][1] = -0.6f;
+
+    Matrix<5, 5> U{0};
+    Matrix<5, 5> Vt{0};
+
+    SVD::AssembleUAndVt(m, n, p, transposeNeeded, QL, QR, U, Vt);
+
+    // U = QL[:,0:2]
+    REQUIRE_THAT(U.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+    REQUIRE_THAT(U.Get(1, 1), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+
+    // Vt = QR[:,0:2]ᵀ
+    REQUIRE_THAT(Vt.Get(0, 0), Catch::Matchers::WithinRel(0.6f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(0, 1), Catch::Matchers::WithinRel(0.8f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(1, 0), Catch::Matchers::WithinRel(0.8f, 1e-6f));
+    REQUIRE_THAT(Vt.Get(1, 1), Catch::Matchers::WithinRel(-0.6f, 1e-6f));
+  }
+}
+
+// ============================================================================
+// TEST 13: Bidiagonalize (Phase 1) - Square matrix
+// ===========================================================================
+TEST_CASE("SVD Phase 1: Bidiagonalize square matrix", "[Matrix][SVD]") {
+  // Test case 1: 3×3 matrix
+  // C++ verified reference:
+  //   W[0] = [-4.123106, -5.335784,  6.548462]
+  //   W[1] = [ 0.000000,  7.037714, -8.107580]
+  //   W[2] = [ 0.000000,  0.000000,  0.620321]
+  // Note: W[0][2]=6.548462 is NOT zeroed because right HH at k=0 has only 1 element
+  {
+    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0.0f, 0.0f,
+                   4.0f, 5.0f, 6.0f, 0.0f, 0.0f,
+                   0.0f, 7.0f, 8.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 3, 3, 3, QL, QR);
+
+    // Subdiagonal elements should be zero: W[1][0], W[2][0], W[2][1]
+    REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(2, 1), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+
+    // Verify orthogonality
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+
+  // Test case 2: 2×2 matrix (simplest non-trivial case)
+  // C++ verified reference:
+  //   W[0] = [-3.162278, -4.427189]
+  //   W[1] = [-0.000000,  0.632456]
+  {
+    Matrix<5, 5> W{3.0f, 4.0f, 0.0f, 0.0f, 0.0f,
+                   1.0f, 2.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 2, 2, 2, QL, QR);
+
+    // For 2×2, bidiagonal form has no elements to zero out
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+
+  // Test case 3: Diagonal matrix (no transformations needed)
+  // C++ verified reference: unchanged
+  {
+    Matrix<5, 5> W{10.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 5.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 2.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> W_orig{10.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                        0.0f, 5.0f, 0.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 2.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 3, 3, 3, QL, QR);
+
+    // Diagonal matrix may have sign flips but absolute values preserved
+    float err = 0.0f;
+    for (uint8_t i = 0; i < 3; i++) {
+        float diff = fabsf(W.Get(i, i)) - fabsf(W_orig.Get(i, i));
+        err += diff * diff;
+    }
+    REQUIRE_THAT(sqrtf(err), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+
+    // QL and QR may have sign flips but should remain orthogonal
+    // Check that |QL[i][j]| and |QR[i][j]| match identity pattern
+    for (uint8_t i = 0; i < 5; i++) {
+      for (uint8_t j = 0; j < 5; j++) {
+        float expected = (i == j) ? 1.0f : 0.0f;
+        REQUIRE_THAT(fabsf(QL.Get(i, j)), Catch::Matchers::WithinAbs(expected, 1e-6f));
+        REQUIRE_THAT(fabsf(QR.Get(i, j)), Catch::Matchers::WithinAbs(expected, 1e-6f));
+      }
+    }
+  }
+}
+
+// ============================================================================
+// TEST 14: Bidiagonalize (Phase 1) — Tall matrix (m > n)
+// ===========================================================================
+TEST_CASE("SVD Phase 1: Bidiagonalize tall matrix", "[Matrix][SVD]") {
+  // Test case 1: 4×3 matrix
+  // C++ verified reference (partial - subdiagonal zeros):
+  //   W[0] = [-4.123106, -5.335784,  6.548462]
+  //   W[1] = [-0.000000,  12.228222, 12.026182]
+  //   W[2] = [ 0.000000,  0.000000, -1.577527]
+  //   W[3] = [ 0.000000,  0.000000,  0.000000]
+  {
+    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0.0f, 0.0f,
+                   4.0f, 5.0f, 6.0f, 0.0f, 0.0f,
+                   0.0f, 7.0f, 8.0f, 0.0f, 0.0f,
+                   0.0f, 10.0f, 9.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 4, 3, 3, QL, QR);
+
+    // Zero below subdiagonal: W[2][0], W[3][0], W[3][1]
+    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(3, 1), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+
+    // Verify orthogonality
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+
+  // Test case 2: 3×2 matrix
+  // C++ verified reference:
+  //   W[0] = [-5.916080, -7.437357]
+  //   W[1] = [-0.000001,  0.828077]
+  //   W[2] = [-0.000000, -0.000000]
+  {
+    Matrix<5, 5> W{1.0f, 2.0f, 0.0f, 0.0f, 0.0f,
+                   3.0f, 4.0f, 0.0f, 0.0f, 0.0f,
+                   5.0f, 6.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 3, 2, 2, QL, QR);
+
+    // Zero below subdiagonal: W[2][0] ≈ 0
+    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+
+    // Verify orthogonality
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+
+  // Test case 3: 5×3 matrix (full 5-row tall)
+  {
+    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0.0f, 0.0f,
+                   4.0f, 5.0f, 6.0f, 0.0f, 0.0f,
+                   0.0f, 7.0f, 8.0f, 0.0f, 0.0f,
+                   0.0f, 10.0f, 9.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 5, 3, 3, QL, QR);
+
+    // Zero below subdiagonal: W[2][0], W[3][0], W[4][0], W[3][1], W[4][1]
+    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(4, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(3, 1), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(4, 1), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+
+    // Verify orthogonality
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+}
+
+// ============================================================================
+// TEST 15: Bidiagonalize (Phase 1) — Wide matrix (m < n)
+// ===========================================================================
+TEST_CASE("SVD Phase 1: Bidiagonalize wide matrix", "[Matrix][SVD]") {
+  // Test case 1: 2×4 matrix
+  // C++ verified reference:
+  //   W[0] = [-5.099020, -6.275717, 11.401754,  0.000000]
+  //   W[1] = [ 0.000000, -0.784465,  2.806586, -0.350823]
+  // Note: W[1][2]=2.806586 and W[1][3]=-0.350823 are NOT zeroed because
+  // right HH at k=0 has only 2 elements (cols 2,3), so it zeros col 3 but preserves col 2
+  {
+    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 4.0f, 0.0f,
+                   5.0f, 6.0f, 7.0f, 8.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 2, 4, 2, QL, QR);
+
+    // For 2×4: right HH at k=0 has 2 elements (cols 2,3)
+    // It zeros col 3 but preserves col 2 as the superdiagonal element for row 1
+    // W[0][3] should be zeroed (above superdiagonal in row 0)
+    REQUIRE_THAT(W.Get(0, 3), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    // W[1][2] and W[1][3] are part of the bidiagonal structure for row 1
+    // (superdiagonal at col 2, and right HH preserves first element)
+    REQUIRE_THAT(W.Get(1, 2), !Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+
+    // Verify orthogonality
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+
+  // Test case 2: 3×5 matrix
+  {
+    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 4.0f, 5.0f,
+                   6.0f, 7.0f, 8.0f, 9.0f, 10.0f,
+                   0.0f, 11.0f, 12.0f, 13.0f, 14.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 3, 5, 3, QL, QR);
+
+    // For 3×5: check that subdiagonal elements are zero
+    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(2, 1), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+
+    // Verify orthogonality
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+
+  // Test case 3: 1×3 matrix (row vector)
+  // C++ verified reference: W[0] = [1.0, 2.0, 3.0] (no transformations needed)
+  {
+    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 1, 3, 1, QL, QR);
+
+    // For 1×3, no transformations needed
+    REQUIRE_THAT(W.Get(0, 0), Catch::Matchers::WithinAbs(1.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(0, 1), Catch::Matchers::WithinAbs(2.0f, 1e-4f));
+    REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(3.0f, 1e-4f));
+
+    // Verify orthogonality
+    REQUIRE(isOrthogonal5(QL));
+    REQUIRE(isOrthogonal5(QR));
+  }
+}
+
+// ============================================================================
+// TEST 16: Bidiagonalize — Reconstruction property
+// ===========================================================================
+TEST_CASE("SVD Phase 1: Bidiagonalize reconstruction property", "[Matrix][SVD]") {
+  // Test: QLᵀ · W_original · QR = B (bidiagonal)
+  // This verifies that the accumulated transformations correctly represent
+  // the bidiagonalization.
+  {
+    Matrix<5, 5> W_orig{1.0f, 2.0f, 3.0f, 0.0f, 0.0f,
+                        4.0f, 5.0f, 6.0f, 0.0f, 0.0f,
+                        0.0f, 7.0f, 8.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> W = W_orig;
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 3, 3, 3, QL, QR);
+
+    // Compute QLᵀ · W_orig · QR and verify it equals W (the bidiagonal result)
+    Matrix<5, 5> Qt = QL.Transpose();
+    Matrix<5, 5> QtW_orig{0};
+    Qt.Mult(W_orig, QtW_orig);
+
+    Matrix<5, 5> QtW_origQR{0};
+    QtW_orig.Mult(QR, QtW_origQR);
+
+    // The reconstruction should match the bidiagonal result
+    float err = frobeniusNorm5(W - QtW_origQR);
+    REQUIRE_THAT(err, Catch::Matchers::WithinAbs(1e-3f, 1e-3f));
+  }
+
+  // Test: Tall matrix reconstruction (4×3)
+  {
+    Matrix<5, 5> W_orig{1.0f, 2.0f, 3.0f, 0.0f, 0.0f,
+                        4.0f, 5.0f, 6.0f, 0.0f, 0.0f,
+                        0.0f, 7.0f, 8.0f, 0.0f, 0.0f,
+                        0.0f, 10.0f, 9.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> W = W_orig;
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 4, 3, 3, QL, QR);
+
+    Matrix<5, 5> Qt = QL.Transpose();
+    Matrix<5, 5> QtW_orig{0};
+    Qt.Mult(W_orig, QtW_orig);
+
+    Matrix<5, 5> QtW_origQR{0};
+    QtW_orig.Mult(QR, QtW_origQR);
+
+    float err = frobeniusNorm5(W - QtW_origQR);
+    REQUIRE_THAT(err, Catch::Matchers::WithinAbs(1e-3f, 1e-3f));
+  }
+
+  // Test: Wide matrix reconstruction (2×4)
+  {
+    Matrix<5, 5> W_orig{1.0f, 2.0f, 3.0f, 4.0f, 0.0f,
+                        5.0f, 6.0f, 7.0f, 8.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                        0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+    Matrix<5, 5> W = W_orig;
+
+    Matrix<5, 5> QL{0}, QR{0};
+    for (uint8_t i = 0; i < 5; i++) {
+      QL[i][i] = 1.0f;
+      QR[i][i] = 1.0f;
+    }
+
+    SVD::Bidiagonalize(W, 2, 4, 2, QL, QR);
+
+    Matrix<5, 5> Qt = QL.Transpose();
+    Matrix<5, 5> QtW_orig{0};
+    Qt.Mult(W_orig, QtW_orig);
+
+    Matrix<5, 5> QtW_origQR{0};
+    QtW_orig.Mult(QR, QtW_origQR);
+
+    float err = frobeniusNorm5(W - QtW_origQR);
+    REQUIRE_THAT(err, Catch::Matchers::WithinAbs(1e-3f, 1e-3f));
+  }
+}
+
+// ============================================================================
+// TEST: SolveBidiagonalBlock2x2 — 2×2 upper-bidiagonal block SVD
+// ============================================================================
+// Reference singular values generated with scipy.linalg.svd for
+// B = [[a, b], [0, d]].
+TEST_CASE("SVD Building Block: SolveBidiagonalBlock2x2", "[Matrix][SVD]") {
+  struct Case2x2 {
+    float a, b, d;
+    float refSigma[2];
+  };
+  const Case2x2 cases[] = {
+      {2.5f, -1.3f, 0.8f, {2.84346151f, 0.70336806f}},
+      {3.0f, 0.0f, 1.0f, {3.0f, 1.0f}},
+      {1.0f, 2.0f, 0.0f, {2.23606798f, 0.0f}},
+      {-1.5f, 0.7f, -2.2f, {2.37779179f, 1.38784229f}},
+      {1.0f, 1e-4f, 0.0f, {1.0f, 0.0f}},
+      {-1.770486f, 0.281880f, 0.208573f, {1.79308863f, 0.20594385f}},
+      {0.866025f, 1.0f, 0.5f, {1.37890797f, 0.31402567f}},
+  };
+
+  for (const auto &tc : cases) {
+    float Ublock[2][2] = {{0}}, Vblock[2][2] = {{0}}, sigma[2] = {0};
+    SVD::SolveBidiagonalBlock2x2(tc.a, tc.b, tc.d, Ublock, Vblock, sigma);
+
+    // 1. Singular values match scipy
+    REQUIRE_THAT(sigma[0],
+                 Catch::Matchers::WithinRel(tc.refSigma[0], 1e-3f));
+    if (tc.refSigma[1] > 0.0f) {
+      REQUIRE_THAT(sigma[1],
+                   Catch::Matchers::WithinRel(tc.refSigma[1], 1e-3f));
+    } else {
+      REQUIRE(sigma[1] < 1e-3f);
+    }
+    REQUIRE(sigma[0] >= sigma[1]);
+
+    // 2. Ublock and Vblock are orthogonal (MᵀM = I)
+    for (int i = 0; i < 2; i++) {
+      for (int j = i; j < 2; j++) {
+        float dotU = Ublock[0][i] * Ublock[0][j] + Ublock[1][i] * Ublock[1][j];
+        float dotV = Vblock[0][i] * Vblock[0][j] + Vblock[1][i] * Vblock[1][j];
+        float expected = (i == j) ? 1.0f : 0.0f;
+        REQUIRE_THAT(dotU, Catch::Matchers::WithinAbs(expected, 1e-3f));
+        REQUIRE_THAT(dotV, Catch::Matchers::WithinAbs(expected, 1e-3f));
+      }
+    }
+
+    // 3. Ublock · diag(sigma) · Vblockᵀ reproduces B = [[a,b],[0,d]]
+    //    (C[i][j] = sum_k U[i][k] * sigma[k] * V[j][k])
+    float C[2][2] = {{0}, {0}};
+    for (int i = 0; i < 2; i++)
+      for (int j = 0; j < 2; j++)
+        for (int k = 0; k < 2; k++)
+          C[i][j] += Ublock[i][k] * sigma[k] * Vblock[j][k];
+    REQUIRE_THAT(C[0][0], Catch::Matchers::WithinAbs(tc.a, 1e-2f));
+    REQUIRE_THAT(C[0][1], Catch::Matchers::WithinAbs(tc.b, 1e-2f));
+    REQUIRE_THAT(C[1][0], Catch::Matchers::WithinAbs(0.0f, 1e-2f));
+    REQUIRE_THAT(C[1][1], Catch::Matchers::WithinAbs(tc.d, 1e-2f));
+  }
+}
+
+// ============================================================================
+// TEST: JacobiEigenSymmetric — cyclic Jacobi eigenvalue decomposition
+// ============================================================================
+// Reference eigenvalues generated with scipy.linalg.eigvalsh (desc).
+TEST_CASE("SVD Building Block: JacobiEigenSymmetric", "[Matrix][SVD]") {
+  struct CaseJac {
+    float S[5][5];
+    uint8_t n;
+    float refEig[5];
+  };
+
+  // (i) T = BᵀB from a real bidiagonalization (3×3)
+  float T3[5][5] = {
+      {65.999993f, -124.470864f, 0.0f, 0.0f, 0.0f},
+      {-124.470864f, 237.877008f, -0.499065f, 0.0f, 0.0f},
+      {0.0f, -0.499065f, 0.122959f, 0.0f, 0.0f},
+      {0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
+      {0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
+  };
+  // (ii) random-looking 3×3 symmetric (seed 42)
+  float S3[5][5] = {
+      {0.304717f, -0.04971f, 0.439146f, 0.0f, 0.0f},
+      {-0.04971f, -1.951035f, -0.809211f, 0.0f, 0.0f},
+      {0.439146f, -0.809211f, -0.016801f, 0.0f, 0.0f},
+      {0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
+      {0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
+  };
+  // (iii) random-looking 4×4 symmetric (seed 42)
+  float S4[5][5] = {
+      {-0.853044f, 1.00332f, -0.090545f, -0.307449f, 0.0f},
+      {1.00332f, 0.467509f, 0.009579f, 0.795646f, 0.0f},
+      {-0.090545f, 0.009579f, -0.049926f, -0.169696f, 0.0f},
+      {-0.307449f, 0.795646f, -0.169696f, -0.428328f, 0.0f},
+      {0.0f, 0.0f, 0.0f, 0.0f, 0.0f},
+  };
+
+  float refs[3][5] = {
+      {303.195295f, 0.765908223f, 0.0387564408f, 0, 0},
+      {0.7227162f, -0.13661881f, -2.24921639f, 0, 0},
+      {1.22127596f, -0.01555681f, -0.31307273f, -1.75643542f, 0},
+  };
+  uint8_t ns[3] = {3, 3, 4};
+  float (*mats[3])[5] = {T3, S3, S4};
+  float maxAbs[3] = {237.877008f, 1.951035f, 1.00332f};
+
+  for (int c = 0; c < 3; c++) {
+    float T[5][5];
+    for (int i = 0; i < 5; i++)
+      for (int j = 0; j < 5; j++)
+        T[i][j] = mats[c][i][j];
+    float S_orig[5][5];
+    for (int i = 0; i < 5; i++)
+      for (int j = 0; j < 5; j++)
+        S_orig[i][j] = mats[c][i][j];
+
+    float evals[5] = {0};
+    // JacobiEigenSymmetric operates on Matrix<N,N> — copy the raw test
+    // data in, run the solver, copy the eigenvector matrix back out.
+    Matrix<5, 5> Tm{0};
+    for (int i = 0; i < 5; i++)
+      for (int j = 0; j < 5; j++)
+        Tm[i][j] = T[i][j];
+    Matrix<5, 5> Vm{0};
+    SVD::JacobiEigenSymmetric(Tm, ns[c], evals, Vm);
+    float V[5][5] = {{0}};
+    for (int i = 0; i < 5; i++)
+      for (int j = 0; j < 5; j++)
+        V[i][j] = Vm[i][j];
+
+    // 1. Sorted eigenvalues match scipy
+    float sorted[5] = {0};
+    for (int i = 0; i < ns[c]; i++) sorted[i] = evals[i];
+    // Sort descending to match the scipy reference order
+    for (int i = 0; i < ns[c] - 1; i++) {
+      int maxIdx = i;
+      for (int j = i + 1; j < ns[c]; j++)
+        if (sorted[j] > sorted[maxIdx])
+          maxIdx = j;
+      if (maxIdx != i) {
+        float t = sorted[i];
+        sorted[i] = sorted[maxIdx];
+        sorted[maxIdx] = t;
+      }
+    }
+    for (int i = 0; i < ns[c]; i++) {
+      if (fabsf(refs[c][i]) > 0.01f) {
+        REQUIRE_THAT(sorted[i],
+                     Catch::Matchers::WithinRel(refs[c][i], 1e-3f));
+      } else {
+        REQUIRE_THAT(sorted[i], Catch::Matchers::WithinAbs(refs[c][i], 1e-3f));
+      }
+    }
+
+    // 2. V is orthogonal (VᵀV = I on the n×n part)
+    for (int i = 0; i < ns[c]; i++) {
+      for (int j = i; j < ns[c]; j++) {
+        float dot = 0.0f;
+        for (int k = 0; k < ns[c]; k++) dot += V[k][i] * V[k][j];
+        float expected = (i == j) ? 1.0f : 0.0f;
+        REQUIRE_THAT(dot, Catch::Matchers::WithinAbs(expected, 1e-3f));
+      }
+    }
+
+    // 3. Residual ‖S_orig·V − V·diag(evals)‖ small
+    //    (col i of S_orig·V must equal evals_i · col i of V)
+    float residual = 0.0f;
+    for (int i = 0; i < ns[c]; i++) {
+      for (int r = 0; r < ns[c]; r++) {
+        float Sv = 0.0f;
+        for (int k = 0; k < ns[c]; k++) Sv += S_orig[r][k] * V[k][i];
+        float diff = Sv - evals[i] * V[r][i];
+        residual += diff * diff;
+      }
+    }
+    residual = sqrtf(residual);
+    REQUIRE_THAT(residual,
+                 Catch::Matchers::WithinAbs(0.0f,
+                                            1e-2f * maxAbs[c]));
+  }
+}
+
+// ============================================================================
+// TEST: DeflateBidiagonal / BidiagonalIsDiagonal
+// ============================================================================
+TEST_CASE("SVD Building Block: DeflateBidiagonal and BidiagonalIsDiagonal",
+          "[Matrix][SVD]") {
+  float tol = 1e-8f;
+
+  // IsDiagonal: true on a diagonal matrix
+  {
+    Matrix<5, 5> W{10.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 5.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 2.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+    REQUIRE(SVD::BidiagonalIsDiagonal(W, 5, tol));
+  }
+
+  // IsDiagonal: false when a superdiagonal is significant
+  {
+    Matrix<5, 5> W{10.0f, 1e-3f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 5.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 2.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
+    REQUIRE_FALSE(SVD::BidiagonalIsDiagonal(W, 5, tol));
+  }
+
+  // Deflate: small superdiagonals zeroed, significant ones kept
+  {
+    Matrix<5, 5> W{1.0f, 0.5f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 2.0f, 1e-9f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 3.0f, 0.3f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 4.0f, 1e-12f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 5.0f};
+    SVD::DeflateBidiagonal(W, 5, tol);
+    REQUIRE_THAT(W.Get(0, 1), Catch::Matchers::WithinAbs(0.5f, 1e-6f));
+    REQUIRE(W.Get(1, 2) == 0.0f);
+    REQUIRE_THAT(W.Get(2, 3), Catch::Matchers::WithinAbs(0.3f, 1e-6f));
+    REQUIRE(W.Get(3, 4) == 0.0f);
+    // Diagonal untouched
+    REQUIRE_THAT(W.Get(0, 0), Catch::Matchers::WithinAbs(1.0f, 1e-6f));
+    REQUIRE_THAT(W.Get(4, 4), Catch::Matchers::WithinAbs(5.0f, 1e-6f));
+    // NOT fully diagonal: significant superdiagonals (0.5, 0.3) remain
+    REQUIRE_FALSE(SVD::BidiagonalIsDiagonal(W, 5, tol));
+  }
+
+  // Deflate on an already-diagonal-ish matrix makes IsDiagonal true
+  {
+    Matrix<5, 5> W{1.0f, 1e-9f, 0.0f, 0.0f, 0.0f,
+                   0.0f, 2.0f, 1e-11f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 3.0f, 0.0f, 0.0f,
+                   0.0f, 0.0f, 0.0f, 4.0f, 1e-10f,
+                   0.0f, 0.0f, 0.0f, 0.0f, 5.0f};
+    SVD::DeflateBidiagonal(W, 5, tol);
+    REQUIRE(SVD::BidiagonalIsDiagonal(W, 5, tol));
+  }
+}

unit-tests/svd-integration-test.cpp

@@ -62,7 +62,7 @@ TEST_CASE("SVD Integration: 2x2 [[1,2],[3,4]]", "[Matrix][SVD][Integration]") {
       err += diff * diff;
     }
   err = sqrtf(err);
-  REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 
   std::cout << "SVD 2x2 [[1,2],[3,4]]:\n";
   std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0)
@@ -86,8 +86,8 @@ TEST_CASE("SVD Integration: 3x3 diagonal [10,5,2]",
   // U and Vt should be identity (or close) for diagonal matrix
   float uErr = frobeniusNorm(U - Matrix<3, 3>{1, 0, 0, 0, 1, 0, 0, 0, 1});
   float vtErr = frobeniusNorm(Vt - Matrix<3, 3>{1, 0, 0, 0, 1, 0, 0, 0, 1});
-  REQUIRE_THAT(uErr, Catch::Matchers::WithinRel(0.0f, 1e-2f));
+  REQUIRE_THAT(uErr, Catch::Matchers::WithinAbs(0.0f, 1e-2f));
-  REQUIRE_THAT(vtErr, Catch::Matchers::WithinRel(0.0f, 1e-2f));
+  REQUIRE_THAT(vtErr, Catch::Matchers::WithinAbs(0.0f, 1e-2f));
 }
 
 TEST_CASE("SVD Integration: 3x3 rank-deficient [[1,2,3],[4,5,6],[7,8,9]]",
@@ -120,7 +120,7 @@ TEST_CASE("SVD Integration: 3x3 rank-deficient [[1,2,3],[4,5,6],[7,8,9]]",
       err += diff * diff;
     }
   err = sqrtf(err);
-  REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-2f));
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-2f));
 
   std::cout << "SVD 3x3 rank-deficient:\n";
   std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) << ", "
@@ -155,7 +155,7 @@ TEST_CASE("SVD Integration: tall 4x3 matrix", "[Matrix][SVD][Integration]") {
       err += diff * diff;
     }
   err = sqrtf(err);
-  REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-2f));
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-2f));
 
   std::cout << "SVD tall 4x3:\n";
   std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) << ", "
@@ -175,42 +175,25 @@ TEST_CASE("SVD Integration: wide 3x5 matrix", "[Matrix][SVD][Integration]") {
   REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.46540f, 1e-2f));
   REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-2f));
 
-  // Check reconstruction: U (3x5) * diag(sigma) (5x3) = 3x3, then * Vt (3x5) =
+  // Check reconstruction: A (3x5) = U * Sigma * Vt, where U (3x5) has
-  // 3x5
+  // its meaningful part in the first 3 columns, sigma (5x1) in the
-  Matrix<3, 5> recon{0};
+  // first 3 entries, and Vt (5x5) in its first 3 rows (right
-  Matrix<3, 5> Usig{0};
+  // singular vectors as rows). So:
-  for (int i = 0; i < 3; i++)
+  //   A[i][j] = sum_k U[i][k] * sigma[k] * Vt[k][j]
-    for (int j = 0; j < 5; j++)
-      Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0);
-
-  // For wide matrix: A = U * Sigma * Vt where U is 3x5, Sigma is 5x5
-  // (diagonal), Vt is 5x5 But our implementation returns sigma as 3x1 and Vt as
-  // 3x5 So we need: recon = Usig (3x5) * Vt (3x5)^T ... no that doesn't work
-  // either The SVD for wide matrices is: A = U * Sigma * Vt where:
-  //   U is m×m (3×3), Sigma is m×n (3×5), Vt is n×n (5×5)
-  // But our API returns U as m×n (3×5), sigma as n×1 (3×1), Vt as n×n (3×5)
-  // So: recon = U (3x5) * diag(sigma) (5x5) * Vt (5x5)^T ...
-  // Actually, looking at the implementation, for wide matrices we swap roles.
-  // Let me just check reconstruction using the actual dimensions returned.
-
-  // For wide matrix: A (3x5) = U (3x5) * diag(sigma) (5x5 padded) * Vt (5x5)
-  // But our API returns Vt as 3x5, not 5x5
-  // The implementation stores: U = QR[:,0:p]^T (3x5), sigma (3x1), Vt =
-  // QL[:,0:p]^T (3x5) Reconstruction: A[i][j] = sum_k U[i][k]*sigma[k]*Vt[j][k]
 
   float err2 = 0.0f;
   for (int i = 0; i < 3; i++) {
     for (int j = 0; j < 5; j++) {
       float recon_val = 0.0f;
       for (int k = 0; k < 3; k++) {
-        recon_val += U.Get(i, k) * sigma.Get(k, 0) * Vt.Get(j, k);
+        recon_val += U.Get(i, k) * sigma.Get(k, 0) * Vt.Get(k, j);
       }
       float diff = recon_val - A.Get(i, j);
       err2 += diff * diff;
     }
   }
   err2 = sqrtf(err2);
-  REQUIRE_THAT(err2, Catch::Matchers::WithinRel(0.0f, 1e-2f));
+  REQUIRE_THAT(err2, Catch::Matchers::WithinAbs(0.0f, 1e-2f));
 
   std::cout << "SVD wide 3x5:\n";
   std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0) << ", "
@@ -230,7 +213,7 @@ TEST_CASE("SVD Integration: identity 3x3", "[Matrix][SVD][Integration]") {
   REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(1.0f, 1e-3f));
 
   float err = frobeniusNorm(U - Matrix<3, 3>{1, 0, 0, 0, 1, 0, 0, 0, 1});
-  REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-2f));
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-2f));
 }
 
 TEST_CASE("SVD Integration: symmetric positive definite 2x2 [[5,3],[3,5]]",
@@ -261,9 +244,120 @@ TEST_CASE("SVD Integration: symmetric positive definite 2x2 [[5,3],[3,5]]",
       err += diff * diff;
     }
   err = sqrtf(err);
-  REQUIRE_THAT(err, Catch::Matchers::WithinRel(0.0f, 1e-3f));
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
 
   std::cout << "SVD SPD 2x2 [[5,3],[3,5]]:\n";
   std::cout << "Sigma: [" << sigma.Get(0, 0) << ", " << sigma.Get(1, 0)
             << "]\n";
 }
+
+// ----------------------------------------------------------------------------
+// Matrix::SVD member wrapper (delegates to SVD::SVD)
+// ----------------------------------------------------------------------------
+
+/**
+ * Reconstruction error ‖U·diag(sigma)·Vᵀ − A‖_F. Zero-padded entries of
+ * U/sigma/Vt (wide/tall cases) are zero by the output conventions, so the
+ * full product equals U[:, :k]·diag(sigma[:k])·Vt[:k, :].
+ */
+template <uint8_t rows, uint8_t columns>
+static float svdReconstructionError(const Matrix<rows, columns> &A,
+                                    const Matrix<rows, columns> &U,
+                                    const Matrix<columns, 1> &sigma,
+                                    const Matrix<columns, columns> &Vt) {
+  Matrix<rows, columns> recon{0};
+  Matrix<rows, columns> Usig{0};
+  for (int i = 0; i < rows; i++)
+    for (int j = 0; j < columns; j++)
+      Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0);
+  Usig.Mult(Vt, recon);
+  float err = 0.0f;
+  for (int i = 0; i < rows; i++)
+    for (int j = 0; j < columns; j++) {
+      float diff = recon.Get(i, j) - A.Get(i, j);
+      err += diff * diff;
+    }
+  return sqrtf(err);
+}
+
+/**
+ * Orthonormality of the first k columns of M: the k×k leading block of
+ * MᵀM must equal I_k. (For a tall SVD, U has k = min(rows, cols)
+ * meaningful columns and this is the full UᵀU.)
+ */
+template <uint8_t r, uint8_t c>
+static bool leadingColumnsOrthonormal(const Matrix<r, c> &M, uint8_t k,
+                                      float tol = 1e-4f) {
+  Matrix<c, r> Mt = M.Transpose();
+  Matrix<c, c> MtM{0};
+  Mt.Mult(M, MtM);
+  for (int i = 0; i < k; i++)
+    for (int j = 0; j < k; j++) {
+      float expected = (i == j) ? 1.0f : 0.0f;
+      if (fabsf(MtM.Get(i, j) - expected) > tol)
+        return false;
+    }
+  return true;
+}
+
+/**
+ * Orthonormality of the first k rows of M: the k×k leading block of
+ * M·Mᵀ must equal I_k. (Vᵀ may have zero-padded trailing rows in the
+ * wide case, so check only the meaningful leading block.)
+ */
+template <uint8_t r, uint8_t c>
+static bool leadingRowsOrthonormal(const Matrix<r, c> &M, uint8_t k,
+                                   float tol = 1e-4f) {
+  Matrix<c, r> Mt = M.Transpose();
+  Matrix<r, r> MMt{0};
+  M.Mult(Mt, MMt);
+  for (int i = 0; i < k; i++)
+    for (int j = 0; j < k; j++) {
+      float expected = (i == j) ? 1.0f : 0.0f;
+      if (fabsf(MMt.Get(i, j) - expected) > tol)
+        return false;
+    }
+  return true;
+}
+
+TEST_CASE("Matrix::SVD wrapper: 3x2 tall [[1,2],[3,4],[5,6]]",
+          "[Matrix][SVD][Wrapper]") {
+  Matrix<3, 2> A{1, 2, 3, 4, 5, 6};
+  Matrix<3, 2> U{0};
+  Matrix<2, 1> sigma{0};
+  Matrix<2, 2> Vt{0};
+
+  A.SVD(U, sigma, Vt);
+
+  // Reference singular values from numpy: [9.52552, 0.514301]
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(9.52552f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(0.514301f, 1e-3f));
+
+  REQUIRE(leadingColumnsOrthonormal(U, 2));
+  REQUIRE(leadingRowsOrthonormal(Vt, 2));
+
+  float err = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("Matrix::SVD wrapper: 2x3 wide [[1,2,3],[4,5,6]]",
+          "[Matrix][SVD][Wrapper]") {
+  Matrix<2, 3> A{1, 2, 3, 4, 5, 6};
+  Matrix<2, 3> U{0};
+  Matrix<3, 1> sigma{0};
+  Matrix<3, 3> Vt{0};
+
+  A.SVD(U, sigma, Vt);
+
+  // Reference singular values from numpy: [9.50803, 0.77287]; the third
+  // entry (wide-matrix padding) must be zero.
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(9.50803f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(0.77287f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+
+  REQUIRE(leadingColumnsOrthonormal(U, 2));
+  REQUIRE(leadingRowsOrthonormal(Vt, 2));
+
+  float err = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-3f));
+}

unit-tests/svd-reference-values.py

@@ -401,6 +401,37 @@ def main():
         ("Zero 3x3", np.zeros((3,3))),
         ("Col vector 2x1", np.array([[3],[4]], dtype=np.float64)),
         ("Row vector 1x2", np.array([[3,4]], dtype=np.float64)),
+        # Large-size instantiation cases (N > 5). Literals MUST match the
+        # C++ test matrices in unit-tests/matrix-tests.cpp exactly, and the
+        # C++ references use float32 inputs: cast to float32 before svd().
+        ("Tall 7x5", np.array([
+            [-0.7528, 2.7043, 1.392, 0.592, -2.0639],
+            [-2.064, -2.6515, 2.1971, 0.6067, 1.2484],
+            [-2.8765, 2.8195, 1.9947, -1.726, -1.9091],
+            [-1.8996, -1.1745, 0.1485, -0.4083, -1.2526],
+            [0.6711, -2.163, -1.2471, -0.8018, -0.2636],
+            [1.7111, -1.802, 0.0854, 0.5545, -2.7213],
+            [0.6453, -1.9769, -2.6097, 2.6933, 2.7938]], dtype=np.float32)),
+        ("Square 6x6", np.array([
+            [1.2336, -0.7815, -1.6093, 0.7369, -0.2394, -1.5118],
+            [-0.0193, -1.8624, 1.6373, -0.9649, 0.6501, -0.7532],
+            [0.0803, 0.1868, -1.2606, 1.8783, 1.1005, 1.758],
+            [1.5793, 0.3916, 1.6875, -1.646, -1.2161, -1.8191],
+            [-0.6987, -0.4453, -0.9146, 1.315, -0.573, -0.8763],
+            [0.1708, -1.4363, 1.2088, -1.7018, 1.089, 1.9475]], dtype=np.float32)),
+        ("Wide 5x8", np.array([
+            [-1.5064, -2.4724, 1.5773, 1.0343, 1.145, 1.3564, -2.1298, -0.7077],
+            [-1.9207, 1.8155, 0.6165, -0.8455, -2.1822, -0.9451, -0.8741, 1.148],
+            [0.6878, 1.9361, -0.1389, -1.902, 1.0662, 1.3039, 0.3064, 1.3548],
+            [-0.031, 0.1137, -0.3623, -2.3729, -1.9605, -2.3429, 0.6821, -0.9282],
+            [0.0429, 2.0378, -1.2535, -0.4481, 1.2778, -1.356, -2.1151, -1.0512]], dtype=np.float32)),
+        ("Tall 6x4 rank-def", np.array([
+            [-0.086904, 1.410225, 1.308323, 2.234762],
+            [0.022123, 0.896751, 0.324176, 0.773607],
+            [-0.473015, 1.555111, 0.290059, 1.157726],
+            [-0.78371, 1.398884, -1.930606, -1.548717],
+            [0.201518, -0.626835, 0.976596, 0.875294],
+            [-1.24206, 1.60595, -3.078089, -2.73695]], dtype=np.float32)),
     ]
     
     for name, A in test_matrices: