Working on breaking up the steps into manageable chunks

Reviewed-on: #7

.gitea/workflows/Merge-Checker.yaml

@@ -0,0 +1,102 @@
+name: Merge-Checker
+
+on:
+  pull_request:
+    branches: ["**"]
+
+jobs:
+  build_and_test:
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Checkout source code
+        uses: actions/checkout@v3
+        with:
+          persist-credentials: true
+          fetch-depth: 0
+
+      - name: Install dependencies (CMake + Ninja + build tools)
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y cmake ninja-build build-essential time git
+
+      - name: Configure project with CMake
+        run: cmake -G Ninja -S . -B build/
+
+      - name: Build with Ninja
+        run: ninja -C build/
+
+      - name: Run all unit tests except matrix-timing-tests
+        run: |
+          for test_exec in build/unit-tests/matrix-tests build/unit-tests/quaternion-tests build/unit-tests/vector-3d-tests; do
+            if [ -x "$test_exec" ]; then
+              echo "Running $test_exec"
+              "$test_exec"
+            else
+              echo "Warning: $test_exec not found or not executable"
+            fi
+          done
+      - name: Run matrix-timing-tests
+        run: |
+          mkdir -p unit-tests/timing-results
+          if [ -x build/unit-tests/matrix-timing-tests ]; then
+            echo "Running matrix-timing-tests with timing"
+            /usr/bin/time -v build/unit-tests/matrix-timing-tests -d yes &> unit-tests/timing-results/matrix-timing-tests.txt
+            cat unit-tests/timing-results/matrix-timing-tests.txt
+          else
+            echo "matrix-timing-tests executable not found or not executable"
+            exit 1
+          fi
+
+      - name: Compare timing results
+        id: check_diff
+        run: |
+          git show origin/${{ github.event.pull_request.head.ref }}:unit-tests/timing-results/matrix-timing-tests.txt > old.txt || echo "" > old.txt
+          cp unit-tests/timing-results/matrix-timing-tests.txt new.txt
+
+          echo "Comparing timing results for changes ≥ 0.1s (ignoring 'Timing Tests' lines)..."
+
+          changed=0
+
+          awk -v changed_ref=/tmp/timings_changed.flag '
+            BEGIN {
+              change_threshold = 0.1
+            }
+            FILENAME == "old.txt" && /^[0-9]+\.[0-9]+ s: / {
+              label = substr($0, index($0, ":") + 2)
+              if (label != "Timing Tests") {
+                label_times[label] = $1
+              }
+            }
+            FILENAME == "new.txt" && /^[0-9]+\.[0-9]+ s: / {
+              new_time = $1
+              label = substr($0, index($0, ":") + 2)
+              if (label == "Timing Tests") next
+
+              old_time = label_times[label]
+              delta = new_time - old_time
+              if (delta < 0) delta = -delta
+
+              if (old_time != "" && delta >= change_threshold) {
+                printf "⚠️  %.3f s → %.3f s: %s (Δ=%.3f s)\n", old_time, new_time, label, delta
+                system("touch " changed_ref)
+              } else if (old_time == "") {
+                printf "🆕 New timing entry: %.3f s: %s\n", new_time, label
+                system("touch " changed_ref)
+              }
+            }
+            END {
+              if (!system("test -f " changed_ref)) {
+                exit 0
+              } else {
+                print "✅ Timings haven’t changed significantly (Δ < 0.1s)."
+                exit 0
+              }
+            }
+          ' old.txt new.txt
+
+          if [ -f /tmp/timings_changed.flag ]; then
+            echo "timings_changed=true" >> $GITHUB_OUTPUT
+          else
+            echo "timings_changed=false" >> $GITHUB_OUTPUT
+          fi

.gitignore

@@ -1,2 +1,2 @@
 build/
-venv/
+.cache/

.vscode/launch.json

@@ -27,16 +27,25 @@
             "internalConsoleOptions": "openOnSessionStart"
         },
         {
-            "name": "Run Matrix Unit Tests",
+            "name": "Debug Quaternion Unit Tests",
-            "type": "cpp",
+            "type": "cppdbg",
             "request": "launch",
-            "program": "${workspaceFolder}/build/unit-tests/matrix-tests",
+            "program": "${workspaceFolder}/build/unit-tests/quaternion-tests",
             "args": [],
             "stopAtEntry": false,
             "cwd": "${workspaceFolder}",
             "environment": [],
             "externalConsole": false,
-            "preLaunchTask": "build_tests",  // Compile unit tests before running
+            "MIMode": "gdb",
+            "miDebuggerPath": "/usr/bin/gdb", // Adjust to your debugger path
+            "setupCommands": [
+                {
+                    "description": "Enable pretty-printing for gdb",
+                    "text": "-enable-pretty-printing",
+                    "ignoreFailures": true
+                }
+            ],
+            "preLaunchTask": "build_tests", // Task to compile unit tests
             "internalConsoleOptions": "openOnSessionStart"
         }
     ]

.vscode/settings.json

@@ -1,8 +1,5 @@
 {
     "C_Cpp.intelliSenseEngine": "default",
-    "clangd.arguments": [
-        "--include-directory=build/unit-tests"
-    ],
     "C_Cpp.default.intelliSenseMode": "linux-gcc-x64",
     "files.associations": {
         "*.h": "cpp",
@@ -71,8 +68,16 @@
         "typeinfo": "cpp",
         "variant": "cpp",
         "shared_mutex": "cpp",
-        "complex": "cpp"
+        "charconv": "cpp",
+        "format": "cpp",
+        "csignal": "cpp",
+        "span": "cpp"
     },
-    "clangd.enable": false,
+    "clangd.enable": true,
-    "C_Cpp.dimInactiveRegions": false
+    "C_Cpp.dimInactiveRegions": false,
+    "editor.defaultFormatter": "xaver.clang-format",
+    "clangd.inactiveRegions.useBackgroundHighlight": false,
+    "clangd.arguments": [
+        "--compile-commands-dir=${workspaceFolder}/build"
+    ],
 }

.vscode/tasks.json

@@ -4,12 +4,14 @@
         {
             "label": "build_tests",
             "type": "shell",
-            "command": "cd build && ninja matrix-tests",
+            "command": "cd build && ninja",
             "group": {
                 "kind": "build",
                 "isDefault": true
             },
-            "problemMatcher": ["$gcc"],
+            "problemMatcher": [
+                "$gcc"
+            ],
             "detail": "Generated task to build unit test executable"
         }
     ]

CMakeLists.txt

@@ -1,40 +1,21 @@
-cmake_minimum_required(VERSION 3.6)
+cmake_minimum_required (VERSION 3.11)
 project(Vector3D)
 
+add_subdirectory(src)
 add_subdirectory(unit-tests)
 
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 17)
 
+add_compile_options(-Wall -Wextra -Wpedantic)
 add_compile_options (-fdiagnostics-color=always)
+set(CMAKE_COLOR_DIAGNOSTICS ON)
 
-# Vector3d
+include(FetchContent)
-add_library(Vector3D 
+ 
-    STATIC
+FetchContent_Declare(
-    Vector3D.hpp
+    Catch2
+    GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+    GIT_TAG v3.8.0 # or a later release
 )
  
-set_target_properties(Vector3D
+FetchContent_MakeAvailable(Catch2)
-    PROPERTIES
-    LINKER_LANGUAGE CXX
-)
-
-target_include_directories(Vector3D PUBLIC
-    include
-)
-
-# Matrix
-add_library(Matrix 
-    STATIC
-    Matrix.hpp
-    Matrix.cpp
-)
-
-set_target_properties(Matrix
-    PROPERTIES
-    LINKER_LANGUAGE CXX
-)
-
-target_include_directories(Matrix 
-    PUBLIC
-    .
-)

README.md

@@ -1 +1,12 @@
-A Simple matrix math library focused on embedded development which avoids and heap memory allocation unless you explicitly ask for it.
+# Introduction
+This matrix math library is focused on embedded development and avoids any heap memory allocation unless you explicitly ask for it.
+It uses templates to pre-allocate matrices on the stack.
+
+# Building
+1. Initialize the repositiory with the command:
+```bash
+cmake -S . -B build -G Ninja
+```
+
+2. Go into the build folder and run `ninja`
+3. That's it. You can test out the build by running `./unit-tests/matrix-tests`

Vector3D.hpp

@@ -1,81 +0,0 @@
-#pragma once
-
-#include <cstdint>
-#include <cmath>
-#include <type_traits>
-
-template <typename Type>
-class V3D{
-    public:
-    constexpr V3D(const V3D& other):
-    x(other.x),
-    y(other.y),
-    z(other.z){
-        static_assert(std::is_arithmetic<Type>::value, "Type must be a number");
-    }
-    
-    constexpr V3D(Type x=0, Type y=0, Type z=0): 
-    x(x), 
-    y(y), 
-    z(z){
-        static_assert(std::is_arithmetic<Type>::value, "Type must be a number");
-    }
-
-    template <typename OtherType>
-    constexpr V3D(const V3D<OtherType> other):
-    x(static_cast<Type>(other.x)),
-    y(static_cast<Type>(other.y)),
-    z(static_cast<Type>(other.z)){
-        static_assert(std::is_arithmetic<Type>::value, "Type must be a number");
-        static_assert(std::is_arithmetic<OtherType>::value, "OtherType must be a number");
-    }
-
-    V3D& operator=(const V3D &other){
-        this->x = other.x;
-        this->y = other.y;
-        this->z = other.z;
-        return *this;
-    }
-
-    V3D& operator+=(const V3D &other){
-        this->x += other.x;
-        this->y += other.y;
-        this->z += other.z;
-        return *this;
-    }
-
-    V3D& operator-=(const V3D &other){
-        this->x -= other.x;
-        this->y -= other.y;
-        this->z -= other.z;
-        return *this;
-    }
-
-    V3D& operator/=(const Type scalar){
-        if(scalar == 0){
-            return *this;
-        }
-        this->x /= scalar;
-        this->y /= scalar;
-        this->z /= scalar;
-        return *this;
-    }
-
-    V3D& operator*=(const Type scalar){
-        this->x *= scalar;
-        this->y *= scalar;
-        this->z *= scalar;
-        return *this;
-    }
-
-    bool operator==(const V3D &other){
-        return this->x == other.x && this->y == other.y && this->z == other.z;
-    }
-
-    float magnitude(){
-        return std::sqrt(static_cast<float>(this->x * this->x + this->y * this->y + this->z * this->z));
-    }
-    Type x;
-    Type y;
-    Type z;
-};

library.json

@@ -0,0 +1,20 @@
+{
+    "name": "Vector3D",
+    "version": "1.0.0",
+    "description": "Contains a V3D object for easy 3d vector math and a Matrix object for more complicated linear algebra operations.",
+    "keywords": "linear algebra, vector, matrix, 3D",
+    "repository": {
+        "type": "git",
+        "url": "https://github.com/Cynopolis/Vector3D.git"
+    },
+    "authors": [
+        {
+            "name": "Cynopolis",
+            "email": "megaveganzombie@gmail.com",
+            "url": "https://github.com/Cynopolis"
+        }
+    ],
+    "license": "None Yet",
+    "frameworks": "*",
+    "platforms": "*"
+}

qr-decom.py

@@ -1,159 +0,0 @@
-import numpy as np
-
-# QR decomposition using the householder reflection method
-def householder_reflection(A):
-    """
-    Perform QR decomposition using Householder reflection.
-    
-    Arguments:
-    A -- A matrix to be decomposed (m x n).
-    
-    Returns:
-    Q -- Orthogonal matrix (m x m).
-    R -- Upper triangular matrix (m x n).
-    """
-    A = A.astype(float)  # Ensure the matrix is of type float
-    m, n = A.shape
-    Q = np.eye(m)  # Initialize Q as an identity matrix
-    R = A.copy()  # R starts as a copy of A
-
-    # Apply Householder reflections for each column
-    for k in range(n):
-        # Step 1: Compute the Householder vector
-        x = R[k:m, k]
-        e1 = np.zeros_like(x)
-        e1[0] = np.linalg.norm(x) if x[0] >= 0 else -np.linalg.norm(x)
-        v = x + e1
-        v = v / np.linalg.norm(v)  # Normalize v
-        
-        # Step 2: Apply the reflection to the matrix
-        R[k:m, k:n] = R[k:m, k:n] - 2 * np.outer(v, v.T @ R[k:m, k:n])
-        
-        # Step 3: Apply the reflection to Q
-        Q[:, k:m] = Q[:, k:m] - 2 * np.outer(Q[:, k:m] @ v, v)
-
-    # The resulting Q and R are the QR decomposition
-    return Q, R
-
-# Example usage
-A = np.array([[12, -51, 4], 
-              [6, 167, -68], 
-              [-4, 24, -41]])
-
-Q, R = householder_reflection(A)
-print("Q matrix:")
-print(Q)
-print("\nR matrix:")
-print(R)
-print("Multiplied Together:")
-print(Q@R)
-
-def svd_decomposition(A):
-    """
-    Perform Singular Value Decomposition (SVD) from scratch.
-    
-    Arguments:
-    A -- The matrix to be decomposed (m x n).
-    
-    Returns:
-    U -- Orthogonal matrix of left singular vectors (m x m).
-    Sigma -- Diagonal matrix of singular values (m x n).
-    Vt -- Orthogonal matrix of right singular vectors (n x n).
-    """
-    # Step 1: Compute A^T A
-    AtA = np.dot(A.T, A)  # A transpose multiplied by A
-    
-    # Step 2: Compute the eigenvalues and eigenvectors of A^T A
-    eigenvalues, V = np.linalg.eig(AtA)
-    
-    # Step 3: Sort eigenvalues in descending order and sort V accordingly
-    sorted_indices = np.argsort(eigenvalues)[::-1]  # Indices to sort eigenvalues in descending order
-    eigenvalues = eigenvalues[sorted_indices]
-    V = V[:, sorted_indices]
-    
-    # Step 4: Compute the singular values (sqrt of eigenvalues)
-    singular_values = np.sqrt(eigenvalues)
-    
-    # Step 5: Construct the Sigma matrix
-    m, n = A.shape
-    Sigma = np.zeros((m, n))  # Initialize Sigma as a zero matrix
-    for i in range(min(m, n)):
-        Sigma[i, i] = singular_values[i]  # Place the singular values on the diagonal
-    
-    # Step 6: Compute the U matrix using A * V = U * Sigma
-    U = np.dot(A, V)  # A * V gives us the unnormalized U
-    # Normalize the columns of U
-    for i in range(U.shape[1]):
-        U[:, i] = U[:, i] / singular_values[i]  # Normalize each column by the corresponding singular value
-    
-    # Step 7: Return U, Sigma, Vt
-    return U, Sigma, V.T  # V.T is the transpose of V
-
-# Example usage
-A = np.array([[12, -51, 4],
-              [6, 167, -68],
-              [-4, 24, -41]])
-
-U, Sigma, Vt = svd_decomposition(A)
-
-print("\nSVD DECOMPOSITION\nU matrix:")
-print(U)
-print("\nSigma matrix:")
-print(Sigma)
-print("\nVt matrix:")
-print(Vt)
-print("Multiplied together:")
-print(U@Sigma@Vt)
-
-def eigen_decomposition_qr(A, max_iter=1000, tol=1e-9):
-    """
-    Compute the eigenvalues and eigenvectors of a matrix A using the QR algorithm
-    with QR decomposition.
-
-    Arguments:
-    A -- A square matrix (n x n).
-    max_iter -- Maximum number of iterations for convergence (default 1000).
-    tol -- Tolerance for convergence (default 1e-9).
-    
-    Returns:
-    eigenvalues -- List of eigenvalues.
-    eigenvectors -- Matrix of eigenvectors.
-    """
-    # Make a copy of A to perform the iteration
-    A_copy = A.copy()
-    n = A_copy.shape[0]
-    
-    # Initialize the matrix for eigenvectors (this will accumulate the Q matrices)
-    eigenvectors = np.eye(n)
-    
-    # Perform QR iterations
-    for _ in range(max_iter):
-        # Perform QR decomposition on A_copy
-        Q, R = householder_reflection(A_copy)
-        
-        # Update A_copy to be R * Q (QR algorithm step)
-        A_copy = R @ Q
-        
-        # Accumulate the eigenvectors
-        eigenvectors = eigenvectors @ Q
-        
-        # Check for convergence: if the off-diagonal elements are small enough, we stop
-        off_diagonal_norm = np.linalg.norm(np.tril(A_copy, -1))  # Norm of the lower triangle (off-diagonal)
-        if off_diagonal_norm < tol:
-            break
-    
-    # The eigenvalues are the diagonal elements of the matrix A_copy
-    eigenvalues = np.diag(A_copy)
-    
-    return eigenvalues, eigenvectors
-
-# Example usage
-A = np.array([[12, -51, 4], 
-              [6, 167, -68], 
-              [-4, 24, -41]])
-
-eigenvalues, eigenvectors = eigen_decomposition_qr(A)
-
-
-print("\n\nEigenvalues:", eigenvalues)
-print("Eigenvectors:\n", eigenvectors)

src/CMakeLists.txt

@@ -0,0 +1,76 @@
+# Quaternion Interface
+add_library(vector-3d-intf
+    INTERFACE
+)
+
+target_include_directories(vector-3d-intf
+    INTERFACE
+    .
+)
+
+target_link_libraries(vector-3d-intf
+    INTERFACE
+)
+
+# Quaternion
+add_library(quaternion 
+    STATIC
+    Quaternion.cpp
+)
+
+target_link_libraries(quaternion
+    PUBLIC
+    vector-3d-intf
+    PRIVATE
+)
+
+set_target_properties(quaternion
+    PROPERTIES
+    LINKER_LANGUAGE CXX
+)
+
+# Vector3d
+add_library(vector-3d 
+    STATIC
+    Vector3D.cpp
+)
+
+target_link_libraries(vector-3d
+    PUBLIC
+    vector-3d-intf
+    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
+    SVD.cpp
+)
+
+target_link_libraries(svd
+    PUBLIC
+    vector-3d-intf
+    PRIVATE
+)
+
+set_target_properties(svd
+    PROPERTIES
+    LINKER_LANGUAGE CXX
+)

src/Matrix.cpp

@@ -1,3 +1,10 @@
+// This #ifndef section makes clangd happy so that it can properly do type hints
+// in this file
+#ifndef MATRIX_H_
+#define MATRIX_H_
+#include "Matrix.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"
@@ -5,18 +12,45 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdlib>
-#include <type_traits>
+#include <cstring>
-
-template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns>::Matrix(float value) {
-  this->Fill(value);
-}
 
 template <uint8_t rows, uint8_t columns>
 Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) {
   this->setMatrixToArray(array);
 }
 
+template <uint8_t rows, uint8_t columns>
+template <typename... Args,
+          std::enable_if_t<(std::is_arithmetic_v<Args> && ...), int>>
+Matrix<rows, columns>::Matrix(Args... args) {
+  constexpr uint16_t arraySize{static_cast<uint16_t>(rows) *
+                               static_cast<uint16_t>(columns)};
+
+  std::initializer_list<float> initList{static_cast<float>(args)...};
+  // if there is only one value, we actually want to do a fill
+  if (sizeof...(args) == 1) {
+    this->Fill(*initList.begin());
+  }
+  static_assert(sizeof...(args) == arraySize || sizeof...(args) == 1,
+                "You did not provide the right amount of initializers for this "
+                "matrix size");
+
+  // choose whichever buffer size is smaller for the copy length
+  uint32_t minSize =
+      std::min(arraySize, static_cast<uint16_t>(initList.size()));
+  memcpy(this->matrix.begin(), initList.begin(), minSize * sizeof(float));
+}
+
+template <uint8_t rows, uint8_t columns>
+Matrix<rows, columns> Matrix<rows, columns>::Identity() {
+  Matrix<rows, columns> identityMatrix{0};
+  uint32_t minDimension = std::min(rows, columns);
+  for (uint8_t idx{0}; idx < minDimension; idx++) {
+    identityMatrix[idx][idx] = 1;
+  }
+  return identityMatrix;
+}
+
 template <uint8_t rows, uint8_t columns>
 Matrix<rows, columns>::Matrix(const Matrix<rows, columns> &other) {
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
@@ -27,19 +61,6 @@ Matrix<rows, columns>::Matrix(const Matrix<rows, columns> &other) {
   }
 }
 
-template <uint8_t rows, uint8_t columns>
-template <typename... Args>
-Matrix<rows, columns>::Matrix(Args... args) {
-  constexpr uint16_t arraySize{static_cast<uint16_t>(rows) *
-                               static_cast<uint16_t>(columns)};
-
-  std::initializer_list<float> initList{static_cast<float>(args)...};
-  // choose whichever buffer size is smaller for the copy length
-  uint32_t minSize =
-      std::min(arraySize, static_cast<uint16_t>(initList.size()));
-  memcpy(this->matrix.begin(), initList.begin(), minSize * sizeof(float));
-}
-
 template <uint8_t rows, uint8_t columns>
 void Matrix<rows, columns>::setMatrixToArray(
     const std::array<float, rows * columns> &array) {
@@ -91,21 +112,18 @@ Matrix<rows, columns>::Mult(const Matrix<columns, other_columns> &other,
                             Matrix<rows, other_columns> &result) const {
   // allocate some buffers for all of our dot products
   Matrix<1, columns> this_row;
-  Matrix<rows, 1> other_column;
+  Matrix<columns, 1> other_column;
-  Matrix<1, rows> other_column_t;
 
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     // get our row
     this->GetRow(row_idx, this_row);
-    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
+    for (uint8_t column_idx{0}; column_idx < other_columns; column_idx++) {
       // get the other matrix'ss column
       other.GetColumn(column_idx, other_column);
-      // transpose the other matrix's column
-      other_column.Transpose(other_column_t);
 
       // the result's index is equal to the dot product of these two vectors
       result[row_idx][column_idx] =
-          Matrix<rows, columns>::dotProduct(this_row, other_column_t);
+          Matrix<rows, columns>::DotProduct(this_row, other_column.Transpose());
     }
   }
 
@@ -125,13 +143,13 @@ Matrix<rows, columns>::Mult(float scalar, Matrix<rows, columns> &result) const {
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns> &
+Matrix<rows, columns> Matrix<rows, columns>::Invert() const {
-Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
   // since all matrix sizes have to be statically specified at compile time we
   // can do this
   static_assert(rows == columns,
                 "Your matrix isn't square and can't be inverted");
 
+  Matrix<rows, columns> result{};
   // unfortunately we can't calculate this at compile time so we'll just reurn
   // zeros
   float determinant{this->Det()};
@@ -160,8 +178,8 @@ Matrix<rows, columns>::Invert(Matrix<rows, columns> &result) const {
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<columns, rows> &
+Matrix<columns, rows> Matrix<rows, columns>::Transpose() const {
-Matrix<rows, columns>::Transpose(Matrix<columns, rows> &result) const {
+  Matrix<columns, rows> result{};
   for (uint8_t column_idx{0}; column_idx < rows; column_idx++) {
     for (uint8_t row_idx{0}; row_idx < columns; row_idx++) {
       result[row_idx][column_idx] = this->Get(column_idx, row_idx);
@@ -173,9 +191,10 @@ Matrix<rows, columns>::Transpose(Matrix<columns, rows> &result) const {
 
 // explicitly define the determinant for a 2x2 matrix because it is definitely
 // the fastest way to calculate a 2x2 matrix determinant
-template <> float Matrix<0, 0>::Det() const { return 1e+6; }
+// template <>
-template <> float Matrix<1, 1>::Det() const { return this->matrix[0]; }
+// inline float Matrix<0, 0>::Det() const { return 1e+6; }
-template <> float Matrix<2, 2>::Det() const {
+template <> inline float Matrix<1, 1>::Det() const { return this->matrix[0]; }
+template <> inline float Matrix<2, 2>::Det() const {
   return this->matrix[0] * this->matrix[3] - this->matrix[1] * this->matrix[2];
 }
 
@@ -271,8 +290,13 @@ void Matrix<rows, columns>::ToString(std::string &stringBuffer) const {
 }
 
 template <uint8_t rows, uint8_t columns>
-std::array<float, columns> &Matrix<rows, columns>::
+const float *Matrix<rows, columns>::ToArray() const {
-operator[](uint8_t row_index) {
+  return this->matrix.data();
+}
+
+template <uint8_t rows, uint8_t columns>
+std::array<float, columns> &
+Matrix<rows, columns>::operator[](uint8_t row_index) {
   if (row_index > rows - 1) {
     // TODO: We should throw something here instead of failing quietly.
     row_index = 0;
@@ -284,38 +308,36 @@ operator[](uint8_t row_index) {
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns> &Matrix<rows, columns>::
+Matrix<rows, columns> &
-operator=(const Matrix<rows, columns> &other) {
+Matrix<rows, columns>::operator=(const Matrix<rows, columns> &other) {
-  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
+  memcpy(this->matrix.begin(), other.matrix.begin(),
-    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
+         rows * columns * sizeof(float));
-      this->matrix[row_idx * columns + column_idx] =
+
-          other.Get(row_idx, column_idx);
-    }
-  }
   // return a reference to ourselves so you can chain together these functions
   return *this;
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns> Matrix<rows, columns>::
+Matrix<rows, columns>
-operator+(const Matrix<rows, columns> &other) const {
+Matrix<rows, columns>::operator+(const Matrix<rows, columns> &other) const {
   Matrix<rows, columns> buffer{};
   this->Add(other, buffer);
   return buffer;
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns> Matrix<rows, columns>::
+Matrix<rows, columns>
-operator-(const Matrix<rows, columns> &other) const {
+Matrix<rows, columns>::operator-(const Matrix<rows, columns> &other) const {
   Matrix<rows, columns> buffer{};
   this->Sub(other, buffer);
   return buffer;
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns> Matrix<rows, columns>::
+template <uint8_t other_columns>
-operator*(const Matrix<rows, columns> &other) const {
+Matrix<rows, other_columns> Matrix<rows, columns>::operator*(
-  Matrix<rows, columns> buffer{};
+    const Matrix<columns, other_columns> &other) const {
+  Matrix<rows, other_columns> buffer{};
   this->Mult(other, buffer);
   return buffer;
 }
@@ -327,9 +349,25 @@ Matrix<rows, columns> Matrix<rows, columns>::operator*(float scalar) const {
   return buffer;
 }
 
+template <uint8_t rows, uint8_t columns>
+Matrix<rows, columns> Matrix<rows, columns>::operator/(float scalar) const {
+  Matrix<rows, columns> buffer = *this;
+  if (scalar == 0) {
+    buffer.Fill(1e+10);
+    return buffer;
+  }
+
+  for (uint8_t row = 0; row < rows; row++) {
+    for (uint8_t column = 0; column < columns; column++) {
+      buffer[row][column] /= scalar;
+    }
+  }
+  return buffer;
+}
+
 template <uint8_t rows, uint8_t columns>
 template <uint8_t vector_size>
-float Matrix<rows, columns>::dotProduct(const Matrix<1, vector_size> &vec1,
+float Matrix<rows, columns>::DotProduct(const Matrix<1, vector_size> &vec1,
                                         const Matrix<1, vector_size> &vec2) {
   float sum{0};
   for (uint8_t i{0}; i < vector_size; i++) {
@@ -341,7 +379,7 @@ float Matrix<rows, columns>::dotProduct(const Matrix<1, vector_size> &vec1,
 
 template <uint8_t rows, uint8_t columns>
 template <uint8_t vector_size>
-float Matrix<rows, columns>::dotProduct(const Matrix<vector_size, 1> &vec1,
+float Matrix<rows, columns>::DotProduct(const Matrix<vector_size, 1> &vec1,
                                         const Matrix<vector_size, 1> &vec2) {
   float sum{0};
   for (uint8_t i{0}; i < vector_size; i++) {
@@ -353,7 +391,11 @@ float Matrix<rows, columns>::dotProduct(const Matrix<vector_size, 1> &vec1,
 
 template <uint8_t rows, uint8_t columns>
 void Matrix<rows, columns>::Fill(float value) {
-  this->matrix.fill(value);
+  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
+    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
+      this->matrix[row_idx * columns + column_idx] = value;
+    }
+  }
 }
 
 template <uint8_t rows, uint8_t columns>
@@ -409,8 +451,8 @@ Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns> &
+float Matrix<rows, columns>::EuclideanNorm() const {
-Matrix<rows, columns>::Normalize(Matrix<rows, columns> &result) const {
+
   float sum{0};
   for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
     for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
@@ -419,90 +461,147 @@ Matrix<rows, columns>::Normalize(Matrix<rows, columns> &result) const {
     }
   }
 
-  if (sum == 0) {
+  return sqrt(sum);
-    // this wouldn't do anything anyways
-    result.Fill(1e+6);
-    return result;
-  }
-
-  sum = sqrt(sum);
-
-  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
-    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
-      result[row_idx][column_idx] = this->Get(row_idx, column_idx) / sum;
-    }
-  }
-
-  return result;
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, rows> Matrix<rows, columns>::Eye() {
+template <uint8_t sub_rows, uint8_t sub_columns, uint8_t row_offset,
-  Matrix<rows, rows> i_matrix;
+          uint8_t column_offset>
-  i_matrix.Fill(0);
+Matrix<sub_rows, sub_columns> Matrix<rows, columns>::SubMatrix() const {
-  for (uint8_t i{0}; i < rows; i++) {
+  // static assert that sub_rows + row_offset <= rows
-    i_matrix[i][i] = 1;
+  // static assert that sub_columns + column_offset <= columns
+  static_assert(sub_rows + row_offset <= rows,
+                "The submatrix you're trying to get is out of bounds (rows)");
+  static_assert(
+      sub_columns + column_offset <= columns,
+      "The submatrix you're trying to get is out of bounds (columns)");
+
+  Matrix<sub_rows, sub_columns> buffer{};
+  for (uint8_t row_idx{0}; row_idx < sub_rows; row_idx++) {
+    for (uint8_t column_idx{0}; column_idx < sub_columns; column_idx++) {
+      buffer[row_idx][column_idx] =
+          this->Get(row_idx + row_offset, column_idx + column_offset);
     }
-  return i_matrix;
+  }
+  return buffer;
 }
 
 template <uint8_t rows, uint8_t columns>
-void Matrix<rows, columns>::QR_Decomposition(Matrix<rows, columns> &Q,
+template <uint8_t sub_rows, uint8_t sub_columns>
-                                             Matrix<rows, columns> &R) const {
+void Matrix<rows, columns>::SetSubMatrix(
-  Q = Matrix<rows, columns>::Eye(); // Q starts as the identity matrix
+    uint8_t rowOffset, uint8_t columnOffset,
-  R = *this; // R starts as a copy of this matrix (For this algorithm we'll call
+    const Matrix<sub_rows, sub_columns> &sub_matrix) {
-             // this matrix A)
+  int16_t adjustedSubRows = sub_rows;
+  int16_t adjustedSubColumns = sub_columns;
+  int16_t adjustedRowOffset = rowOffset;
+  int16_t adjustedColumnOffset = columnOffset;
 
-  for (uint8_t row{0}; row < rows; row++) {
+  // a bunch of safety checks to make sure we don't overflow the matrix
-    // compute the householder vector
+  if (sub_rows > rows) {
-    const uint8_t houseHoldVectorSize{rows - row};
+    adjustedSubRows = rows;
-    const uint8_t subMatrixSize{columns - row};
+  }
-    Matrix<houseHoldVectorSize, 1> x{};
+  if (sub_columns > columns) {
-    this->SubMatrix(row, row, x);
+    adjustedSubColumns = columns;
+  }
 
-    Matrix<houseHoldVectorSize, 1> e1{};
+  if (adjustedSubRows + adjustedRowOffset >= rows) {
-    e1.Fill(0);
+    adjustedRowOffset =
-    if (x[0][0] >= 0) {
+        std::max(0, static_cast<int16_t>(rows) - adjustedSubRows);
-      e1[0][0] = x.Norm();
+  }
+
+  if (adjustedSubColumns + adjustedColumnOffset >= columns) {
+    adjustedColumnOffset =
+        std::max(0, static_cast<int16_t>(columns) - adjustedSubColumns);
+  }
+
+  for (uint8_t row_idx{0}; row_idx < adjustedSubRows; row_idx++) {
+    for (uint8_t column_idx{0}; column_idx < adjustedSubColumns; column_idx++) {
+      this->matrix[(row_idx + adjustedRowOffset) * columns + column_idx +
+                   adjustedColumnOffset] = sub_matrix.Get(row_idx, column_idx);
+    }
+  }
+}
+
+// QR decomposition: decomposes this matrix A into Q and R
+// Assumes square matrix
+template <uint8_t rows, uint8_t columns>
+void Matrix<rows, columns>::QRDecomposition(Matrix<rows, columns> &Q,
+                                            Matrix<columns, columns> &R) const {
+  static_assert(columns <= rows, "QR decomposition requires columns <= rows");
+
+  Q.Fill(0);
+  R.Fill(0);
+  Matrix<rows, 1> a_col, e, u, Q_column_k{};
+  Matrix<1, rows> e_T{};
+
+  for (uint8_t column = 0; column < columns; column++) {
+    this->GetColumn(column, a_col);
+    u = a_col;
+    // -----------------------
+    // ----- CALCULATE Q -----
+    // -----------------------
+    for (uint8_t k = 0; k <= column; k++) {
+      Q.GetColumn(k, Q_column_k);
+      Matrix<1, rows> Q_column_k_T = Q_column_k.Transpose();
+      u = u - Q_column_k * (Q_column_k_T * a_col);
+    }
+    float norm = u.EuclideanNorm();
+    if (norm > 1e-4) {
+      u = u / norm;
     } else {
-      e1[0][0] = -x.Norm();
+      u.Fill(0);
+    }
+    Q.SetSubMatrix(0, column, u);
+
+    // -----------------------
+    // ----- CALCULATE R -----
+    // -----------------------
+    for (uint8_t k = 0; k <= column; k++) {
+      Q.GetColumn(k, e);
+      R[k][column] = (a_col.Transpose() * e).Get(0, 0);
+    }
+  }
 }
 
-    Matrix<houseHoldVectorSize, 1> v = x + e1;
+template <uint8_t rows, uint8_t columns>
-    v = v * (1 / v.Norm()); // normalize V
+void Matrix<rows, columns>::EigenQR(Matrix<rows, rows> &eigenVectors,
+                                    Matrix<rows, 1> &eigenValues,
+                                    uint32_t maxIterations,
+                                    float tolerance) const {
+  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
-    // Apply the reflection to the R matrix
+  Matrix<rows, rows> QQ{Matrix<rows, rows>::Identity()};
-    // ************************************
+  Matrix<rows, rows> shift{0};
-    // initialize R's submatrix
-    Matrix<houseHoldVectorSize, subMatrixSize> R_subMatrix{};
-    R.SubMatrix(row, row, R_subMatrix);
-    // create some temporary buffers
-    Matrix<1, subMatrixSize> vR{};
-    Matrix<1, houseHoldVectorSize> v_T{};
-    v.Transpose(v_T);
-    Matrix<houseHoldVectorSize, subMatrixSize> vR_outer{};
-    // calculate the reflection
-    R_subMatrix =
-        R_subMatrix - 2 * Matrix<rows, columns>::OuterProduct(
-                              v_T, v_T.Mult(R_subMatrix, vR), vR_outer);
-    // save the reflection back to R
-    R.CopySubMatrixInto(row, row, R_subMatrix);
 
-    // ************************************
+  for (uint32_t iter = 0; iter < maxIterations; ++iter) {
-    // Apply the reflection to the Q matrix
+    Matrix<rows, rows> Q, R;
-    // ************************************
-    // initialize Q's submatrix
-    Matrix<rows, houseHoldVectorSize> Q_subMatrix{};
-    Q.SubMatrix(0, row, Q_subMatrix);
-    // create some temporary buffers
-    Matrix<rows, 1> Qv{};
-    Matrix<rows, houseHoldVectorSize> Qv_outer{};
 
-    Q_subMatrix = Q_subMatrix - 2 * Matrix<rows, columns>::OuterProduct(
+    // // QR shift lets us "attack" the first diagonal to speed up the algorithm
-                                        Q_subMatrix.Mult(v, Qv), v, Qv_outer);
+    // shift = Matrix<rows, rows>::Identity() * Ak[rows - 1][rows - 1];
-    Q.CopySubMatrixInto(0, row, Q_subMatrix);
+    (Ak - shift).QRDecomposition(Q, R);
+    Ak = R * Q + shift;
+    QQ = QQ * Q;
+
+    // Check convergence: off-diagonal norm
+    float offDiagSum = 0.0f;
+    for (uint32_t row = 1; row < rows; row++) {
+      for (uint32_t column = 0; column < row; column++) {
+        offDiagSum += fabs(Ak[row][column]);
       }
     }
 
+    if (offDiagSum < tolerance) {
+      break;
+    }
+  }
+
+  // Diagonal elements are the eigenvalues
+  for (uint8_t i = 0; i < rows; i++) {
+    eigenValues[i][0] = Ak[i][i];
+  }
+  eigenVectors = QQ;
+}
+
 #endif // MATRIX_H_

src/Matrix.hpp

@@ -1,8 +1,9 @@
-#ifndef MATRIX_H_
+#pragma once
-#define MATRIX_H_
 
 #include <array>
 #include <cstdint>
+#include <string>
+#include <type_traits>
 
 // TODO: Add a function to calculate eigenvalues/vectors
 // TODO: Add a function to compute RREF
@@ -11,16 +12,13 @@
 
 template <uint8_t rows, uint8_t columns> class Matrix {
 public:
+  static_assert(rows > 0, "Template error: rows must be greater than 0.");
+  static_assert(columns > 0, "Template error: columns must be greater than 0.");
   /**
    * @brief create a matrix but leave all of its values unitialized
    */
   Matrix() = default;
 
-  /**
-   * @brief Create a matrix but fill all of its entries with one value
-   */
-  Matrix(float value);
-
   /**
    * @brief Initialize a matrix with an array
    */
@@ -32,9 +30,17 @@ public:
   Matrix(const Matrix<rows, columns> &other);
 
   /**
-   * @brief Initialize a matrix directly with any number of arguments
+   * @brief Initialize a matrix directly with scalar values
+   * Uses SFINAE to only accept arithmetic types (int, float, double, etc.)
    */
-  template <typename... Args> Matrix(Args... args);
+  template <typename... Args,
+            std::enable_if_t<(std::is_arithmetic_v<Args> && ...), int> = 0>
+  Matrix(Args... args);
+
+  /**
+   * @brief Create an identity matrix
+   */
+  static Matrix<rows, columns> Identity();
 
   /**
    * @brief Set all elements in this to value
@@ -112,79 +118,20 @@ public:
    * @param result A buffer to store the result into
    * @warning this is super slow! Only call it if you absolutely have to!!!
    */
-  Matrix<rows, columns> &Invert(Matrix<rows, columns> &result) const;
+  Matrix<rows, columns> Invert() const;
 
   /**
    * @brief Transpose this matrix
    * @param result A buffer to store the result into
    */
-  Matrix<columns, rows> &Transpose(Matrix<columns, rows> &result) const;
+  Matrix<columns, rows> Transpose() const;
 
   /**
-   * @brief reduce the matrix so the sum of its elements equal 1
+   * @brief Returns the euclidean magnitude of the matrix. Also known as the L2
+   * norm
    * @param result a buffer to store the result into
    */
-  Matrix<rows, columns> &Normalize(Matrix<rows, columns> &result) const;
+  float EuclideanNorm() const;
-
-  /**
-   * @brief return an identity matrix of the specified size
-   */
-  static Matrix<rows, rows> Eye();
-
-  /**
-   * @brief write a copy of a sub matrix into the given result matrix.
-   * @param rowIndex The row index to start the copy from
-   * @param columnIndex the column index to start the copy from
-   * @param result the matrix buffer to write the sub matrix into. The size of
-   * the matrix buffer allows the function to determine the end indices of the
-   * sub matrix
-   */
-  template <uint8_t subRows, uint8_t subColumns>
-  Matrix<subRows, subColumns> &
-  SubMatrix(uint8_t rowIndex, uint8_t columnIndex,
-            Matrix<subRows, subColumns> &result) const {
-    return result;
-  }
-
-  /**
-   * @brief write a copy of a sub matrix into this matrix starting at the given
-   * idnex.
-   * @param rowIndex The row index to start the copy from
-   * @param columnIndex the column index to start the copy from
-   * @param subMatrix The submatrix to copy into this matrix. The size of
-   * the matrix buffer allows the function to determine the end indices of the
-   * sub matrix
-   */
-  template <uint8_t subRows, uint8_t subColumns>
-  void CopySubMatrixInto(uint8_t rowIndex, uint8_t columnIndex,
-                         const Matrix<subRows, subColumns> &subMatrix) {}
-
-  /**
-   * @brief Returns the norm of the matrix
-   */
-  float Norm() { return 0; }
-
-  template <uint8_t vec1Length, uint8_t vec2Length>
-  static Matrix<vec1Length, vec2Length> &
-  OuterProduct(const Matrix<1, vec1Length> &vec1,
-               const Matrix<1, vec2Length> &vec2,
-               Matrix<vec1Length, vec2Length> &result) {
-    return result;
-  }
-
-  template <uint8_t vec1Length, uint8_t vec2Length>
-  static Matrix<vec1Length, vec2Length> &
-  OuterProduct(const Matrix<vec1Length, 1> &vec1,
-               const Matrix<vec2Length, 1> &vec2,
-               Matrix<vec1Length, vec2Length> &result) {
-    return result;
-  }
-  /**
-   * @brief Calulcate the QR decomposition of a matrix
-   * @param Q the
-   */
-  void QR_Decomposition(Matrix<rows, columns> &Q,
-                        Matrix<rows, columns> &R) const;
 
   /**
    * @brief Get a row from the matrix
@@ -211,8 +158,16 @@ public:
    */
   constexpr uint8_t GetColumnSize() { return columns; }
 
+  /**
+   * @brief Write a string representation of the matrix into the buffer
+   */
   void ToString(std::string &stringBuffer) const;
 
+  /**
+   * @brief Returns the internal representation of the matrix as an array
+   */
+  const float *ToArray() const;
+
   /**
    * @brief Get an element from the matrix
    * @param row the row index of the element
@@ -221,10 +176,6 @@ public:
    */
   float Get(uint8_t row_index, uint8_t column_index) const;
 
-  // *******************************************************
-  // ************** OPERATOR OVERRIDES *********************
-  // *******************************************************
-
   /**
    * @brief get the specified row of the matrix returned as a reference to the
    * internal array
@@ -243,29 +194,68 @@ public:
 
   Matrix<rows, columns> operator-(const Matrix<rows, columns> &other) const;
 
-  Matrix<rows, columns> operator*(const Matrix<rows, columns> &other) const;
+  template <uint8_t other_columns>
+  Matrix<rows, other_columns>
+  operator*(const Matrix<columns, other_columns> &other) const;
 
   Matrix<rows, columns> operator*(float scalar) const;
 
-private:
+  Matrix<rows, columns> operator/(float scalar) const;
+
+  template <uint8_t sub_rows, uint8_t sub_columns, uint8_t row_offset,
+            uint8_t column_offset>
+  Matrix<sub_rows, sub_columns> SubMatrix() const;
+
+  template <uint8_t sub_rows, uint8_t sub_columns>
+  void SetSubMatrix(uint8_t rowOffset, uint8_t columnOffset,
+                    const Matrix<sub_rows, sub_columns> &sub_matrix);
+
   /**
    * @brief take the dot product of the two vectors
    */
   template <uint8_t vector_size>
-  static float dotProduct(const Matrix<1, vector_size> &vec1,
+  static float DotProduct(const Matrix<1, vector_size> &vec1,
                           const Matrix<1, vector_size> &vec2);
 
   template <uint8_t vector_size>
-  static float dotProduct(const Matrix<vector_size, 1> &vec1,
+  static float DotProduct(const Matrix<vector_size, 1> &vec1,
                           const Matrix<vector_size, 1> &vec2);
 
+  static float DotProduct(const Matrix<1, 1> &vec1, const Matrix<1, 1> &vec2) {
+    return vec1.Get(0, 0) * vec2.Get(0, 0);
+  }
+
+  /**
+   * @brief Performs QR decomposition on this matrix
+   * @param Q a buffer that will contain Q after the function completes
+   * @param R a buffer that will contain R after the function completes
+   */
+  void QRDecomposition(Matrix<rows, columns> &Q,
+                       Matrix<columns, columns> &R) const;
+
+  /**
+   * @brief Uses QR decomposition to efficiently calculate the eigenvectors
+   * and values of this matrix
+   * @param eigenVectors a buffer that will contain the eigenvectors fo this
+   * matrix
+   * @param eigenValues a buffer that will contain the eigenValues fo this
+   * matrix
+   * @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.
+   */
+  void EigenQR(Matrix<rows, rows> &eigenVectors, Matrix<rows, 1> &eigenValues,
+               uint32_t maxIterations = 1000, float tolerance = 1e-6f) const;
+
+protected:
+  std::array<float, rows * columns> matrix;
+
+private:
   Matrix<rows, columns> &adjugate(Matrix<rows, columns> &result) const;
 
   void setMatrixToArray(const std::array<float, rows * columns> &array);
-
-  std::array<float, rows * columns> matrix;
 };
 
+#ifndef MATRIX_H_
 #include "Matrix.cpp"
-
 #endif // MATRIX_H_

src/Quaternion.cpp

@@ -0,0 +1,120 @@
+#include "Quaternion.h"
+#include <cmath>
+
+/**
+ * @brief Create a quaternion from an angle and axis
+ * @param angle The angle to rotate by
+ * @param axis The axis to rotate around
+ */
+Quaternion Quaternion::FromAngleAndAxis(float angle, const Matrix<1, 3> &axis) {
+  const float halfAngle = angle / 2;
+  const float sinHalfAngle = sin(halfAngle);
+  Matrix<1, 3> normalizedAxis = axis / axis.EuclideanNorm();
+  return Quaternion{static_cast<float>(cos(halfAngle)),
+                    normalizedAxis.Get(0, 0) * sinHalfAngle,
+                    normalizedAxis.Get(0, 1) * sinHalfAngle,
+                    normalizedAxis.Get(0, 2) * sinHalfAngle};
+}
+
+float Quaternion::operator[](uint8_t index) const {
+  if (index < 4) {
+    return this->matrix[index];
+  }
+
+  // index out of bounds
+  return 1e+6;
+}
+
+void Quaternion::operator=(const Quaternion &other) {
+  memcpy(&(this->matrix), &(other.matrix), 4 * sizeof(float));
+}
+
+Quaternion Quaternion::operator*(const Quaternion &other) const {
+  Quaternion result{};
+  this->Q_Mult(other, result);
+  return result;
+}
+
+Quaternion Quaternion::operator*(float scalar) const {
+  return Quaternion{this->w * scalar, this->v1 * scalar, this->v2 * scalar,
+                    this->v3 * scalar};
+}
+
+Quaternion Quaternion::operator+(const Quaternion &other) const {
+  return Quaternion{this->w + other.w, this->v1 + other.v1, this->v2 + other.v2,
+                    this->v3 + other.v3};
+}
+
+Quaternion &Quaternion::Q_Mult(const Quaternion &other,
+                               Quaternion &buffer) const {
+
+  // eq. 6
+  buffer.w = (other.w * this->w - other.v1 * this->v1 - other.v2 * this->v2 -
+              other.v3 * this->v3);
+  buffer.v1 = (other.w * this->v1 + other.v1 * this->w - other.v2 * this->v3 +
+               other.v3 * this->v2);
+  buffer.v2 = (other.w * this->v2 + other.v1 * this->v3 + other.v2 * this->w -
+               other.v3 * this->v1);
+  buffer.v3 = (other.w * this->v3 - other.v1 * this->v2 + other.v2 * this->v1 +
+               other.v3 * this->w);
+  return buffer;
+}
+
+Quaternion &Quaternion::Rotate(Quaternion &other, Quaternion &buffer) const {
+  Quaternion prime{this->w, -this->v1, -this->v2, -this->v3};
+  buffer.v1 = other.v1;
+  buffer.v2 = other.v2;
+  buffer.v3 = other.v3;
+  buffer.w = 0;
+
+  Quaternion temp{};
+  this->Q_Mult(buffer, temp);
+  temp.Q_Mult(prime, buffer);
+  return buffer;
+}
+
+void Quaternion::Normalize() {
+  float magnitude = sqrt(this->v1 * this->v1 + this->v2 * this->v2 +
+                         this->v3 * this->v3 + this->w * this->w);
+  if (magnitude == 0) {
+    return;
+  }
+  this->v1 /= magnitude;
+  this->v2 /= magnitude;
+  this->v3 /= magnitude;
+  this->w /= magnitude;
+}
+
+Matrix<3, 3> Quaternion::ToRotationMatrix() const {
+  float xx = this->v1 * this->v1;
+  float yy = this->v2 * this->v2;
+  float zz = this->v3 * this->v3;
+  Matrix<3, 3> rotationMatrix{1 - 2 * (yy - zz),
+                              2 * (this->v1 * this->v2 - this->v3 * this->w),
+                              2 * (this->v1 * this->v3 + this->v2 * this->w),
+                              2 * (this->v1 * this->v2 + this->v3 * this->w),
+                              1 - 2 * (xx - zz),
+                              2 * (this->v2 * this->v3 - this->v1 * this->w),
+                              2 * (this->v1 * this->v3 - this->v2 * this->w),
+                              2 * (this->v2 * this->v3 + this->v1 * this->w),
+                              1 - 2 * (xx - yy)};
+  return rotationMatrix;
+};
+
+Matrix<3, 1> Quaternion::ToEulerAngle() const {
+  float sqv1 = this->v1 * this->v1;
+  float sqv2 = this->v2 * this->v2;
+  float sqv3 = this->v3 * this->v3;
+  float sqw = this->w * this->w;
+
+  Matrix<3, 1> eulerAngle;
+  {
+    atan2(2.0 * (this->v1 * this->v2 + this->v3 * this->w),
+          (sqv1 - sqv2 - sqv3 + sqw));
+    asin(-2.0 * (this->v1 * this->v3 - this->v2 * this->w) /
+         (sqv1 + sqv2 + sqv3 + sqw));
+    atan2(2.0 * (this->v2 * this->v3 + this->v1 * this->w),
+          (-sqv1 - sqv2 + sqv3 + sqw));
+  };
+  return eulerAngle;
+}

src/Quaternion.h

@@ -0,0 +1,90 @@
+#ifndef QUATERNION_H_
+#define QUATERNION_H_
+
+#include "Matrix.hpp"
+class Quaternion : public Matrix<1, 4> {
+public:
+  Quaternion() : Matrix<1, 4>() {}
+  Quaternion(float w, float v1, float v2, float v3)
+      : Matrix<1, 4>(w, v1, v2, v3) {}
+  Quaternion(const Quaternion &q) : Matrix<1, 4>(q.w, q.v1, q.v2, q.v3) {}
+  Quaternion(const Matrix<1, 4> &matrix) : Matrix<1, 4>(matrix) {}
+  Quaternion(const std::array<float, 4> &array) : Matrix<1, 4>(array) {}
+
+  /**
+   * @brief Create a quaternion from an angle and axis
+   * @param angle The angle to rotate by
+   * @param axis The axis to rotate around
+   */
+  static Quaternion FromAngleAndAxis(float angle, const Matrix<1, 3> &axis);
+
+  /**
+   * @brief Access the elements of the quaternion
+   * @param index The index of the element to access
+   * @return The value of the element at the index
+   */
+  float operator[](uint8_t index) const;
+
+  /**
+   * @brief Assign one quaternion to another
+   */
+  void operator=(const Quaternion &other);
+
+  /**
+   * @brief Do quaternion multiplication
+   */
+  Quaternion operator*(const Quaternion &other) const;
+
+  /**
+   * @brief Multiply the quaternion by a scalar
+   */
+  Quaternion operator*(float scalar) const;
+
+  /**
+   * @brief Add two quaternions together
+   * @param other The quaternion to add to this one
+   * @return The net quaternion
+   */
+  Quaternion operator+(const Quaternion &other) const;
+
+  /**
+   * @brief Q_Mult a quaternion by another quaternion
+   * @param other The quaternion to rotate by
+   * @param buffer The buffer to store the result in
+   * @return A reference to the buffer
+   */
+  Quaternion &Q_Mult(const Quaternion &other, Quaternion &buffer) const;
+
+  /**
+   * @brief Rotate a quaternion by this quaternion
+   * @param other The quaternion to rotate
+   * @param buffer The buffer to store the result in
+   *
+   */
+  Quaternion &Rotate(Quaternion &other, Quaternion &buffer) const;
+
+  /**
+   * @brief Normalize the quaternion to a magnitude of 1
+   */
+  void Normalize();
+
+  /**
+   * @brief Convert the quaternion to a rotation matrix
+   * @return The rotation matrix
+   */
+  Matrix<3, 3> ToRotationMatrix() const;
+
+  /**
+   * @brief Convert the quaternion to an Euler angle representation
+   * @return The Euler angle representation of the quaternion
+   */
+  Matrix<3, 1> ToEulerAngle() const;
+
+  // Give people an easy way to access the elements
+  float &w{matrix[0]};
+  float &v1{matrix[1]};
+  float &v2{matrix[2]};
+  float &v3{matrix[3]};
+};
+
+#endif // QUATERNION_H_

src/SVD.cpp

@@ -0,0 +1,846 @@
+// This #ifndef section makes clangd happy so that it can properly do type hints
+// in this file
+#ifndef SVD_H_
+#define SVD_H_
+#include "SVD.hpp"
+#endif
+
+#ifdef SVD_H_ // since the .cpp file has to be included by the .hpp file this
+              // will evaluate to true
+#include "SVD.hpp"
+#include <cstdint>
+
+// ============================================================================
+// SVD Building Block Implementations
+// ============================================================================
+
+float SVD::ComputeHouseholder(const float *x, uint8_t len, float *v,
+                              float &alpha) {
+  // Compute ||x||
+  float norm = 0.0f;
+  for (uint8_t i = 0; i < len; i++) {
+    norm += x[i] * x[i];
+  }
+  norm = sqrtf(norm);
+
+  if (norm < 1e-30f) {
+    alpha = 0.0f;
+    for (uint8_t i = 0; i < len; i++) {
+      v[i] = 0.0f;
+    }
+    return 0.0f;
+  }
+
+  // Choose sign to avoid cancellation: alpha has opposite sign of x[0]
+  alpha = (x[0] >= 0.0f) ? -norm : norm;
+
+  // v = x - alpha * e1, then normalize
+  float v0 = x[0] - alpha;
+
+  // Compute ||v||² directly: v0² + x₁² + ... + xₙ₋₁²
+  float vv = v0 * v0;
+  for (uint8_t i = 1; i < len; i++) {
+    vv += x[i] * x[i];
+  }
+
+  if (vv < 1e-30f) {
+    // Already aligned with e1
+    for (uint8_t i = 0; i < len; i++) {
+      v[i] = (i == 0) ? 1.0f : 0.0f;
+    }
+    return norm;
+  }
+
+  float scale = 1.0f / sqrtf(vv);
+  for (uint8_t i = 0; i < len; i++) {
+    v[i] = (i == 0) ? v0 * scale : x[i] * scale;
+  }
+
+  return norm;
+}
+
+void SVD::ApplyHouseholderLeft(Matrix<5, 5> &W, const float *v,
+                               uint8_t startRow, uint8_t endRow) {
+  uint8_t len = endRow - startRow + 1;
+
+  // Compute vᵀv (should be 2.0 for our normalized vectors, but compute
+  // explicitly)
+  float vv = 0.0f;
+  for (uint8_t i = 0; i < len; i++) {
+    vv += v[i] * v[i];
+  }
+  if (vv < 1e-30f)
+    return;
+
+  float twoOverVv = 2.0f / vv;
+
+  // W = (I - 2vvᵀ) · W
+  for (uint8_t col = 0; col < 5; col++) {
+    float dot = 0.0f;
+    for (uint8_t i = 0; i < len; i++) {
+      dot += v[i] * W[startRow + i][col];
+    }
+    dot *= twoOverVv;
+    for (uint8_t i = 0; i < len; i++) {
+      W[startRow + i][col] -= dot * v[i];
+    }
+  }
+}
+
+void SVD::ApplyHouseholderRight(Matrix<5, 5> &W, const float *v,
+                                uint8_t startCol, uint8_t endCol) {
+  uint8_t len = endCol - startCol + 1;
+
+  float vv = 0.0f;
+  for (uint8_t i = 0; i < len; i++) {
+    vv += v[i] * v[i];
+  }
+  if (vv < 1e-30f)
+    return;
+
+  float twoOverVv = 2.0f / vv;
+
+  // W = W · (I - 2vvᵀ)
+  for (uint8_t row = 0; row < 5; row++) {
+    float dot = 0.0f;
+    for (uint8_t i = 0; i < len; i++) {
+      dot += W[row][startCol + i] * v[i];
+    }
+    dot *= twoOverVv;
+    for (uint8_t i = 0; i < len; i++) {
+      W[row][startCol + i] -= dot * v[i];
+    }
+  }
+}
+
+[[gnu::unused]] void SVD::ComputeGivens(float x, float y, float &c, float &s) {
+  float r = sqrtf(x * x + y * y);
+
+  if (r < 1e-30f) {
+    c = 1.0f;
+    s = 0.0f;
+    return;
+  }
+
+  c = x / r;
+  s = y / r;
+}
+
+[[gnu::unused]] void SVD::ApplyGivensLeft(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
+                          float s, uint8_t startCol, uint8_t endCol) {
+  // [c  s] [row_i]   = [new_row_i]
+  // [-s c] [row_j]     [new_row_j]
+  for (uint8_t col = startCol; col <= endCol && col < 5; col++) {
+    float t1 = W[i][col];
+    float t2 = W[j][col];
+    W[i][col] = c * t1 + s * t2;
+    W[j][col] = -s * t1 + c * t2;
+  }
+}
+
+[[gnu::unused]] void SVD::ApplyGivensRight(Matrix<5, 5> &W, uint8_t i, uint8_t j, float c,
+                           float s, uint8_t startRow, uint8_t endRow) {
+  // [col_i  col_j] · [c -s] = [new_col_i  new_col_j]
+  //                  [s  c]
+  for (uint8_t row = startRow; row <= endRow && row < 5; row++) {
+    float t1 = W[row][i];
+    float t2 = W[row][j];
+    W[row][i] = c * t1 + s * t2;
+    W[row][j] = -s * t1 + c * t2;
+  }
+}
+
+// ============================================================================
+// Phase 1: Householder Bidiagonalization
+// ============================================================================
+
+void SVD::Bidiagonalize(Matrix<5, 5> &W,
+                        uint8_t m, uint8_t q, uint8_t p,
+                        Matrix<5, 5> &QL,
+                        Matrix<5, 5> &QR) {
+  // Working matrix W is m×q (padded to 5×5).
+  // QL and QR are initialized to identity by the caller.
+  // We reduce W to upper bidiagonal form B using Householder reflections.
+
+  float hhVec[5]; // Householder vector storage
+
+  for (uint8_t k = 0; k < p; k++) {
+    // --- Left Householder on column k, rows k..m-1 ---
+    // Zero out subdiagonal elements below B[k+1][k]
+    {
+      uint8_t len = m - k;
+      if (len <= 1)
+        continue;
+
+      // Extract the column segment W[k..k+len-1][k]
+      float x[5];
+      for (uint8_t i = 0; i < len; i++) {
+        x[i] = W[k + i][k];
+      }
+
+      // Compute Householder reflector
+      float alpha;
+      SVD::ComputeHouseholder(x, len, hhVec, alpha);
+
+      if (alpha == 0.0f)
+        continue;
+
+      // Apply H from left to W: W = H·W (columns k..q-1)
+      SVD::ApplyHouseholderLeft(W, hhVec, k, k + len - 1);
+
+      // Apply H from right to QL: QL = QL · H
+      SVD::ApplyHouseholderRight(QL, hhVec, k, k + len - 1);
+    }
+
+    // --- Right Householder on row k, columns k+1..q-1 ---
+    // Zero out elements above the first superdiagonal in row k.
+    // The Householder maps [W[k][k+1], ..., W[k][q-1]] to [gamma, 0, ..., 0],
+    // preserving the first superdiagonal element (now gamma) and zeroing the rest.
+    {
+      int len = static_cast<int>(q) - 1 - k;
+      if (len <= 1)
+        continue; // Need at least 2 elements to zero something out
+
+      // Extract the row segment starting from column k+1
+      float x[5];
+      for (uint8_t i = 0; i < len; i++) {
+        x[i] = W[k][k + 1 + i];
+      }
+
+      // Compute Householder reflector
+      float alpha;
+      SVD::ComputeHouseholder(x, len, hhVec, alpha);
+
+      if (alpha == 0.0f)
+        continue;
+
+      // Apply H from right to W: W = W·H (columns k+1..k+len-1)
+      SVD::ApplyHouseholderRight(W, hhVec, k + 1, k + len);
+
+      // Apply H from right to QR: QR = QR · H
+      SVD::ApplyHouseholderRight(QR, hhVec, k + 1, k + len);
+    }
+  }
+}
+
+// ============================================================================
+// Phase 2 helpers: block solving of the bidiagonal matrix
+// ============================================================================
+
+void SVD::DeflateBidiagonal(Matrix<5, 5> &W, uint8_t p, float tol) {
+  // Zero out superdiagonal elements that are negligible relative to the
+  // local diagonal scale. This deflates the bidiagonal matrix into
+  // independent unreduced blocks, each of which can be solved on its own.
+  if (p < 2)
+    return;
+  for (uint8_t i = 0; i < p - 1; i++) {
+    float test = fabsf(W[i][i + 1]);
+    float scale = fabsf(W[i][i]) + fabsf(W[i + 1][i + 1]);
+    // Use absolute threshold for small scales to avoid division issues
+    if (test < tol * fmaxf(scale, 1e-10f)) {
+      W[i][i + 1] = 0;
+    }
+  }
+}
+
+bool SVD::BidiagonalIsDiagonal(const Matrix<5, 5> &W, uint8_t p, float tol) {
+  // True when every superdiagonal element of the p×p bidiagonal matrix
+  // has been reduced to (numerically) zero, i.e. the diagonal holds the
+  // singular values and no unreduced blocks remain.
+  if (p < 2)
+    return true;
+  for (uint8_t i = 0; i < p - 1; i++) {
+    if (fabsf(W.Get(i, i + 1)) > tol * 1e-30f) {
+      return false;
+    }
+  }
+  return true;
+}
+
+void SVD::SolveBidiagonalBlock2x2(float a, float b, float d, float Ublock[2][2],
+                                  float Vblock[2][2], float sigma[2]) {
+  // Full SVD of the 2×2 upper-bidiagonal block B = [[a, b], [0, d]]:
+  //   B = Ublock · diag(sigma[0], sigma[1]) · Vblockᵀ
+  // where:
+  //   - sigma[0] ≥ sigma[1] ≥ 0
+  //   - columns of Ublock are the left singular vectors
+  //   - columns of Vblock are the right singular vectors (Vblock = scipy Vᵀᵀ)
+  //
+  // Uses eigen-decomposition of BᵀB = [[a², ab], [ab, b²+d²]] (symmetric
+  // 2×2, closed form), then uᵢ = B·vᵢ/σᵢ.
+
+  // Singular values = sqrt of eigenvalues of BᵀB (trace/det closed form)
+  float trace = a * a + b * b + d * d;
+  float det = a * a * d * d;
+  float disc = trace * trace - 4.0f * det;
+  if (disc < 0)
+    disc = 0;
+  float sqrtDisc = sqrtf(disc);
+  float hi = sqrtf((trace + sqrtDisc) / 2.0f);
+  float lo = sqrtf((trace - sqrtDisc) / 2.0f);
+  if (lo > hi) {
+    float tmp = hi;
+    hi = lo;
+    lo = tmp;
+  }
+  sigma[0] = hi;
+  sigma[1] = lo;
+
+  // Right singular vector v1: eigenvector of BᵀB for λ1 = hi².
+  // Null-space vector of (BᵀB − λ1·I) is [ab, λ1 − a²].
+  float a2 = a * a;
+  float ab_val = a * b;
+  float e1x = ab_val;
+  float e1y = hi * hi - a2;
+  float normE1 = sqrtf(e1x * e1x + e1y * e1y);
+  float v1x, v1y;
+  if (normE1 > 1e-30f) {
+    v1x = e1x / normE1;
+    v1y = e1y / normE1;
+  } else {
+    // Degenerate (e.g. b = 0 and |a| ≥ |d|): e₁ is already an eigenvector
+    v1x = 1.0f;
+    v1y = 0.0f;
+  }
+
+  // v2 is the unit vector orthogonal to v1 (completes the 2D basis)
+  float v2x = -v1y;
+  float v2y = v1x;
+
+  // Vblock columns = right singular vectors
+  Vblock[0][0] = v1x;
+  Vblock[1][0] = v1y;
+  Vblock[0][1] = v2x;
+  Vblock[1][1] = v2y;
+
+  // Ublock columns: uᵢ = B·vᵢ / σᵢ, with a rank-deficiency guard.
+  // When σᵢ ≈ 0, dividing produces inf/NaN; instead fill the U column with
+  // the signed orthogonal complement of the other U column (keeps Ublock
+  // orthogonal, and B·vᵢ ≈ 0 so any unit complement satisfies the SVD).
+  float u1x, u1y, u2x, u2y;
+  if (hi > 1e-30f) {
+    u1x = (a * v1x + b * v1y) / hi;
+    u1y = d * v1y / hi;
+  } else {
+    u1x = 1.0f;
+    u1y = 0.0f;
+  }
+  if (lo > 1e-30f) {
+    u2x = (a * v2x + b * v2y) / lo;
+    u2y = d * v2y / lo;
+  } else {
+    u2x = -u1y;
+    u2y = u1x;
+  }
+
+  Ublock[0][0] = u1x;
+  Ublock[1][0] = u1y;
+  Ublock[0][1] = u2x;
+  Ublock[1][1] = u2y;
+}
+
+void SVD::JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5],
+                               float V[5][5]) {
+  // Cyclic Jacobi eigenvalue algorithm on symmetric n×n matrix T (in place).
+  // On return:
+  //   - T is (near-)diagonal; its diagonal entries are the eigenvalues
+  //   - evals[i] = T[i][i] (unsorted)
+  //   - columns of V are the corresponding eigenvectors (V is accumulated
+  //     as V ← V·J so that T·V = V·Λ)
+  float jacTol = 1e-10f;
+
+  // V starts as the identity: eigenvector accumulator
+  for (uint8_t i = 0; i < n; i++)
+    for (uint8_t j = 0; j < n; j++)
+      V[i][j] = (i == j) ? 1.0f : 0.0f;
+
+  for (uint32_t jacIter = 0; jacIter < 100; jacIter++) {
+    // Check convergence over ALL off-diagonal entries, not just the
+    // tridiagonal band: cyclic Jacobi on a 3x3+ block fills non-band
+    // entries (e.g. T[0][2]) during sweeps, so a band-only test can
+    // declare convergence too early.
+    bool converged = true;
+    for (uint8_t i = 0; i < n - 1 && converged; i++) {
+      for (uint8_t j = i + 1; j < n; j++) {
+        float scale = fabsf(T[i][i]) + fabsf(T[j][j]);
+        if (fabsf(T[i][j]) > jacTol * fmaxf(scale, 1e-30f)) {
+          converged = false;
+          break;
+        }
+      }
+    }
+    if (converged)
+      break;
+
+    // Cyclic Jacobi: zero out T[p][q] for p < q
+    for (uint8_t p = 0; p < n - 1; p++) {
+      for (uint8_t q = p + 1; q < n; q++) {
+        float tPQ = T[p][q];
+        if (fabsf(tPQ) < jacTol * 1e-30f)
+          continue;
+
+        float tPP = T[p][p];
+        float tQQ = T[q][q];
+        float theta = (tQQ - tPP) / (2.0f * tPQ);
+        float t;
+        if (theta >= 0.0f)
+          t = 1.0f / (theta + sqrtf(1.0f + theta * theta));
+        else
+          t = -1.0f / (-theta + sqrtf(1.0f + theta * theta));
+
+        float c = 1.0f / sqrtf(1.0f + t * t);
+        float s = t * c;
+
+        // Update T
+        T[p][p] = tPP - t * tPQ;
+        T[q][q] = tQQ + t * tPQ;
+        T[p][q] = 0.0f;
+        T[q][p] = 0.0f;
+
+        // Update other elements
+        for (uint8_t k = 0; k < n; k++) {
+          if (k == p || k == q)
+            continue;
+          float tPK = T[k][p];
+          float tQK = T[k][q];
+          T[k][p] = c * tPK - s * tQK;
+          T[p][k] = T[k][p];
+          T[k][q] = s * tPK + c * tQK;
+          T[q][k] = T[k][q];
+        }
+
+        // Accumulate eigenvectors
+        for (uint8_t k = 0; k < n; k++) {
+          float vKP = V[k][p];
+          float vKQ = V[k][q];
+          V[k][p] = c * vKP - s * vKQ;
+          V[k][q] = s * vKP + c * vKQ;
+        }
+      }
+    }
+  }
+
+  for (uint8_t i = 0; i < n; i++) {
+    evals[i] = fabsf(T[i][i]);
+  }
+}
+
+void SVD::ApplyBlockFactorsToAccumulators(uint8_t blockStart, uint8_t blockSize,
+                                          const float Ublock[5][5],
+                                          const float Vblock[5][5],
+                                          uint8_t rowsQL, uint8_t rowsQR,
+                                          Matrix<5, 5> &QL,
+                                          Matrix<5, 5> &QR) {
+  // Fold the block SVD factors into the accumulated Householder
+  // transformation matrices:
+  //   QL[:, blockStart..blockStart+blockSize-1] ← QL[:, ...] · Ublock
+  //     (over rows 0..rowsQL−1)
+  //   QR[:, blockStart..blockStart+blockSize-1] ← QR[:, ...] · Vblock
+  //     (over rows 0..rowsQR−1)
+  //
+  // rowsQL / rowsQR are the meaningful row extents of the accumulators:
+  // for a transposed (wide) problem W = Aᵀ has n rows, so QL carries n
+  // meaningful rows while in the normal case it carries m.
+
+  for (uint8_t j = 0; j < rowsQL; j++) {
+    for (uint8_t i = 0; i < blockSize; i++) {
+      float sum = 0.0f;
+      for (uint8_t k = 0; k < blockSize; k++) {
+        sum += QL[j][blockStart + k] * Ublock[k][i];
+      }
+      QL[j][blockStart + i] = sum;
+    }
+  }
+
+  for (uint8_t j = 0; j < rowsQR; j++) {
+    for (uint8_t i = 0; i < blockSize; i++) {
+      float sum = 0.0f;
+      for (uint8_t k = 0; k < blockSize; k++) {
+        sum += QR[j][blockStart + k] * Vblock[k][i];
+      }
+      QR[j][blockStart + i] = sum;
+    }
+  }
+}
+
+void SVD::SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart,
+                                     uint8_t blockSize, uint8_t rowsQL,
+                                     uint8_t rowsQR, Matrix<5, 5> &QL,
+                                     Matrix<5, 5> &QR, float tol) {
+  // Full SVD of an unreduced upper-bidiagonal block of size > 2 via
+  // eigen-decomposition of the tridiagonal T = BᵀB:
+  //   1. Snapshot the ORIGINAL block diagonal/superdiagonal from W
+  //   2. Form T = BᵀB (tridiagonal symmetric)
+  //   3. JacobiEigenSymmetric → eigenvalues + eigenvector matrix V
+  //   4. Sort eigenvalues descending, reordering V
+  //   5. Ublock = B_orig · V · Σ⁻¹ (computed from the SNAPSHOT so that
+  //      overwriting W's diagonal does not corrupt it)
+  //   6. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators
+  //   7. Only now write sqrt(eigenvalues) into W's diagonal and zero the
+  //      block's superdiagonals
+  (void)tol; // Jacobi convergence tolerance is internal
+
+  // Step 1: snapshot original block values (diagonal d[i], superdiag e[i])
+  float d[5], e[4];
+  for (uint8_t i = 0; i < blockSize; i++) {
+    d[i] = W[blockStart + i][blockStart + i];
+  }
+  for (uint8_t i = 0; i < blockSize - 1; i++) {
+    e[i] = W[blockStart + i][blockStart + i + 1];
+  }
+
+  // Step 2: form T = BᵀB (tridiagonal)
+  //   T[i][i]   = d[i]² + e[i−1]²  (e[−1] = 0)
+  //   T[i][i+1] = d[i] · e[i]
+  float T[5][5] = {{0}};
+  for (uint8_t i = 0; i < blockSize; i++) {
+    float diag = d[i] * d[i];
+    if (i > 0) {
+      diag += e[i - 1] * e[i - 1];
+    }
+    T[i][i] = diag;
+    if (i < blockSize - 1) {
+      float off = d[i] * e[i];
+      T[i][i + 1] = off;
+      T[i + 1][i] = off;
+    }
+  }
+
+  // Step 3: Jacobi eigenvalue algorithm
+  float evals[5] = {0};
+  float V[5][5] = {{0}};
+  SVD::JacobiEigenSymmetric(T, blockSize, evals, V);
+
+  // Step 4: sort eigenvalues descending, reordering eigenvector columns
+  for (uint8_t i = 0; i < blockSize - 1; i++) {
+    for (uint8_t j = i + 1; j < blockSize; j++) {
+      if (evals[j] > evals[i]) {
+        float tmpE = evals[i];
+        evals[i] = evals[j];
+        evals[j] = tmpE;
+        for (uint8_t k = 0; k < blockSize; k++) {
+          float tmpV = V[k][i];
+          V[k][i] = V[k][j];
+          V[k][j] = tmpV;
+        }
+      }
+    }
+  }
+
+  // Step 5: Ublock = B_orig · V · Σ⁻¹, from the SNAPSHOT values.
+  // Column i of Ublock is u_i = (B_orig · v_i) / σ_i.
+  float Ublock[5][5] = {{0}};
+  for (uint8_t i = 0; i < blockSize; i++) {
+    float sigmaI = sqrtf(evals[i]);
+    for (uint8_t r = 0; r < blockSize; r++) {
+      float result = d[r] * V[r][i];
+      if (r + 1 < blockSize) {
+        result += e[r] * V[r + 1][i];
+      }
+      Ublock[r][i] = (sigmaI > 1e-30f) ? result / sigmaI : 0.0f;
+    }
+  }
+
+  // Step 6: fold the factors into the accumulators
+  SVD::ApplyBlockFactorsToAccumulators(blockStart, blockSize, Ublock, V,
+                                       rowsQL, rowsQR, QL, QR);
+
+  // Step 7: W last — write singular values onto the diagonal and zero
+  // the block's superdiagonals
+  for (uint8_t i = 0; i < blockSize; i++) {
+    W[blockStart + i][blockStart + i] = sqrtf(evals[i]);
+    if (i < blockSize - 1) {
+      W[blockStart + i][blockStart + i + 1] = 0.0f;
+    }
+  }
+}
+
+// ============================================================================
+// Phase 3: Extract and Sort Singular Values
+// ============================================================================
+
+void SVD::ExtractAndSortSingularValues(Matrix<5, 5> &W,
+                                        Matrix<5, 1> &sigma,
+                                        uint8_t p,
+                                        Matrix<5, 5> &QL,
+                                        Matrix<5, 5> &QR) {
+  // Extract singular values as absolute values of diagonal elements.
+  // If a diagonal element is negative, flip the sign of the corresponding
+  // column in QL to maintain U * Sigma * Vt = A.
+  for (uint8_t i = 0; i < p; i++) {
+    if (W[i][i] < 0.0f) {
+      // Flip sign of column i in QL
+      for (uint8_t k = 0; k < 5; k++) {
+        QL[k][i] = -QL[k][i];
+      }
+    }
+    sigma[i][0] = fabsf(W[i][i]);
+  }
+
+  // Sort singular values in descending order and reorder U, V accordingly
+  for (uint8_t i = 0; i < p - 1; i++) {
+    for (uint8_t j = i + 1; j < p; j++) {
+      if (sigma[j][0] > sigma[i][0]) {
+        // Swap singular values
+        float tmpS = sigma[i][0];
+        sigma[i][0] = sigma[j][0];
+        sigma[j][0] = tmpS;
+
+        // Swap columns of QL
+        for (uint8_t k = 0; k < 5; k++) {
+          float tmpQ = QL[k][i];
+          QL[k][i] = QL[k][j];
+          QL[k][j] = tmpQ;
+        }
+
+        // Swap columns of QR
+        for (uint8_t k = 0; k < 5; k++) {
+          float tmpQ = QR[k][i];
+          QR[k][i] = QR[k][j];
+          QR[k][j] = tmpQ;
+        }
+      }
+    }
+  }
+}
+
+// ============================================================================
+// Phase 4: Assemble Final U and Vt Matrices
+// ============================================================================
+
+void SVD::AssembleUAndVt(uint8_t m, uint8_t n, uint8_t p,
+                         bool transposeNeeded,
+                         const Matrix<5, 5> &QL,
+                         const Matrix<5, 5> &QR,
+                         Matrix<5, 5> &U,
+                         Matrix<5, 5> &Vt) {
+  // Initialize output matrices to zero
+  for (uint8_t i = 0; i < 5; i++)
+    for (uint8_t j = 0; j < 5; j++) {
+      U[i][j] = 0;
+      Vt[i][j] = 0;
+    }
+
+  // ---- Compute Final U and Vt ----
+  // After bidiagonalization, A = QL · W · QRᵀ (QL = product of left
+  // Householders in application order, QR = product of right Householders),
+  // and the block solvers fold the block SVD factors in: QL <- QL·U_block,
+  // QR <- QR·V_block. Hence:
+  //
+  // Non-transposed (m ≥ n): A = (QL)·Σ·(QR)ᵀ
+  //   U = QL (first m rows, first p columns)
+  //   Vt = QRᵀ (first p rows of the n×n matrix)
+  //
+  // Transposed (wide, m < n): we computed SVD of Aᵀ = QL·W·QRᵀ, so
+  //   A = (QR)·Σ·(QL)ᵀ
+  //   U = QR (first m rows, first p columns)  — NOT QRᵀ
+  //   Vt = QLᵀ (the FULL n×n transpose: QL is the left singular-vector
+  //     matrix of Aᵀ and has n = rows(W) meaningful rows, so Vt needs all
+  //     n rows, not just p)
+
+  for (uint8_t i = 0; i < m; i++) {
+    for (uint8_t j = 0; j < n; j++) {
+      if (j < p) {
+        if (transposeNeeded) {
+          U[i][j] = QR.Get(i, j);
+        } else {
+          U[i][j] = QL.Get(i, j);
+        }
+      } else {
+        U[i][j] = 0;
+      }
+    }
+  }
+
+  for (uint8_t i = 0; i < n; i++) {
+    for (uint8_t j = 0; j < n; j++) {
+      if (transposeNeeded) {
+        Vt[i][j] = QL.Get(j, i);
+      } else if (i < p) {
+        Vt[i][j] = QR.Get(j, i);
+      } else {
+        Vt[i][j] = 0;
+      }
+    }
+  }
+}
+
+// ============================================================================
+// SVD Implementation - Golub-Kahan-Reinsch Algorithm
+// ============================================================================
+/**
+ * @brief SVD for any m×n matrix using Householder bidiagonalization +
+ *        implicit QR iteration on the bidiagonal form.
+ *
+ * Given A (m×n), computes U (m×k), Σ (k×k diagonal), Vᵀ (k×n) where
+ * k = min(m,n) and A = U·Σ·Vᵀ.
+ *
+ * We store results as:
+ *   - U: Matrix<m, n> — first k columns are meaningful
+ *   - sigma: Matrix<n, 1> — first k entries are non-zero singular values
+ *   - Vt: Matrix<n, n> — first k rows are meaningful
+ *
+ * For m < n (wide matrices), we work with Aᵀ and swap roles of U and V.
+ */
+template <uint8_t rows, uint8_t columns>
+void SVD::SVD(Matrix<rows, columns> &matrixToDecompose,
+              Matrix<rows, columns> &U, Matrix<columns, 1> &sigma,
+              Matrix<columns, columns> &Vt) {
+  static_assert(rows <= 5 && columns <= 5,
+                "SVD currently supports matrices up to 5×5");
+
+  uint8_t m = rows;
+  uint8_t n = columns;
+  uint8_t p = (m < n) ? m : n; // rank = min(m,n)
+
+  // For wide matrices (m < n), work with Aᵀ instead.
+  // SVD(A) = U·Σ·Vᵀ ⟺ SVD(Aᵀ) = V·Σ·Uᵀ
+  // So if we compute SVD(Aᵀ) = Ũ·Σ·Ṽᵀ, then U = Ṽ and Vt = Ũᵀ.
+  bool transposeNeeded = (m < n);
+
+  // Working matrix: always p×p or larger square
+  Matrix<5, 5> W{0};
+  for (uint8_t i = 0; i < m; i++) {
+    for (uint8_t j = 0; j < n; j++) {
+      float val = matrixToDecompose.Get(i, j);
+      if (transposeNeeded) {
+        W[j][i] = val; // store Aᵀ
+      } else {
+        W[i][j] = val;
+      }
+    }
+  }
+
+  // After bidiagonalization, W holds the bidiagonal matrix B.
+  // Q_L and Q_R accumulate the Householder transformations.
+  Matrix<5, 5> QL{0}, QR{0};
+  for (uint8_t i = 0; i < 5; i++) {
+    QL[i][i] = 1;
+    QR[i][i] = 1;
+  }
+
+  // ---- Phase 1: Householder Bidiagonalization ----
+  // For non-transpose (m ≥ n): W is m×n, bidiagonalize to get B (m×n)
+  // For transpose (m < n): W is n×m (= Aᵀ), bidiagonalize to get B (n×m)
+  // Pass the ACTUAL dimensions of W: Bidiagonalize needs the full row count
+  // so that the last left Householder (k = p-1, len = rowsW - k) folds the
+  // extra rows into the last diagonal element and zeros them out.
+  {
+    uint8_t rowsW = transposeNeeded ? n : m; // rows of W
+    uint8_t colsW = transposeNeeded ? m : n; // columns of W
+    SVD::Bidiagonalize(W, rowsW, colsW, p, QL, QR);
+  }
+
+  // ---- Phase 2: QR Iteration on Bidiagonal Matrix -->
+  // W now contains the upper bidiagonal matrix B.
+  // We apply QR iterations to converge superdiagonal elements to zero,
+  // leaving singular values on the diagonal.
+  //
+  // rowsQL / rowsQR are the meaningful row extents of the accumulators.
+  // QL is the LEFT factor of the bidiagonalized working matrix W, so it
+  // carries rowsW = (transposeNeeded ? n : m) meaningful rows: in the wide
+  // (transposed) case Vt = QLᵀ needs ALL n rows, so block factors must be
+  // applied over 0..n−1. QR is only ever read back over its first m rows
+  // (as U), but applying factors over all n rows is harmless and matches
+  // the full-row Householder application in Bidiagonalize.
+  uint8_t rowsQL = transposeNeeded ? n : m;
+  uint8_t rowsQR = n;
+
+  uint32_t maxIter = 1000;
+  float tol = 1e-8f;
+
+  for (uint32_t iter = 0; iter < maxIter; iter++) {
+    // Deflate: zero out negligible SUPERDIAGONAL elements
+    SVD::DeflateBidiagonal(W, p, tol);
+
+    // If all superdiagonal elements are zero, we're done
+    if (SVD::BidiagonalIsDiagonal(W, p, tol))
+      break;
+
+    // Process all unreduced blocks in the matrix
+    bool processedAny = false;
+    uint8_t blockStart = 0;
+
+    while (blockStart < p - 1) {
+      // Find end of current unreduced block
+      uint8_t blockEnd = blockStart;
+      while (blockEnd < p - 1 &&
+             fabsf(W[blockEnd][blockEnd + 1]) >
+                 tol * fmaxf(fabsf(W[blockEnd][blockEnd]) +
+                                 fabsf(W[blockEnd + 1][blockEnd + 1]),
+                             1e-10f)) {
+        blockEnd++;
+      }
+
+      // blockStart..blockEnd is an unreduced block of size
+      // (blockEnd - blockStart + 1)
+      uint8_t blockSize = blockEnd - blockStart + 1;
+
+      if (blockSize == 2) {
+        // Handle 2×2 block directly using closed-form solution
+        float Ub2[2][2] = {{0}}, Vb2[2][2] = {{0}};
+        float sig[2] = {0, 0};
+        SVD::SolveBidiagonalBlock2x2(W[blockStart][blockStart],
+                                     W[blockStart][blockEnd],
+                                     W[blockEnd][blockEnd], Ub2, Vb2, sig);
+
+        float Ublock[5][5] = {{0}}, Vblock[5][5] = {{0}};
+        for (uint8_t i = 0; i < 2; i++)
+          for (uint8_t j = 0; j < 2; j++) {
+            Ublock[i][j] = Ub2[i][j];
+            Vblock[i][j] = Vb2[i][j];
+          }
+
+        // Apply block factors to the accumulators
+        SVD::ApplyBlockFactorsToAccumulators(blockStart, 2, Ublock, Vblock,
+                                             rowsQL, rowsQR, QL, QR);
+
+        // Store singular values on diagonal, zero the superdiagonal
+        W[blockStart][blockStart] = sig[0];
+        W[blockEnd][blockEnd] = sig[1];
+        W[blockStart][blockEnd] = 0;
+      } else if (blockSize > 2) {
+        // Larger blocks: SVD via eigendecomposition of BᵀB (Jacobi)
+        SVD::SolveBidiagonalBlockJacobi(W, blockStart, blockSize, rowsQL,
+                                        rowsQR, QL, QR, tol);
+      }
+
+      processedAny = true;
+      blockStart = blockEnd + 1; // Move to next block
+    }
+
+    if (!processedAny) {
+      // No unreduced blocks found, but superdiagonal is not all zero
+      // This can happen with numerical issues, just break
+      break;
+    }
+  }
+
+  // ---- Phase 3: Extract and Sort Singular Values ----
+  // Use internal 5×1 buffer for sigma
+  Matrix<5, 1> sigmaInternal{0};
+  ExtractAndSortSingularValues(W, sigmaInternal, p, QL, QR);
+
+  // ---- Phase 4: Assemble Final U and Vt ----
+  // Use internal 5×5 buffers for U and Vt
+  Matrix<5, 5> UInternal{0}, VtInternal{0};
+  AssembleUAndVt(m, n, p, transposeNeeded, QL, QR, UInternal, VtInternal);
+
+  // Copy results to output parameters
+  for (uint8_t i = 0; i < columns; i++) {
+    sigma[i][0] = sigmaInternal.Get(i, 0);
+  }
+  for (uint8_t i = 0; i < rows; i++) {
+    for (uint8_t j = 0; j < columns; j++) {
+      U[i][j] = UInternal.Get(i, j);
+    }
+  }
+  // Vt is a columns×columns (n×n) matrix: BOTH bounds must run over columns.
+  for (uint8_t i = 0; i < columns; i++) {
+    for (uint8_t j = 0; j < columns; j++) {
+      Vt[i][j] = VtInternal.Get(i, j);
+    }
+  }
+}
+
+#endif

src/SVD.hpp

@@ -0,0 +1,345 @@
+#pragma once
+#include "Matrix.hpp"
+
+/**
+ * @brief library that uses Matrix.hpp and performs SVD on a matrix
+ */
+namespace SVD {
+/**
+ * @brief Compute the Singular Value Decomposition (SVD) of this matrix.
+ *
+ * Decomposes A into U × Σ × Vᵀ where:
+ *   - U is an m×k orthogonal matrix (left singular vectors)
+ *   - Σ is a k×k diagonal matrix with non-negative singular values
+ *     (stored as a k×1 column vector)
+ *   - Vᵀ is a k×n orthogonal matrix (right singular vectors, transposed)
+ *   - k = min(m, n)
+ *
+ * The decomposition satisfies: A ≈ U × diag(Σ) × Vᵀ
+ * Singular values are returned in descending order.
+ *
+ * @param U Output: left singular vectors (m×k orthogonal matrix)
+ * @param sigma Output: singular values as k×1 column vector, sorted descending
+ * @param Vt Output: right singular vectors transposed (k×n matrix)
+ *
+ * @note This implementation uses the Golub-Kahan-Reinsch algorithm:
+ *       1. Householder bidiagonalization of A
+ *       2. Implicit QR iteration on the bidiagonal matrix
+ *       3. Accumulation of U and V factors throughout
+ */
+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);
+
+// ========================================================================
+// SVD Building Block Functions (for unit testing)
+// These operate on internal 5×5 working arrays for maximum flexibility.
+// ========================================================================
+
+/**
+ * @brief Compute a Householder reflector vector.
+ *
+ * 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 len Number of valid elements in x
+ * @param v Output: normalized Householder vector (v[0] is the first element)
+ * @param alpha Output: the resulting first element after reflection
+ * @return The norm of the input vector x
+ */
+static float ComputeHouseholder(const float *x, uint8_t len, float *v,
+                                float &alpha);
+
+/**
+ * @brief Apply a Householder reflection from the left.
+ *
+ * Transforms W = (I - 2·v·vᵀ) · W where v operates on rows [startRow..endRow].
+ *
+ * @param W Input/output: matrix to transform (5×5 working array)
+ * @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,
+                                 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].
+ *
+ * @param W Input/output: matrix to transform (5×5 working array)
+ * @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,
+                                  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) 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)
+ *
+ * @param W Input/output: matrix to bidiagonalize (5×5, must be at least p×q)
+ * @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,
+ *         output: QLᵀ such that QLᵀ·W = B)
+ * @param QR Input/output: right Householder accumulation (initialized to identity,
+ *         output: QR such that W·QR = B after left apply)
+ */
+static void Bidiagonalize(Matrix<5, 5> &W,
+                          uint8_t m, uint8_t q, uint8_t p,
+                          Matrix<5, 5> &QL,
+                          Matrix<5, 5> &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.
+ *
+ * @param W Input/output: bidiagonal matrix (5×5 working array, first p×p used)
+ * @param p Size of the bidiagonal matrix (min(rows, columns))
+ * @param tol Relative deflation tolerance (e.g. 1e-8f)
+ */
+static void DeflateBidiagonal(Matrix<5, 5> &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.
+ *
+ * @param W Input: bidiagonal matrix (5×5 working array, 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
+ */
+static bool BidiagonalIsDiagonal(const Matrix<5, 5> &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
+ *   - columns of V are the corresponding eigenvectors (T·V = V·Λ)
+ *
+ * @param T Input/output: symmetric matrix (5×5 storage, first n×n used,
+ *        destroyed in place)
+ * @param n Matrix size (≤ 5)
+ * @param evals Output: eigenvalues, unsorted (evals[i] = T[i][i])
+ * @param V Output: eigenvector matrix, columns are eigenvectors
+ */
+static void JacobiEigenSymmetric(float T[5][5], uint8_t n, float evals[5],
+                                 float V[5][5]);
+
+/**
+ * @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).
+ *
+ * @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 (blockSize×blockSize in 5×5 storage)
+ * @param Vblock Right singular-vector factor of the block (blockSize×blockSize in 5×5 storage)
+ * @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
+ */
+static void ApplyBlockFactorsToAccumulators(uint8_t blockStart,
+                                            uint8_t blockSize,
+                                            const float Ublock[5][5],
+                                            const float Vblock[5][5],
+                                            uint8_t rowsQL, uint8_t rowsQR,
+                                            Matrix<5, 5> &QL,
+                                            Matrix<5, 5> &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. JacobiEigenSymmetric on T → eigenvalues (unsorted) + V
+ *   3. Sort eigenvalues descending, reordering V
+ *   4. Ublock = B_orig · V · Σ⁻¹ (from the snapshot, so W is not
+ *      overwritten before Ublock is computed)
+ *   5. Fold Ublock/Vblock into QL/QR via ApplyBlockFactorsToAccumulators
+ *   6. Only then write sqrt(eigenvalues) into W's diagonal and zero the
+ *      block's superdiagonals
+ *
+ * @param W Input/output: bidiagonal matrix (5×5 working array); 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, ≤ 5)
+ * @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
+ * @param tol (unused: Jacobi convergence tolerance is internal)
+ */
+static void SolveBidiagonalBlockJacobi(Matrix<5, 5> &W, uint8_t blockStart,
+                                       uint8_t blockSize, uint8_t rowsQL,
+                                       uint8_t rowsQR, Matrix<5, 5> &QL,
+                                       Matrix<5, 5> &QR, float tol);
+
+/**
+ * @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.
+ *
+ * @param W Input: bidiagonal matrix (5×5 working array)
+ * @param sigma Output: sorted singular values (5×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)
+ */
+static void ExtractAndSortSingularValues(Matrix<5, 5> &W,
+                                          Matrix<5, 1> &sigma,
+                                          uint8_t p,
+                                          Matrix<5, 5> &QL,
+                                          Matrix<5, 5> &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 = QL[:,0:p]ᵀ
+ *
+ * @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 (5×5)
+ * @param QR Right Householder accumulation (5×5)
+ * @param U Output: left singular vectors (m×n matrix, only first p columns used)
+ * @param Vt Output: right singular vectors transposed (n×n matrix, only first p rows used)
+ */
+static void AssembleUAndVt(uint8_t m, uint8_t n, uint8_t p,
+                           bool transposeNeeded,
+                           const Matrix<5, 5> &QL,
+                           const Matrix<5, 5> &QR,
+                           Matrix<5, 5> &U,
+                           Matrix<5, 5> &Vt);
+
+/**
+ * @brief Compute a Givens rotation that zeros out y.
+ *
+ * Computes c, s such that:
+ *   [c  s] [x] = [r]
+ *   [-s c] [y]   [0]
+ * where r = sqrt(x² + y²).
+ *
+ * @param x First element
+ * @param y Second element (to be zeroed)
+ * @param c Output: cosine of rotation angle
+ * @param s Output: sine of rotation angle
+ */
+static void ComputeGivens(float x, float y, float &c, float &s);
+
+/**
+ * @brief Apply a Givens rotation from the left to rows i and j.
+ *
+ * Applies [c s; -s c] to rows i, j of W (columns startCol..endCol).
+ *
+ * @param W Input/output: matrix to transform
+ * @param i First row index
+ * @param j Second row index
+ * @param c Cosine of rotation angle
+ * @param s Sine of rotation angle
+ * @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,
+                            float s, uint8_t startCol, uint8_t endCol);
+
+/**
+ * @brief Apply a Givens rotation from the right to columns i and j.
+ *
+ * Applies [c -s; s c]ᵀ to columns i, j of W (rows startRow..endRow).
+ *
+ * @param W Input/output: matrix to transform
+ * @param i First column index
+ * @param j Second column index
+ * @param c Cosine of rotation angle
+ * @param s Sine of rotation angle
+ * @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,
+                             float s, uint8_t startRow, uint8_t endRow);
+} // namespace SVD
+
+#ifndef SVD_H_
+#include "SVD.cpp"
+#endif // SVD_H_

src/Vector3D.cpp

@@ -0,0 +1,162 @@
+#ifdef VECTOR3D_H_ // since the .cpp file has to be included by the .hpp file this
+                   // will evaluate to true
+#include <cmath>
+#include <type_traits>
+#include <string>
+
+template <typename Type>
+V3D<Type>::V3D(const Matrix<1, 3> &other)
+{
+    this->x = other.Get(0, 0);
+    this->y = other.Get(0, 1);
+    this->z = other.Get(0, 2);
+}
+
+template <typename Type>
+V3D<Type>::V3D(const Matrix<3, 1> &other)
+{
+    this->x = other.Get(0, 0);
+    this->y = other.Get(1, 0);
+    this->z = other.Get(2, 0);
+}
+
+template <typename Type>
+V3D<Type>::V3D(const V3D &other) : x(other.x),
+                                   y(other.y),
+                                   z(other.z)
+{
+    static_assert(std::is_arithmetic<Type>::value, "Type must be a number");
+}
+
+template <typename Type>
+V3D<Type>::V3D(Type x, Type y, Type z) : x(x),
+                                         y(y),
+                                         z(z)
+{
+    static_assert(std::is_arithmetic<Type>::value, "Type must be a number");
+}
+
+template <typename Type>
+template <typename OtherType>
+V3D<Type>::V3D(const V3D<OtherType> &other)
+{
+    static_assert(std::is_arithmetic<Type>::value, "Type must be a number");
+    static_assert(std::is_arithmetic<OtherType>::value, "OtherType must be a number");
+    this->x = static_cast<Type>(other.x);
+    this->y = static_cast<Type>(other.y);
+    this->z = static_cast<Type>(other.z);
+}
+
+template <typename Type>
+std::array<Type, 3> V3D<Type>::ToArray() const
+{
+    return {this->x, this->y, this->z};
+}
+
+template <typename Type>
+void V3D<Type>::operator=(const V3D<Type> &other)
+{
+    this->x = other.x;
+    this->y = other.y;
+    this->z = other.z;
+}
+
+template <typename Type>
+V3D<Type> V3D<Type>::operator+(Type other) const
+{
+    return V3D<Type>{this->x + other, this->y + other, this->z + other};
+}
+
+template <typename Type>
+V3D<Type> V3D<Type>::operator+(const V3D<Type> &other) const
+{
+    return V3D<Type>{this->x + other.x, this->y + other.y, this->z + other.z};
+}
+
+template <typename Type>
+V3D<Type> V3D<Type>::operator-(Type other) const
+{
+    return V3D<Type>{this->x - other, this->y - other, this->z - other};
+}
+
+template <typename Type>
+V3D<Type> V3D<Type>::operator-(const V3D<Type> &other) const
+{
+    return V3D<Type>{this->x - other.x, this->y - other.y, this->z - other.z};
+}
+
+template <typename Type>
+V3D<Type> V3D<Type>::operator*(Type scalar) const
+{
+    return V3D<Type>{this->x * scalar, this->y * scalar, this->z * scalar};
+}
+
+template <typename Type>
+V3D<Type> V3D<Type>::operator/(Type scalar) const
+{
+    return V3D<Type>{this->x / scalar, this->y / scalar, this->z / scalar};
+}
+
+template <typename Type>
+V3D<Type> &V3D<Type>::operator+=(Type other)
+{
+    *this = *this + other;
+    return *this;
+}
+
+template <typename Type>
+V3D<Type> &V3D<Type>::operator+=(const V3D<Type> &other)
+{
+    *this = *this + other;
+    return *this;
+}
+
+template <typename Type>
+V3D<Type> &V3D<Type>::operator-=(Type other)
+{
+    *this = *this - other;
+    return *this;
+}
+
+template <typename Type>
+V3D<Type> &V3D<Type>::operator-=(const V3D<Type> &other)
+{
+    *this = *this - other;
+    return *this;
+}
+
+template <typename Type>
+V3D<Type> &V3D<Type>::operator/=(Type scalar)
+{
+    if (scalar == 0)
+    {
+        return *this;
+    }
+    this->x /= scalar;
+    this->y /= scalar;
+    this->z /= scalar;
+    return *this;
+}
+
+template <typename Type>
+V3D<Type> &V3D<Type>::operator*=(Type scalar)
+{
+    this->x *= scalar;
+    this->y *= scalar;
+    this->z *= scalar;
+    return *this;
+}
+
+template <typename Type>
+bool V3D<Type>::operator==(const V3D<Type> &other)
+{
+    return this->x == other.x && this->y == other.y && this->z == other.z;
+}
+
+template <typename Type>
+float V3D<Type>::magnitude()
+{
+    return std::sqrt(static_cast<float>(this->x * this->x + this->y * this->y + this->z * this->z));
+}
+
+#endif // VECTOR3D_H_

src/Vector3D.hpp

@@ -0,0 +1,58 @@
+#ifndef VECTOR3D_H_
+#define VECTOR3D_H_
+
+#include <cstdint>
+#include "Matrix.hpp"
+
+template <typename Type>
+class V3D
+{
+public:
+    V3D(const Matrix<1, 3> &other);
+    V3D(const Matrix<3, 1> &other);
+
+    V3D(const V3D &other);
+
+    V3D(Type x = 0, Type y = 0, Type z = 0);
+
+    template <typename OtherType>
+    V3D(const V3D<OtherType> &other);
+
+    template <typename OtherType>
+    operator OtherType() const;
+
+    std::array<Type, 3> ToArray() const;
+
+    V3D<Type> operator+(Type other) const;
+    V3D<Type> operator+(const V3D<Type> &other) const;
+
+    V3D<Type> operator-(Type other) const;
+    V3D<Type> operator-(const V3D<Type> &other) const;
+
+    V3D<Type> operator*(Type scalar) const;
+
+    V3D<Type> operator/(Type scalar) const;
+
+    void operator=(const V3D<Type> &other);
+
+    V3D<Type> &operator+=(Type other);
+    V3D<Type> &operator+=(const V3D<Type> &other);
+
+    V3D<Type> &operator-=(Type other);
+    V3D<Type> &operator-=(const V3D<Type> &other);
+
+    V3D<Type> &operator/=(Type scalar);
+
+    V3D<Type> &operator*=(Type scalar);
+
+    bool operator==(const V3D<Type> &other);
+
+    float magnitude();
+
+    Type x;
+    Type y;
+    Type z;
+};
+
+#include "Vector3D.cpp"
+#endif // VECTOR3D_H_

unit-tests/CMakeLists.txt

@@ -1,21 +1,55 @@
-cmake_minimum_required (VERSION 3.11)
+# Quaternion tests
+add_executable(quaternion-tests quaternion-tests.cpp)
  
-project ("test_driver")
+target_link_libraries(quaternion-tests 
- 
+    PRIVATE 
-include(FetchContent)
+    quaternion
- 
+    Catch2::Catch2WithMain
-FetchContent_Declare(
-    Catch2
-    GIT_REPOSITORY https://github.com/catchorg/Catch2.git
-    GIT_TAG v3.0.1 # or a later release
 )
 
-FetchContent_MakeAvailable(Catch2)
+# matrix tests
- 
 add_executable(matrix-tests matrix-tests.cpp)
  
 target_link_libraries(matrix-tests 
     PRIVATE 
-    Matrix
+    matrix
+    Catch2::Catch2WithMain
+)
+
+# matrix timing tests
+add_executable(matrix-timing-tests matrix-timing-tests.cpp)
+ 
+target_link_libraries(matrix-timing-tests 
+    PRIVATE 
+    matrix
+    Catch2::Catch2WithMain
+)
+
+# Vector 3D Tests
+add_executable(vector-3d-tests vector-tests.cpp)
+
+target_link_libraries(vector-3d-tests
+    PRIVATE
+    vector-3d
+    Catch2::Catch2WithMain
+)
+
+# SVD building block tests
+add_executable(svd-build-blocks-tests svd-build-blocks-tests.cpp)
+
+target_link_libraries(svd-build-blocks-tests
+    PRIVATE
+    matrix
+    svd
+    Catch2::Catch2WithMain
+)
+
+# SVD integration tests
+add_executable(svd-integration-test svd-integration-test.cpp)
+
+target_link_libraries(svd-integration-test
+    PRIVATE
+    matrix
+    svd
     Catch2::Catch2WithMain
 )

unit-tests/matrix-test-timings.txt

@@ -1,14 +0,0 @@
-Addition: 0.419 s
-Subtraction: 0.421 s
-Multiplication: 3.297 s
-Scalar Multiplication: 0.329 s
-Element Multiply: 0.306 s
-Element Divide: 0.302 s
-Minor Matrix: 0.331 s
-Determinant: 0.177 s
-Matrix of Minors: 0.766 s
-Invert: 0.183 s
-Transpose: 0.215 s
-Normalize: 0.315 s
-GET ROW: 0.008 s
-GET COLUMN: 0.43 s

unit-tests/matrix-tests.cpp

@@ -4,47 +4,68 @@
 
 // include the module you're going to test next
 #include "Matrix.hpp"
+#include "SVD.hpp"
 
 // any other libraries
 #include <array>
 #include <cmath>
 #include <iostream>
 
+// Helper functions
+template <uint8_t rows, uint8_t columns>
+float matrixSum(const Matrix<rows, columns> &matrix) {
+  float sum = 0;
+  for (uint32_t i = 0; i < rows * columns; i++) {
+    float number = matrix.ToArray()[i];
+    sum += number * number;
+  }
+  return std::sqrt(sum);
+}
+
+template <uint8_t rows, uint8_t columns>
+void printLabeledMatrix(const std::string &label,
+                        const Matrix<rows, columns> &matrix) {
+  std::string strBuf = "";
+  matrix.ToString(strBuf);
+  std::cout << label << ":\n" << strBuf << std::endl;
+}
+
+TEST_CASE("Initialization", "Matrix") {
+  SECTION("Array Initialization") {
+    std::array<float, 4> arr2{5, 6, 7, 8};
+    Matrix<2, 2> mat2{arr2};
+    // array initialization
+    REQUIRE(mat2.Get(0, 0) == 5);
+    REQUIRE(mat2.Get(0, 1) == 6);
+    REQUIRE(mat2.Get(1, 0) == 7);
+    REQUIRE(mat2.Get(1, 1) == 8);
+  }
+
+  SECTION("Argument Pack Initialization") {
+    Matrix<2, 2> mat1{1, 2, 3, 4};
+    // template pack initialization
+    REQUIRE(mat1.Get(0, 0) == 1);
+    REQUIRE(mat1.Get(0, 1) == 2);
+    REQUIRE(mat1.Get(1, 0) == 3);
+    REQUIRE(mat1.Get(1, 1) == 4);
+  }
+
+  SECTION("Single Argument Pack Initialization") {
+    Matrix<2, 2> mat1{2};
+    // template pack initialization
+    REQUIRE(mat1.Get(0, 0) == 2);
+    REQUIRE(mat1.Get(0, 1) == 2);
+    REQUIRE(mat1.Get(1, 0) == 2);
+    REQUIRE(mat1.Get(1, 1) == 2);
+  }
+}
+
 TEST_CASE("Elementary Matrix Operations", "Matrix") {
   std::array<float, 4> arr2{5, 6, 7, 8};
   Matrix<2, 2> mat1{1, 2, 3, 4};
   Matrix<2, 2> mat2{arr2};
   Matrix<2, 2> mat3{};
 
-  SECTION("Initialization") {
-    // array initialization
-    REQUIRE(mat1.Get(0, 0) == 1);
-    REQUIRE(mat1.Get(0, 1) == 2);
-    REQUIRE(mat1.Get(1, 0) == 3);
-    REQUIRE(mat1.Get(1, 1) == 4);
-
-    // empty initialization
-    REQUIRE(mat3.Get(0, 0) == 0);
-    REQUIRE(mat3.Get(0, 1) == 0);
-    REQUIRE(mat3.Get(1, 0) == 0);
-    REQUIRE(mat3.Get(1, 1) == 0);
-
-    // template pack initialization
-    REQUIRE(mat2.Get(0, 0) == 5);
-    REQUIRE(mat2.Get(0, 1) == 6);
-    REQUIRE(mat2.Get(1, 0) == 7);
-    REQUIRE(mat2.Get(1, 1) == 8);
-
-    // large matrix
-    Matrix<255, 255> mat6{};
-    mat6.Fill(4);
-    for (uint8_t row{0}; row < 255; row++) {
-      for (uint8_t column{0}; column < 255; column++) {
-        REQUIRE(mat6.Get(row, column) == 4);
-      }
-    }
-  }
-
   SECTION("Fill") {
     mat1.Fill(0);
     REQUIRE(mat1.Get(0, 0) == 0);
@@ -66,10 +87,6 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
   }
 
   SECTION("Addition") {
-    std::string strBuf1 = "";
-    mat1.ToString(strBuf1);
-    std::cout << "Matrix 1:\n" << strBuf1 << std::endl;
-
     mat1.Add(mat2, mat3);
 
     REQUIRE(mat3.Get(0, 0) == 6);
@@ -118,6 +135,36 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
     REQUIRE(mat3.Get(0, 1) == 22);
     REQUIRE(mat3.Get(1, 0) == 43);
     REQUIRE(mat3.Get(1, 1) == 50);
+
+    // Non-square multiplication
+    Matrix<2, 4> mat4{1, 2, 3, 4, 5, 6, 7, 8};
+    Matrix<4, 2> mat5{9, 10, 11, 12, 13, 14, 15, 16};
+    Matrix<2, 2> mat6{};
+    mat6 = mat4 * mat5;
+    REQUIRE(mat6.Get(0, 0) == 130);
+    REQUIRE(mat6.Get(0, 1) == 140);
+    REQUIRE(mat6.Get(1, 0) == 322);
+    REQUIRE(mat6.Get(1, 1) == 348);
+
+    // One more non-square multiplicaiton
+    Matrix<4, 4> mat7{};
+    mat7 = mat5 * mat4;
+    REQUIRE(mat7.Get(0, 0) == 59);
+    REQUIRE(mat7.Get(0, 1) == 78);
+    REQUIRE(mat7.Get(0, 2) == 97);
+    REQUIRE(mat7.Get(0, 3) == 116);
+    REQUIRE(mat7.Get(1, 0) == 71);
+    REQUIRE(mat7.Get(1, 1) == 94);
+    REQUIRE(mat7.Get(1, 2) == 117);
+    REQUIRE(mat7.Get(1, 3) == 140);
+    REQUIRE(mat7.Get(2, 0) == 83);
+    REQUIRE(mat7.Get(2, 1) == 110);
+    REQUIRE(mat7.Get(2, 2) == 137);
+    REQUIRE(mat7.Get(2, 3) == 164);
+    REQUIRE(mat7.Get(3, 0) == 95);
+    REQUIRE(mat7.Get(3, 1) == 126);
+    REQUIRE(mat7.Get(3, 2) == 157);
+    REQUIRE(mat7.Get(3, 3) == 188);
   }
 
   SECTION("Scalar Multiplication") {
@@ -182,18 +229,6 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
     REQUIRE(minorMat4.Get(1, 0) == 4);
     REQUIRE(minorMat4.Get(1, 1) == 5);
   }
-  SECTION("Identity Matrix") {
-    Matrix<3, 3> mat4{Matrix<3, 3>::Eye()};
-    REQUIRE(mat4.Get(0, 0) == 1);
-    REQUIRE(mat4.Get(0, 1) == 0);
-    REQUIRE(mat4.Get(0, 2) == 0);
-    REQUIRE(mat4.Get(1, 0) == 0);
-    REQUIRE(mat4.Get(1, 1) == 1);
-    REQUIRE(mat4.Get(1, 2) == 0);
-    REQUIRE(mat4.Get(2, 0) == 0);
-    REQUIRE(mat4.Get(2, 1) == 0);
-    REQUIRE(mat4.Get(2, 2) == 1);
-  }
 
   SECTION("Determinant") {
     float det1 = mat1.Det();
@@ -234,7 +269,7 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
   }
 
   SECTION("Invert") {
-    mat1.Invert(mat3);
+    mat3 = mat1.Invert();
     REQUIRE_THAT(mat3.Get(0, 0), Catch::Matchers::WithinRel(-2.0F, 1e-6f));
     REQUIRE_THAT(mat3.Get(0, 1), Catch::Matchers::WithinRel(1.0F, 1e-6f));
     REQUIRE_THAT(mat3.Get(1, 0), Catch::Matchers::WithinRel(1.5F, 1e-6f));
@@ -243,7 +278,7 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
 
   SECTION("Transpose") {
     // transpose a square matrix
-    mat1.Transpose(mat3);
+    mat3 = mat1.Transpose();
 
     REQUIRE(mat3.Get(0, 0) == 1);
     REQUIRE(mat3.Get(0, 1) == 3);
@@ -254,7 +289,7 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
     Matrix<2, 3> mat4{1, 2, 3, 4, 5, 6};
     Matrix<3, 2> mat5{};
 
-    mat4.Transpose(mat5);
+    mat5 = mat4.Transpose();
 
     REQUIRE(mat5.Get(0, 0) == 1);
     REQUIRE(mat5.Get(0, 1) == 4);
@@ -264,26 +299,6 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
     REQUIRE(mat5.Get(2, 1) == 6);
   }
 
-  SECTION("Normalize") {
-    mat1.Normalize(mat3);
-
-    float sqrt_30{sqrt(30)};
-
-    REQUIRE(mat3.Get(0, 0) == 1 / sqrt_30);
-    REQUIRE(mat3.Get(0, 1) == 2 / sqrt_30);
-    REQUIRE(mat3.Get(1, 0) == 3 / sqrt_30);
-    REQUIRE(mat3.Get(1, 1) == 4 / sqrt_30);
-
-    Matrix<2, 1> mat4{-0.878877044, 2.92092276};
-    Matrix<2, 1> mat5{};
-    mat4.Normalize(mat5);
-
-    REQUIRE_THAT(mat5.Get(0, 0),
-                 Catch::Matchers::WithinRel(-0.288129855179f, 1e-6f));
-    REQUIRE_THAT(mat5.Get(1, 0),
-                 Catch::Matchers::WithinRel(0.957591346325f, 1e-6f));
-  }
-
   SECTION("GET ROW") {
     Matrix<1, 2> mat1Rows{};
     mat1.GetRow(0, mat1Rows);
@@ -305,113 +320,717 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
     REQUIRE(mat1Columns.Get(0, 0) == 2);
     REQUIRE(mat1Columns.Get(1, 0) == 4);
   }
+
+  SECTION("Get Sub-Matrices") {
+    Matrix<3, 3> mat4{1, 2, 3, 4, 5, 6, 7, 8, 9};
+
+    Matrix<2, 2> mat5 = mat4.SubMatrix<2, 2, 0, 0>();
+
+    REQUIRE(mat5.Get(0, 0) == 1);
+    REQUIRE(mat5.Get(0, 1) == 2);
+    REQUIRE(mat5.Get(1, 0) == 4);
+    REQUIRE(mat5.Get(1, 1) == 5);
+
+    mat5 = mat4.SubMatrix<2, 2, 1, 1>();
+    REQUIRE(mat5.Get(0, 0) == 5);
+    REQUIRE(mat5.Get(0, 1) == 6);
+    REQUIRE(mat5.Get(1, 0) == 8);
+    REQUIRE(mat5.Get(1, 1) == 9);
+
+    Matrix<3, 1> mat6 = mat4.SubMatrix<3, 1, 0, 0>();
+    REQUIRE(mat6.Get(0, 0) == 1);
+    REQUIRE(mat6.Get(1, 0) == 4);
+    REQUIRE(mat6.Get(2, 0) == 7);
+
+    Matrix<1, 3> mat7 = mat4.SubMatrix<1, 3, 0, 0>();
+    REQUIRE(mat7.Get(0, 0) == 1);
+    REQUIRE(mat7.Get(0, 1) == 2);
+    REQUIRE(mat7.Get(0, 2) == 3);
   }
 
-// basically re-run all of the previous tests with huge matrices and time the
+  SECTION("Set Sub-Matrices") {
-// results.
+    Matrix<3, 3> startMatrix{1, 2, 3, 4, 5, 6, 7, 8, 9};
-TEST_CASE("Timing Tests", "Matrix") {
+    Matrix<3, 3> mat4 = startMatrix;
-  std::array<float, 50 * 50> arr1{};
-  for (uint16_t i{0}; i < 50 * 50; i++) {
-    arr1[i] = i;
-  }
-  std::array<float, 50 * 50> arr2{5, 6, 7, 8};
-  for (uint16_t i{50 * 50}; i < 2 * 50 * 50; i++) {
-    arr2[i] = i;
-  }
-  Matrix<50, 50> mat1{arr1};
-  Matrix<50, 50> mat2{arr2};
-  Matrix<50, 50> mat3{};
 
-  // A smaller matrix to use for really badly optimized operations
+    Matrix<2, 2> mat5{10, 11, 12, 13};
-  Matrix<4, 4> mat4{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16};
+    mat4.SetSubMatrix(0, 0, mat5);
-  Matrix<4, 4> mat5{};
+    REQUIRE(mat4.Get(0, 0) == 10);
+    REQUIRE(mat4.Get(0, 1) == 11);
+    REQUIRE(mat4.Get(1, 0) == 12);
+    REQUIRE(mat4.Get(1, 1) == 13);
 
-  SECTION("Addition") {
+    mat4 = startMatrix;
-    for (uint32_t i{0}; i < 10000; i++) {
+    mat4.SetSubMatrix(1, 1, mat5);
-      mat3 = mat1 + mat2;
+    REQUIRE(mat4.Get(1, 1) == 10);
+    REQUIRE(mat4.Get(1, 2) == 11);
+    REQUIRE(mat4.Get(2, 1) == 12);
+    REQUIRE(mat4.Get(2, 2) == 13);
+
+    Matrix<3, 1> mat6{10, 11, 12};
+    mat4.SetSubMatrix(0, 0, mat6);
+    REQUIRE(mat4.Get(0, 0) == 10);
+    REQUIRE(mat4.Get(1, 0) == 11);
+    REQUIRE(mat4.Get(2, 0) == 12);
+
+    Matrix<1, 3> mat7{10, 11, 12};
+    mat4.SetSubMatrix(0, 0, mat7);
+    REQUIRE(mat4.Get(0, 0) == 10);
+    REQUIRE(mat4.Get(0, 1) == 11);
+    REQUIRE(mat4.Get(0, 2) == 12);
   }
 }
 
-  SECTION("Subtraction") {
+TEST_CASE("Identity Matrix", "Matrix") {
-    for (uint32_t i{0}; i < 10000; i++) {
+  SECTION("Square Matrix") {
-      mat3 = mat1 - mat2;
+    Matrix<5, 5> matrix = Matrix<5, 5>::Identity();
+    uint32_t oneColumnIndex{0};
+    for (uint32_t row = 0; row < 5; row++) {
+      for (uint32_t column = 0; column < 5; column++) {
+        float value = matrix[row][column];
+        if (oneColumnIndex == column) {
+          REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
+        } else {
+          REQUIRE_THAT(value, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+        }
+      }
+      oneColumnIndex++;
     }
   }
 
-  SECTION("Multiplication") {
+  SECTION("Wide Matrix") {
-    for (uint32_t i{0}; i < 1000; i++) {
+    Matrix<2, 5> matrix = Matrix<2, 5>::Identity();
-      mat3 = mat1 * mat2;
+
+    uint32_t oneColumnIndex{0};
+    for (uint32_t row = 0; row < 2; row++) {
+      for (uint32_t column = 0; column < 5; column++) {
+        float value = matrix[row][column];
+        if (oneColumnIndex == column && row < 3) {
+          REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
+        } else {
+          REQUIRE_THAT(value, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+        }
+      }
+      oneColumnIndex++;
     }
   }
 
-  SECTION("Scalar Multiplication") {
+  SECTION("Tall Matrix") {
-    for (uint32_t i{0}; i < 10000; i++) {
+    Matrix<5, 2> matrix = Matrix<5, 2>::Identity();
-      mat3 = mat1 * 3;
+    uint32_t oneColumnIndex{0};
+    for (uint32_t row = 0; row < 5; row++) {
+      for (uint32_t column = 0; column < 2; column++) {
+        float value = matrix[row][column];
+        if (oneColumnIndex == column) {
+          REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
+        } else {
+          REQUIRE_THAT(value, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+        }
+      }
+      oneColumnIndex++;
+    }
   }
 }
 
-  SECTION("Element Multiply") {
+// TODO: Add test for scalar division
-    for (uint32_t i{0}; i < 10000; i++) {
+TEST_CASE("Euclidean Norm", "Matrix") {
-      mat1.ElementMultiply(mat2, mat3);
+
+  SECTION("2x2 Normalize") {
+    Matrix<2, 2> mat1{1, 2, 3, 4};
+    Matrix<2, 2> mat2{};
+
+    mat2 = mat1 / mat1.EuclideanNorm();
+
+    float sqrt_30{static_cast<float>(sqrt(30.0f))};
+
+    REQUIRE(mat2.Get(0, 0) == 1 / sqrt_30);
+    REQUIRE(mat2.Get(0, 1) == 2 / sqrt_30);
+    REQUIRE(mat2.Get(1, 0) == 3 / sqrt_30);
+    REQUIRE(mat2.Get(1, 1) == 4 / sqrt_30);
+
+    REQUIRE_THAT(matrixSum(mat2), Catch::Matchers::WithinRel(1.0f, 1e-6f));
+  }
+
+  SECTION("2x1 (Vector) Normalize") {
+    Matrix<2, 1> mat1{-0.878877044, 2.92092276};
+    Matrix<2, 1> mat2{};
+    mat2 = mat1 / mat1.EuclideanNorm();
+
+    REQUIRE_THAT(mat2.Get(0, 0),
+                 Catch::Matchers::WithinRel(-0.288129855179f, 1e-6f));
+    REQUIRE_THAT(mat2.Get(1, 0),
+                 Catch::Matchers::WithinRel(0.957591346325f, 1e-6f));
+
+    float sum = matrixSum(mat2);
+    REQUIRE_THAT(sum, Catch::Matchers::WithinRel(1.0f, 1e-6f));
+  }
+
+  SECTION("Normalized vectors sum to 1") {
+    Matrix<9, 1> mat1{1, 2, 3, 4, 5, 6, 7, 8, 9};
+    Matrix<9, 1> mat2;
+    mat2 = mat1 / mat1.EuclideanNorm();
+    float sum = matrixSum(mat2);
+    REQUIRE_THAT(sum, Catch::Matchers::WithinRel(1.0f, 1e-6f));
+
+    Matrix<2, 3> mat3{1, 2, 3, 4, 5, 6};
+    Matrix<2, 3> mat4{};
+    mat4 = mat3 / mat3.EuclideanNorm();
+    sum = matrixSum(mat4);
+    REQUIRE_THAT(sum, Catch::Matchers::WithinRel(1.0f, 1e-6f));
   }
 }
 
-  SECTION("Element Divide") {
+TEST_CASE("QR Decompositions", "Matrix") {
-    for (uint32_t i{0}; i < 10000; i++) {
+  SECTION("2x2 QRDecomposition") {
-      mat1.ElementDivide(mat2, mat3);
+    Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
+    Matrix<2, 2> Q{}, R{};
+    A.QRDecomposition(Q, R);
+
+    // Check that Q * R ≈ A
+    Matrix<2, 2> QR{};
+    Q.Mult(R, QR);
+    for (int i = 0; i < 2; ++i) {
+      for (int j = 0; j < 2; ++j) {
+        REQUIRE_THAT(QR[i][j], Catch::Matchers::WithinRel(A[i][j], 1e-4f));
       }
     }
 
-  SECTION("Minor Matrix") {
+    // Check that Q is orthonormal: Qᵀ * Q ≈ I
-    // what about matrices of 0,0 or 1,1?
+    Matrix<2, 2> Qt = Q.Transpose();
-    // minor matrix for 2x2 matrix
+    Matrix<2, 2> QtQ{};
-    Matrix<49, 49> minorMat1{};
+    Qt.Mult(Q, QtQ);
-    for (uint32_t i{0}; i < 10000; i++) {
+    for (int i = 0; i < 2; ++i) {
-      mat1.MinorMatrix(minorMat1, 0, 0);
+      for (int j = 0; j < 2; ++j) {
+        if (i == j)
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f));
+        else
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
       }
     }
 
-  SECTION("Determinant") {
+    // Optional: R should be upper triangular
-    for (uint32_t i{0}; i < 100000; i++) {
+    REQUIRE(std::fabs(R[1][0]) < 1e-4f);
-      float det1 = mat4.Det();
+
+    // check that all Q values are correct
+    REQUIRE_THAT(Q[0][0], Catch::Matchers::WithinRel(0.3162f, 1e-4f));
+    REQUIRE_THAT(Q[0][1], Catch::Matchers::WithinRel(0.94868f, 1e-4f));
+    REQUIRE_THAT(Q[1][0], Catch::Matchers::WithinRel(0.94868f, 1e-4f));
+    REQUIRE_THAT(Q[1][1], Catch::Matchers::WithinRel(-0.3162f, 1e-4f));
+
+    // 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::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(R[1][1], Catch::Matchers::WithinRel(0.63246f, 1e-4f));
+  }
+
+  SECTION("3x3 QRDecomposition") {
+    // this symmetrix tridiagonal matrix is well behaved for testing
+    Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
+
+    Matrix<3, 3> Q{}, R{};
+    A.QRDecomposition(Q, R);
+
+    // Check that Q * R ≈ A
+    Matrix<3, 3> QR{};
+    QR = Q * R;
+    for (int i = 0; i < 3; ++i) {
+      for (int j = 0; j < 3; ++j) {
+        REQUIRE_THAT(QR[i][j], Catch::Matchers::WithinRel(A[i][j], 1e-4f));
       }
     }
 
-  SECTION("Matrix of Minors") {
+    // Check that Qᵀ * Q ≈ I
-    for (uint32_t i{0}; i < 100000; i++) {
+    // Since the rank of this matrix is 2, only the top left 2x2 sub-matrix will
-      mat4.MatrixOfMinors(mat5);
+    // equal I.
+    Matrix<3, 3> Qt = Q.Transpose();
+    Matrix<3, 3> QtQ{};
+    QtQ = Qt * Q;
+    for (int i = 0; i < 2; ++i) {
+      for (int j = 0; j < 2; ++j) {
+        if (i == j)
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f));
+        else
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
       }
     }
 
-  SECTION("Invert") {
+    // Optional: Check R is upper triangular
-    for (uint32_t i{0}; i < 100000; i++) {
+    for (int i = 1; i < 3; ++i) {
-      mat4.Invert(mat5);
+      for (int j = 0; j < i; ++j) {
+        REQUIRE(std::fabs(R[i][j]) < 1e-4f);
       }
-  };
-
-  SECTION("Transpose") {
-    for (uint32_t i{0}; i < 10000; i++) {
-      mat1.Transpose(mat3);
     }
   }
 
-  SECTION("Normalize") {
+  SECTION("4x2 QRDecomposition") {
-    for (uint32_t i{0}; i < 10000; i++) {
+    // A simple 4x2 matrix
-      mat1.Normalize(mat3);
+    Matrix<4, 2> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f};
+
+    Matrix<4, 2> Q{};
+    Matrix<2, 2> R{};
+    A.QRDecomposition(Q, R);
+
+    // Check that Q * R ≈ A
+    Matrix<4, 2> QR{};
+    Q.Mult(R, QR);
+    for (int i = 0; i < 4; ++i) {
+      for (int j = 0; j < 2; ++j) {
+        REQUIRE_THAT(QR[i][j], Catch::Matchers::WithinRel(A[i][j], 1e-4f));
       }
     }
 
-  SECTION("GET ROW") {
+    // Check that Qᵀ * Q ≈ I₂
-    Matrix<1, 50> mat1Rows{};
+    Matrix<2, 4> Qt = Q.Transpose();
-    for (uint32_t i{0}; i < 1000000; i++) {
+    Matrix<2, 2> QtQ{};
-      mat1.GetRow(0, mat1Rows);
+    Qt.Mult(Q, QtQ);
+    for (int i = 0; i < 2; ++i) {
+      for (int j = 0; j < 2; ++j) {
+        if (i == j)
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f));
+        else
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
       }
     }
 
-  SECTION("GET COLUMN") {
+    // Check R is upper triangular (i > j ⇒ R[i][j] ≈ 0)
-    Matrix<50, 1> mat1Columns{};
+    for (int i = 1; i < 2; ++i) {
-    for (uint32_t i{0}; i < 1000000; i++) {
+      for (int j = 0; j < i; ++j) {
-      mat1.GetColumn(0, mat1Columns);
+        REQUIRE(std::fabs(R[i][j]) < 1e-4f);
       }
     }
   }
+}
+
+TEST_CASE("Eigenvalues and Vectors", "Matrix") {
+  SECTION("2x2 Eigen") {
+    Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
+    Matrix<2, 2> vectors{};
+    Matrix<2, 1> values{};
+
+    A.EigenQR(vectors, values, 1000000, 1e-20f);
+
+    REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.41597f, 1e-4f));
+    REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.90938f, 1e-4f));
+    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(5.372282f, 1e-4f));
+    REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(-0.372281f, 1e-4f));
+  }
+
+  SECTION("3x3 Rank Defficient Eigen") {
+    SKIP("Skipping this because QR decomposition isn't ready for it");
+    // this symmetrix tridiagonal matrix is well behaved for testing
+    Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
+
+    Matrix<3, 3> vectors{};
+    Matrix<3, 1> values{};
+    A.EigenQR(vectors, values, 1000000, 1e-8f);
+
+    std::string strBuf1 = "";
+    vectors.ToString(strBuf1);
+    std::cout << "Vectors:\n" << strBuf1 << std::endl;
+    strBuf1 = "";
+    values.ToString(strBuf1);
+    std::cout << "Values:\n" << strBuf1 << std::endl;
+
+    REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.23197f, 1e-4f));
+    REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.525322f, 1e-4f));
+    REQUIRE_THAT(vectors[2][0], Catch::Matchers::WithinRel(0.81867f, 1e-4f));
+    REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(-1.11684f, 1e-4f));
+    REQUIRE_THAT(values[1][0], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+    REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(16.1168f, 1e-4f));
+  }
+}
+
+// ============================================================================
+// SVD Tests — Reference values computed via scipy.linalg.svd (Python)
+// ============================================================================
+
+/**
+ * @brief Helper: compute Frobenius norm of a matrix.
+ */
+template <uint8_t rows, uint8_t columns>
+static float frobeniusNorm(const Matrix<rows, columns> &M) {
+  float sum = 0;
+  for (uint8_t i = 0; i < rows; i++) {
+    for (uint8_t j = 0; j < columns; j++) {
+      float v = M.Get(i, j);
+      sum += v * v;
+    }
+  }
+  return sqrtf(sum);
+}
+
+/**
+ * @brief Helper: compute reconstruction error ||A - UΣVᵀ||_F.
+ *
+ * Verifies the fundamental SVD identity A = U × diag(σ) × Vᵀ.
+ * For non-square matrices, only the first min(rows,cols) singular values
+ * contribute to the reconstruction.
+ */
+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) {
+  // Compute U × diag(σ): only first min(rows,cols) columns of U are used
+  constexpr uint8_t k = (rows < columns) ? rows : columns;
+  Matrix<rows, columns> USigma{0};
+  for (uint8_t i = 0; i < rows; i++) {
+    for (uint8_t j = 0; j < k; j++) {
+      USigma[i][j] = U.Get(i, j) * sigma.Get(j, 0);
+    }
+  }
+
+  // Compute (UΣ) × Vᵀ: only first k rows of Vt are used
+  Matrix<rows, columns> UVt{0};
+  for (uint8_t i = 0; i < rows; i++) {
+    for (uint8_t j = 0; j < columns; j++) {
+      float sum = 0;
+      for (uint8_t p = 0; p < k; p++) {
+        sum += USigma[i][p] * Vt.Get(p, j);
+      }
+      UVt[i][j] = sum;
+    }
+  }
+
+  // Compute ||A - UVᵀ||_F
+  Matrix<rows, columns> diff{0};
+  A.Sub(UVt, diff);
+  return frobeniusNorm(diff);
+}
+
+/**
+ * @brief Helper: check orthogonality of the first k columns of M.
+ * Verifies M[:,0:k]ᵀ × M[:,0:k] ≈ I_k.
+ */
+template <uint8_t rows, uint8_t columns>
+static float orthogonalityError(const Matrix<rows, columns> &M) {
+  constexpr uint8_t k = (rows < columns) ? rows : columns;
+
+  // Compute Mᵀ × M (should be I_k in top-left)
+  Matrix<columns, rows> Mt = M.Transpose();
+  Matrix<columns, columns> MtM{0};
+  Mt.Mult(M, MtM);
+
+  float err = 0;
+  for (uint8_t i = 0; i < k; i++) {
+    for (uint8_t j = 0; j < k; j++) {
+      float expected = (i == j) ? 1.0f : 0.0f;
+      err += (MtM.Get(i, j) - expected) * (MtM.Get(i, j) - expected);
+    }
+  }
+  return sqrtf(err);
+}
+
+/**
+ * @brief Helper: check that singular values are sorted in descending order.
+ */
+template <uint8_t maxCols>
+static bool isSortedDescending(const Matrix<maxCols, 1> &sigma, uint8_t count) {
+  for (uint8_t i = 0; i < count - 1; i++) {
+    if (sigma.Get(i + 1, 0) > sigma.Get(i, 0) + 1e-6f) {
+      return false;
+    }
+  }
+  return true;
+}
+
+TEST_CASE("SVD: Simple 2x2 Matrix", "Matrix") {
+  // Reference: scipy.linalg.svd([[1,2],[3,4]])
+  // σ = [5.4649857042, 0.3659661906]
+  Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
+  Matrix<2, 2> U{}, Vt{};
+  Matrix<2, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // Verify singular values (verified with Python scipy.linalg.svd)
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(5.4649857042f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(0.3659661906f, 1e-4f));
+
+  // Verify descending order
+  REQUIRE(isSortedDescending(sigma, 2));
+
+  // Verify U is orthogonal: UᵀU ≈ I
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+
+  // Verify Vt is orthogonal: VtVᵀ ≈ I
+  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::WithinAbs(0.0f, 1e-4f));
+}
+
+TEST_CASE("SVD: Symmetric Positive Definite 2x2", "Matrix") {
+  // Reference: scipy.linalg.svd([[5,3],[3,5]])
+  // σ = [8.0, 2.0] (eigenvalues since symmetric PD)
+  Matrix<2, 2> A{5.0f, 3.0f, 3.0f, 5.0f};
+  Matrix<2, 2> U{}, Vt{};
+  Matrix<2, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(8.0f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.0f, 1e-4f));
+
+  // For symmetric PD matrices, U ≈ V (up to sign)
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+}
+
+TEST_CASE("SVD: Full-Rank 3x3 Matrix", "Matrix") {
+  // Reference: scipy.linalg.svd([[1,2,3],[4,5,6],[7,8,10]])
+  // σ = [17.4125051668, 0.8751613501, 0.1968665211]
+  Matrix<3, 3> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 10.0f};
+  Matrix<3, 3> U{}, Vt{};
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(17.4125051668f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(0.8751613501f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0),
+               Catch::Matchers::WithinRel(0.1968665211f, 1e-4f));
+
+  REQUIRE(isSortedDescending(sigma, 3));
+  REQUIRE_THAT(orthogonalityError(U), Catch::Matchers::WithinAbs(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::WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("SVD: Rank-Deficient 3x3 Matrix", "Matrix") {
+  // Reference: scipy.linalg.svd([[1,2,3],[4,5,6],[7,8,9]])
+  // σ = [16.8481033526, 1.0683695146, ~0] (rank 2)
+  Matrix<3, 3> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f};
+  Matrix<3, 3> U{}, Vt{};
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(16.8481033526f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(1.0683695146f, 1e-4f));
+  // Third singular value should be ~0 (rank deficiency)
+  REQUIRE(sigma.Get(2, 0) < 1e-3f);
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  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};
+  Matrix<3, 3> U{}, Vt{};
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(10.0f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(5.0f, 1e-4f));
+  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::WithinAbs(0.0f, 1e-4f));
+}
+
+TEST_CASE("SVD: Tall Matrix (4×3)", "Matrix") {
+  // Reference: scipy.linalg.svd with full_matrices=False
+  // σ = [25.4624074360, 1.2906616758, ~0] (rank 2)
+  Matrix<4, 3> A{1.0f, 2.0f, 3.0f, 4.0f,  5.0f,  6.0f,
+                 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f};
+  Matrix<4, 3> U{};
+  Matrix<3, 3> Vt{}; // Vt is always n×n
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(25.4624074360f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(1.2906616758f, 1e-4f));
+  REQUIRE(sigma.Get(2, 0) < 1e-3f);
+
+  // U should be 4×3 with orthonormal columns
+  REQUIRE_THAT(orthogonalityError(U), 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 Matrix (3×5)", "Matrix") {
+  // Reference: scipy.linalg.svd with full_matrices=False
+  // σ = [35.1272233336, 2.4653966969, ~0] (rank 2)
+  Matrix<3, 5> A{1.0f, 2.0f,  3.0f,  4.0f,  5.0f,  6.0f,  7.0f, 8.0f,
+                 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f};
+  Matrix<3, 5> U{};
+  Matrix<5, 5> Vt{}; // Vt is always n×n
+  Matrix<5, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(35.1272233336f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(2.4653966969f, 1e-4f));
+  REQUIRE(sigma.Get(2, 0) < 1e-3f);
+
+  // Vt should be 5×5 with orthonormal rows (first k)
+  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 Symmetric Tridiagonal", "Matrix") {
+  // Reference: scipy.linalg.svd for discrete Laplacian-like matrix
+  // σ = [3.7320508076, 3.0, 2.0, 1.0, 0.2679491924]
+  Matrix<5, 5> A{2.0f, -1.0f, 0.0f,  0.0f, 0.0f,  -1.0f, 2.0f, -1.0f, 0.0f,
+                 0.0f, 0.0f,  -1.0f, 2.0f, -1.0f, 0.0f,  0.0f, 0.0f,  -1.0f,
+                 2.0f, -1.0f, 0.0f,  0.0f, 0.0f,  -1.0f, 2.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(3.7320508076f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(3.0f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.0f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(3, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(4, 0),
+               Catch::Matchers::WithinRel(0.2679491924f, 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") {
+  // Reference: scipy.linalg.svd([[0.5,-0.3,0.8],[-0.2,0.7,0.1]])
+  // σ = [1.0384009867, 0.6646227432]
+  Matrix<2, 3> A{0.5f, -0.3f, 0.8f, -0.2f, 0.7f, 0.1f};
+  Matrix<2, 3> U{};
+  Matrix<3, 3> Vt{}; // Vt is always n×n
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0),
+               Catch::Matchers::WithinRel(1.0384009867f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0),
+               Catch::Matchers::WithinRel(0.6646227432f, 1e-4f));
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+}
+
+TEST_CASE("SVD: Near-Singular 2×2 Matrix", "Matrix") {
+  // Condition number ≈ 1e6 — tests numerical stability
+  // Reference: scipy.linalg.svd([[1,0],[0,1e-6]])
+  // σ = [1.0, 1e-6]
+  Matrix<2, 2> A{1.0f, 0.0f, 0.0f, 1e-6f};
+  Matrix<2, 2> U{}, Vt{};
+  Matrix<2, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
+  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::WithinAbs(0.0f, 1e-6f));
+}
+
+TEST_CASE("SVD: Orthogonal Matrix (3×3)", "Matrix") {
+  // For an orthogonal matrix, all singular values should be 1.
+  // Rotation matrix about z-axis by 45°
+  float c = sqrtf(0.5f); // cos(45°)
+  float s = sqrtf(0.5f); // sin(45°)
+  Matrix<3, 3> A{c, -s, 0.0f, s, c, 0.0f, 0.0f, 0.0f, 1.0f};
+  Matrix<3, 3> U{}, Vt{};
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // All singular values should be 1 for an orthogonal matrix
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
+  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::WithinAbs(0.0f, 1e-4f));
+}
+
+TEST_CASE("SVD: Identity Matrix", "Matrix") {
+  // For I, σ = [1, 1, 1], U = V = I
+  Matrix<3, 3> A{1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f};
+  Matrix<3, 3> U{}, Vt{};
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.0f, 1e-4f));
+  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::WithinAbs(0.0f, 1e-6f));
+}
+
+TEST_CASE("SVD: Zero Matrix", "Matrix") {
+  // All singular values should be zero
+  Matrix<3, 3> A{0.0f};
+  Matrix<3, 3> U{}, Vt{};
+  Matrix<3, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE(sigma.Get(0, 0) < 1e-6f);
+  REQUIRE(sigma.Get(1, 0) < 1e-6f);
+  REQUIRE(sigma.Get(2, 0) < 1e-6f);
+
+  float reconErr = svdReconstructionError(A, U, sigma, Vt);
+  REQUIRE_THAT(reconErr, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
+}
+
+TEST_CASE("SVD: 2×1 Column Vector", "Matrix") {
+  // For a column vector v, σ = ||v||, U = v/||v|| (with padding)
+  Matrix<2, 1> A{3.0f, 4.0f};
+  Matrix<2, 1> U{};
+  Matrix<1, 1> Vt{}; // Vt is always n×n
+  Matrix<1, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // σ should be the Euclidean norm: ||[3,4]|| = 5
+  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::WithinAbs(0.0f, 1e-4f));
+}
+
+TEST_CASE("SVD: 1×2 Row Vector", "Matrix") {
+  // For a row vector vᵀ, σ = ||v||, Vt = v/||v|| (with padding)
+  Matrix<1, 2> A{3.0f, 4.0f};
+  Matrix<1, 2> U{};
+  Matrix<2, 2> Vt{}; // Vt is always n×n
+  Matrix<2, 1> sigma{};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // σ should be the Euclidean norm: ||[3,4]|| = 5
+  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::WithinAbs(0.0f, 1e-4f));
+}

unit-tests/matrix-timing-tests.cpp

@@ -0,0 +1,128 @@
+// 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"
+
+// any other libraries
+#include <array>
+#include <cmath>
+#include <cstdint>
+
+// basically re-run all of the matrix tests with huge matrices and time the
+// results.
+TEST_CASE("Timing Tests", "Matrix") {
+  std::array<float, 50 * 50> arr1{};
+  for (uint16_t i{0}; i < 50 * 50; i++) {
+    arr1[i] = i;
+  }
+  std::array<float, 50 * 50> arr2{5, 6, 7, 8};
+  for (uint16_t i{50 * 50}; i < 2 * 50 * 50; i++) {
+    arr2[i] = i;
+  }
+  Matrix<50, 50> mat1{arr1};
+  Matrix<50, 50> mat2{arr2};
+  Matrix<50, 50> mat3{};
+
+  // A smaller matrix to use for really badly optimized operations
+  Matrix<4, 4> mat4{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16};
+  Matrix<4, 4> mat5{};
+
+  SECTION("Addition") {
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat3 = mat1 + mat2;
+    }
+  }
+
+  SECTION("Subtraction") {
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat3 = mat1 - mat2;
+    }
+  }
+
+  SECTION("Multiplication") {
+    for (uint32_t i{0}; i < 1000; i++) {
+      mat3 = mat1 * mat2;
+    }
+  }
+
+  SECTION("Scalar Multiplication") {
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat3 = mat1 * 3;
+    }
+  }
+
+  SECTION("Element Multiply") {
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat1.ElementMultiply(mat2, mat3);
+    }
+  }
+
+  SECTION("Element Divide") {
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat1.ElementDivide(mat2, mat3);
+    }
+  }
+
+  SECTION("Minor Matrix") {
+    // what about matrices of 0,0 or 1,1?
+    // minor matrix for 2x2 matrix
+    Matrix<49, 49> minorMat1{};
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat1.MinorMatrix(minorMat1, 0, 0);
+    }
+  }
+
+  SECTION("Determinant") {
+    for (uint32_t i{0}; i < 1000000; i++) {
+      float det = mat4.Det();
+      (void)det;
+    }
+  }
+
+  SECTION("Matrix of Minors") {
+    for (uint32_t i{0}; i < 1000000; i++) {
+      mat4.MatrixOfMinors(mat5);
+    }
+  }
+
+  SECTION("Invert") {
+    for (uint32_t i{0}; i < 1000000; i++) {
+      mat5 = mat4.Invert();
+    }
+  };
+
+  SECTION("Transpose") {
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat3 = mat1.Transpose();
+    }
+  }
+
+  SECTION("Normalize") {
+    for (uint32_t i{0}; i < 100000; i++) {
+      mat3 = mat1 / mat1.EuclideanNorm();
+    }
+  }
+
+  SECTION("GET ROW") {
+    Matrix<1, 50> mat1Rows{};
+    for (uint32_t i{0}; i < 100000000; i++) {
+      mat1.GetRow(0, mat1Rows);
+    }
+  }
+
+  SECTION("GET COLUMN") {
+    Matrix<50, 1> mat1Columns{};
+    for (uint32_t i{0}; i < 100000000; i++) {
+      mat1.GetColumn(0, mat1Columns);
+    }
+  }
+
+  SECTION("QR Decomposition") {
+    Matrix<50, 50> Q, R{};
+    for (uint32_t i{0}; i < 500; i++) {
+      mat1.QRDecomposition(Q, R);
+    }
+  }
+}

unit-tests/quaternion-tests.cpp

@@ -0,0 +1,103 @@
+// 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 "Quaternion.h"
+
+// any other libraries
+#include <array>
+#include <cmath>
+#include <iostream>
+
+TEST_CASE("Vector Math", "Vector")
+{
+    Quaternion q1{1, 2, 3, 4};
+    Quaternion q2{5, 6, 7, 8};
+
+    SECTION("Initialization")
+    {
+        // explicit initialization
+        REQUIRE(q1.w == 1);
+        REQUIRE(q1.v1 == 2);
+        REQUIRE(q1.v2 == 3);
+        REQUIRE(q1.v3 == 4);
+
+        // fill initialization
+        Quaternion q3{0};
+        REQUIRE(q3.w == 0);
+        REQUIRE(q3.v1 == 0);
+        REQUIRE(q3.v2 == 0);
+        REQUIRE(q3.v3 == 0);
+
+        // copy initialization
+        Quaternion q4{q1};
+        REQUIRE(q4.w == 1);
+        REQUIRE(q4.v1 == 2);
+        REQUIRE(q4.v2 == 3);
+        REQUIRE(q4.v3 == 4);
+
+        // matrix initialization
+        Matrix<1, 4> m1{1, 2, 3, 4};
+        Quaternion q5{m1};
+        REQUIRE(q5.w == 1);
+        REQUIRE(q5.v1 == 2);
+        REQUIRE(q5.v2 == 3);
+        REQUIRE(q5.v3 == 4);
+
+        // array initialization
+        Quaternion q6{std::array<float, 4>{1, 2, 3, 4}};
+        REQUIRE(q6.w == 1);
+        REQUIRE(q6.v1 == 2);
+        REQUIRE(q6.v2 == 3);
+        REQUIRE(q6.v3 == 4);
+    }
+
+    SECTION("Equals")
+    {
+        Quaternion q3{0, 0, 0, 0};
+        q3 = q1;
+        REQUIRE(q3.w == 1);
+        REQUIRE(q3.v1 == 2);
+        REQUIRE(q3.v2 == 3);
+        REQUIRE(q3.v3 == 4);
+    }
+
+    SECTION("Array access")
+    {
+        REQUIRE(q1[0] == 1);
+        REQUIRE(q1[1] == 2);
+        REQUIRE(q1[2] == 3);
+        REQUIRE(q1[3] == 4);
+    }
+
+    SECTION("Addition")
+    {
+        Quaternion q3 = q1 + q2;
+        REQUIRE(q3.w == 6);
+        REQUIRE(q3.v1 == 8);
+        REQUIRE(q3.v2 == 10);
+        REQUIRE(q3.v3 == 12);
+    }
+
+    SECTION("Multiplication")
+    {
+        Quaternion q3;
+        q1.Q_Mult(q2, q3);
+        REQUIRE(q3.w == -60);
+        REQUIRE(q3.v1 == 12);
+        REQUIRE(q3.v2 == 30);
+        REQUIRE(q3.v3 == 24);
+    }
+
+    SECTION("Rotation")
+    {
+        Quaternion q3{Quaternion::FromAngleAndAxis(M_PI / 2, Matrix<1, 3>{0, 0, 1})};
+        Quaternion q4{0, 1, 0, 0};
+        Quaternion q5;
+        q3.Rotate(q4, q5);
+        REQUIRE_THAT(q5.v1, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+        REQUIRE_THAT(q5.v2, Catch::Matchers::WithinRel(1.0f, 1e-6f));
+        REQUIRE_THAT(q5.v3, Catch::Matchers::WithinRel(0.0f, 1e-6f));
+    }
+}

unit-tests/run-matrix-timing-test.sh

@@ -1,7 +0,0 @@
-# be in the root folder of this project when you run this
-cd build/
-ninja matrix-tests
-echo "Running tests. This will take a while."
-./unit-tests/matrix-tests -n "Timing Tests" -d yes > ../unit-tests/matrix-test-timings-temp.txt
-cd ../unit-tests/
-python3 test-timing-post-process.py

unit-tests/svd-integration-test.cpp

@@ -0,0 +1,252 @@
+#include "Matrix.hpp"
+#include "SVD.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <iostream>
+
+// Generic helper functions for any matrix size
+template <uint8_t rows, uint8_t columns>
+static float frobeniusNorm(const Matrix<rows, columns> &M) {
+  float sum = 0.0f;
+  for (int i = 0; i < rows; i++)
+    for (int j = 0; j < columns; j++) {
+      float v = M.Get(i, j);
+      sum += v * v;
+    }
+  return sqrtf(sum);
+}
+
+template <uint8_t n>
+static bool isOrthogonal(const Matrix<n, n> &M, float tol = 1e-4f) {
+  Matrix<n, n> Mt = M.Transpose();
+  Matrix<n, n> MtM{0};
+  Mt.Mult(M, MtM);
+  for (int i = 0; i < n; i++)
+    for (int 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;
+}
+
+TEST_CASE("SVD Integration: 2x2 [[1,2],[3,4]]", "[Matrix][SVD][Integration]") {
+  Matrix<2, 2> A{1, 2, 3, 4};
+  Matrix<2, 2> U{0};
+  Matrix<2, 1> sigma{0};
+  Matrix<2, 2> Vt{0};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // Reference singular values from scipy: [5.464985704219, 0.365966190626]
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(5.4649857f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(0.3659662f, 1e-3f));
+
+  // Check orthogonality of U and Vt (first 2x2 blocks)
+  REQUIRE(isOrthogonal<2>(U));
+  REQUIRE(isOrthogonal<2>(Vt));
+
+  // Check reconstruction: A ≈ U · diag(sigma) · Vt
+  Matrix<2, 2> recon{0};
+  Matrix<2, 2> Usig{0};
+  for (int i = 0; i < 2; i++)
+    for (int j = 0; j < 2; 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 < 2; i++)
+    for (int j = 0; j < 2; j++) {
+      float diff = recon.Get(i, j) - A.Get(i, j);
+      err += diff * diff;
+    }
+  err = sqrtf(err);
+  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)
+            << "]\n";
+}
+
+TEST_CASE("SVD Integration: 3x3 diagonal [10,5,2]",
+          "[Matrix][SVD][Integration]") {
+  Matrix<3, 3> A{10, 0, 0, 0, 5, 0, 0, 0, 2};
+  Matrix<3, 3> U{0};
+  Matrix<3, 1> sigma{0};
+  Matrix<3, 3> Vt{0};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // Singular values should be [10, 5, 2] (already diagonal)
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(10.0f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(5.0f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinRel(2.0f, 1e-3f));
+
+  // 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::WithinAbs(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]]",
+          "[Matrix][SVD][Integration]") {
+  Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
+  Matrix<3, 3> U{0};
+  Matrix<3, 1> sigma{0};
+  Matrix<3, 3> Vt{0};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // Reference: [16.848103352614, 1.068369514555, 0.0]
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(16.8481f, 1e-2f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.06837f, 1e-2f));
+  // Third singular value should be ~0 (rank-deficient)
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-2f));
+
+  // Check reconstruction
+  Matrix<3, 3> recon{0};
+  Matrix<3, 3> Usig{0};
+  for (int i = 0; i < 3; i++)
+    for (int j = 0; j < 3; 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 < 3; i++)
+    for (int j = 0; j < 3; j++) {
+      float diff = recon.Get(i, j) - A.Get(i, j);
+      err += diff * diff;
+    }
+  err = sqrtf(err);
+  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) << ", "
+            << sigma.Get(2, 0) << "]\n";
+}
+
+TEST_CASE("SVD Integration: tall 4x3 matrix", "[Matrix][SVD][Integration]") {
+  Matrix<4, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
+  Matrix<4, 3> U{0};
+  Matrix<3, 1> sigma{0};
+  Matrix<3, 3> Vt{0};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // Reference: [25.462407436036, 1.290661675761, 0.0]
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(25.4624f, 1e-2f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.29066f, 1e-2f));
+  REQUIRE_THAT(sigma.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-2f));
+
+  // Check reconstruction
+  Matrix<4, 3> recon{0};
+  Matrix<4, 3> Usig{0};
+  for (int i = 0; i < 4; i++)
+    for (int j = 0; j < 3; 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 < 4; i++)
+    for (int j = 0; j < 3; j++) {
+      float diff = recon.Get(i, j) - A.Get(i, j);
+      err += diff * diff;
+    }
+  err = sqrtf(err);
+  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) << ", "
+            << sigma.Get(2, 0) << "]\n";
+}
+
+TEST_CASE("SVD Integration: wide 3x5 matrix", "[Matrix][SVD][Integration]") {
+  Matrix<3, 5> A{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
+  Matrix<3, 5> U{0};
+  Matrix<5, 1> sigma{0}; // sigma is columns x 1 = 5x1 for wide matrix
+  Matrix<5, 5> Vt{0};    // Vt is columns x columns = 5x5
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // Reference: [35.127223333575, 2.465396696917, 0.0]
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(35.1272f, 1e-2f));
+  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: A (3x5) = U * Sigma * Vt, where U (3x5) has
+  // its meaningful part in the first 3 columns, sigma (5x1) in the
+  // first 3 entries, and Vt (5x5) in its first 3 rows (right
+  // singular vectors as rows). So:
+  //   A[i][j] = sum_k U[i][k] * sigma[k] * Vt[k][j]
+
+  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(k, j);
+      }
+      float diff = recon_val - A.Get(i, j);
+      err2 += diff * diff;
+    }
+  }
+  err2 = sqrtf(err2);
+  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) << ", "
+            << sigma.Get(2, 0) << "]\n";
+}
+
+TEST_CASE("SVD Integration: identity 3x3", "[Matrix][SVD][Integration]") {
+  Matrix<3, 3> A{1, 0, 0, 0, 1, 0, 0, 0, 1};
+  Matrix<3, 3> U{0};
+  Matrix<3, 1> sigma{0};
+  Matrix<3, 3> Vt{0};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(1.0f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(1.0f, 1e-3f));
+  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::WithinAbs(0.0f, 1e-2f));
+}
+
+TEST_CASE("SVD Integration: symmetric positive definite 2x2 [[5,3],[3,5]]",
+          "[Matrix][SVD][Integration]") {
+  Matrix<2, 2> A{5, 3, 3, 5};
+  Matrix<2, 2> U{0};
+  Matrix<2, 1> sigma{0};
+  Matrix<2, 2> Vt{0};
+
+  SVD::SVD(A, U, sigma, Vt);
+
+  // For SPD matrix, singular values = eigenvalues: [8, 2]
+  REQUIRE_THAT(sigma.Get(0, 0), Catch::Matchers::WithinRel(8.0f, 1e-3f));
+  REQUIRE_THAT(sigma.Get(1, 0), Catch::Matchers::WithinRel(2.0f, 1e-3f));
+
+  // Check reconstruction
+  Matrix<2, 2> recon{0};
+  Matrix<2, 2> Usig{0};
+  for (int i = 0; i < 2; i++)
+    for (int j = 0; j < 2; 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 < 2; i++)
+    for (int j = 0; j < 2; j++) {
+      float diff = recon.Get(i, j) - A.Get(i, j);
+      err += diff * diff;
+    }
+  err = sqrtf(err);
+  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";
+}

unit-tests/svd-reference-values.py

@@ -0,0 +1,482 @@
+#!/usr/bin/env python3
+"""
+Generate reference values for SVD building block unit tests.
+Run this to verify/implement the C++ SVD implementation against scipy/numpy.
+
+Usage: python3 svd-reference-values.py
+"""
+
+import numpy as np
+from scipy.linalg import svd, qr as scipy_qr
+import json
+
+def compute_householder(x):
+    """Compute Householder reflector: H*x = [alpha, 0, 0, ...]^T.
+    
+    Returns (v_normalized, alpha) where v is the normalized Householder vector.
+    H = I - 2*v*v^T / (v^T*v)
+    """
+    x = np.array(x, dtype=np.float64)
+    norm_x = np.linalg.norm(x)
+    
+    if norm_x < 1e-30:
+        return x.copy(), 0.0
+    
+    alpha = -np.sign(x[0]) * norm_x if x[0] != 0 else -norm_x
+    v = x.copy()
+    v[0] -= alpha
+    v_norm = np.linalg.norm(v)
+    
+    if v_norm < 1e-30:
+        return np.zeros_like(x), alpha
+    
+    v /= v_norm
+    return v, alpha
+
+def apply_householder_left(A, v, start_row):
+    """Apply Householder reflection from the left: A = (I - 2vv^T) @ A.
+    
+    v is the normalized Householder vector operating on rows [start_row:].
+    The length of v must match the number of rows affected.
+    """
+    A = A.copy()
+    k = len(v)
+    
+    for col in range(A.shape[1]):
+        dot = np.dot(v, A[start_row:start_row+k, col])
+        A[start_row:start_row+k, col] -= 2.0 * dot * v
+    
+    return A
+
+def apply_householder_right(A, v, start_col):
+    """Apply Householder reflection from the right: A = A @ (I - 2vv^T).
+    
+    v is the normalized Householder vector operating on columns [start_col:].
+    The length of v must match the number of columns affected.
+    """
+    A = A.copy()
+    k = len(v)
+    
+    for row in range(A.shape[0]):
+        dot = np.dot(A[row, start_col:start_col+k], v)
+        A[row, start_col:start_col+k] -= 2.0 * dot * v
+    
+    return A
+
+def compute_givens(x, y):
+    """Compute Givens rotation that zeros out y.
+    
+    Returns (c, s) such that [c s; -s c] @ [x; y] = [r; 0].
+    """
+    r = np.sqrt(x*x + y*y)
+    if r < 1e-30:
+        return 1.0, 0.0
+    c = x / r
+    s = y / r
+    return c, s
+
+def apply_givens_left(A, i, j, c, s):
+    """Apply Givens rotation from the left to rows i and j of A.
+    
+    [c  s]   [row_i]
+    [-s  c] @ [row_j] = [new_row_i]
+                          [new_row_j]
+    """
+    A = A.copy()
+    new_i = c * A[i] + s * A[j]
+    new_j = -s * A[i] + c * A[j]
+    A[i] = new_i
+    A[j] = new_j
+    return A
+
+def apply_givens_right(A, i, j, c, s):
+    """Apply Givens rotation from the right to columns i and j of A.
+    
+    [col_i  col_j] @ [c -s] = [new_col_i  new_col_j]
+                       [s  c]
+    """
+    A = A.copy()
+    new_i = c * A[:, i] + s * A[:, j]
+    new_j = -s * A[:, i] + c * A[:, j]
+    A[:, i] = new_i
+    A[:, j] = new_j
+    return A
+
+def householder_bidiagonalization(A):
+    """Full Householder bidiagonalization: A = Q_L @ B @ Q_R^T.
+    
+    Returns (B, Q_L, Q_R) where B is upper bidiagonal.
+    """
+    m, n = A.shape
+    p = min(m, n)
+    
+    QL = np.eye(m, dtype=np.float64)
+    QR = np.eye(n, dtype=np.float64)
+    W = A.copy()
+    
+    for k in range(p):
+        # Left HH: zero out W[k+1:, k]
+        if k < m - 1:
+            x = W[k+1:, k].copy()
+            v, alpha = compute_householder(x)
+            if np.linalg.norm(v) > 1e-30:
+                W = apply_householder_left(W, v, k + 1)
+                QL = apply_householder_right(QL, v, k + 1)
+        
+        # Right HH: zero out W[k, k+2:] (superdiagonal)
+        if k < p - 1 and k + 2 <= n:
+            x = W[k, k+2:].copy()
+            v, alpha = compute_householder(x)
+            if np.linalg.norm(v) > 1e-30:
+                W = apply_householder_right(W, v, k + 2)
+                QR = apply_householder_right(QR, v, k + 2)
+    
+    return W, QL, QR
+
+def implicit_qr_iteration(B, QR_acc):
+    """Implicit QR iteration on a bidiagonal matrix.
+    
+    Returns (Sigma, QR_acc) where Sigma is diagonal with singular values
+    and QR_acc contains the accumulated right transformations.
+    """
+    m, n = B.shape
+    p = min(m, n)
+    W = B.copy()
+    
+    max_iter = 1000
+    tol = 1e-10
+    
+    for iteration in range(max_iter):
+        # Deflate negligible subdiagonal elements
+        for i in range(p - 1, 0, -1):
+            if abs(W[i, i-1]) < tol * (abs(W[i-1, i-1]) + abs(W[i, i])):
+                W[i, i-1] = 0.0
+        
+        # Find smallest unreduced block [start, end]
+        start = 0
+        for i in range(p - 1):
+            if abs(W[i+1, i]) >= tol * (abs(W[i, i]) + abs(W[i+1, i+1])):
+                start = i + 1
+        
+        end = p - 1
+        for i in range(p - 2, -1, -1):
+            if abs(W[i+1, i]) >= tol * (abs(W[i, i]) + abs(W[i+1, i+1])):
+                end = i
+                break
+        
+        if start >= end:
+            continue
+        
+        # Wilkinson shift from bottom 2x2 corner
+        a, b = W[end-1, end-1], W[end-1, end]
+        c_val, d = W[end, end-1], W[end, end]
+        trace = a + d
+        det = a * d - b * c_val
+        disc = trace**2 - 4 * det
+        
+        if disc >= 0:
+            sqrt_disc = np.sqrt(disc)
+            e1, e2 = (trace + sqrt_disc) / 2, (trace - sqrt_disc) / 2
+            shift = e1 if abs(e1 - d) < abs(e2 - d) else e2
+        else:
+            shift = d
+        
+        # Implicit QR step using Givens rotations
+        # Process from top to bottom within the block
+        x = W[start, start] - shift
+        y = W[start + 1, start]
+        
+        for i in range(start, end):
+            r = np.sqrt(x*x + y*y)
+            if r < 1e-30:
+                x = W[i + 1, i]
+                y = W[i + 1, i + 1] if i + 2 <= end else 0.0
+                continue
+            
+            c_rot = x / r
+            s_rot = y / r
+            
+            # Apply from left to rows i, i+1 (columns i..n-1)
+            for j in range(i, n):
+                t1, t2 = W[i, j], W[i + 1, j]
+                W[i, j] = c_rot * t1 + s_rot * t2
+                W[i + 1, j] = -s_rot * t1 + c_rot * t2
+            
+            # Apply from right to columns i, i+1 (rows 0..i)
+            if i > start:
+                for j in range(i + 1):
+                    t1, t2 = W[j, i], W[j, i + 1]
+                    W[j, i] = c_rot * t1 + s_rot * t2
+                    W[j, i + 1] = -s_rot * t1 + c_rot * t2
+            
+            # Accumulate into QR_acc
+            for j in range(QR_acc.shape[0]):
+                t1, t2 = QR_acc[j, i], QR_acc[j, i + 1]
+                QR_acc[j, i] = c_rot * t1 + s_rot * t2
+                QR_acc[j, i + 1] = -s_rot * t1 + c_rot * t2
+            
+            # Prepare for next rotation
+            x = W[i + 1, i]
+            y = W[i + 1, i + 1] if i + 2 <= end else 0.0
+    
+    return W, QR_acc
+
+
+def main():
+    print("=" * 70)
+    print("SVB BUILDING BLOCK REFERENCE VALUES")
+    print("Generated with scipy/numpy for C++ unit test verification")
+    print("=" * 70)
+    
+    # ------------------------------------------------------------------
+    # Test 1: Householder Vector Computation
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 1: computeHouseholderVector")
+    print("=" * 70)
+    
+    test_vectors = [
+        ("2D [1,3]", [1.0, 3.0]),
+        ("2D [3,4] (norm=5)", [3.0, 4.0]),
+        ("3D [1,2,3]", [1.0, 2.0, 3.0]),
+        ("3D [0,0,1]", [0.0, 0.0, 1.0]),
+        ("4D [5,-3,2,1]", [5.0, -3.0, 2.0, 1.0]),
+    ]
+    
+    for name, vec in test_vectors:
+        v, alpha = compute_householder(vec)
+        x = np.array(vec)
+        Hx = x - 2 * np.dot(v, x) * v
+        
+        print(f"\n{name}:")
+        print(f"  Input:            {list(x)}")
+        print(f"  ||x||:            {np.linalg.norm(x):.15f}")
+        print(f"  alpha:            {alpha:.15f}")
+        print(f"  v (normalized):   {[round(float(vi), 12) for vi in v]}")
+        print(f"  H*x = [alpha,0..]: {[round(float(xi), 12) for xi in Hx]}")
+        print(f"  Off-diagonal ~0:  {np.allclose(Hx[1:], 0, atol=1e-12)}")
+    
+    # ------------------------------------------------------------------
+    # Test 2: Householder Apply Left
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 2: applyHouseholderLeft")
+    print("=" * 70)
+    
+    A_test = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], dtype=np.float64)
+    x_col = A_test[1:, 0].copy()
+    v_left, _ = compute_householder(x_col)
+    
+    print(f"\nInput matrix:\n{A_test}")
+    print(f"Householder vector (rows 1:3): {[round(float(vi), 12) for vi in v_left]}")
+    
+    A_result = apply_householder_left(A_test, v_left, 1)
+    print(f"\nAfter applyHouseholderLeft:\n{A_result}")
+    print(f"  A[1,0] = {A_result[1,0]:.2e}, A[2,0] = {A_result[2,0]:.2e} (should be ~0)")
+    
+    # ------------------------------------------------------------------
+    # Test 3: Householder Apply Right
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 3: applyHouseholderRight")
+    print("=" * 70)
+    
+    A_test = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], dtype=np.float64)
+    x_row = A_test[0, 1:].copy()
+    v_right, _ = compute_householder(x_row)
+    
+    print(f"\nInput matrix:\n{A_test}")
+    print(f"Householder vector (cols 1:3): {[round(float(vi), 12) for vi in v_right]}")
+    
+    A_result = apply_householder_right(A_test, v_right, 1)
+    print(f"\nAfter applyHouseholderRight:\n{A_result}")
+    print(f"  A[0,1] = {A_result[0,1]:.2e}, A[0,2] = {A_result[0,2]:.2e} (should be ~0)")
+    
+    # ------------------------------------------------------------------
+    # Test 4: Givens Rotation Computation
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 4: computeGivens")
+    print("=" * 70)
+    
+    givens_tests = [
+        ("3-4-5 triangle", 3.0, 4.0),
+        ("y already zero", 1.0, 0.0),
+        ("x is zero", 0.0, 5.0),
+        ("Both negative", -3.0, -4.0),
+        ("45 degree case", 1.0, -1.0),
+    ]
+    
+    for name, x, y in givens_tests:
+        c, s = compute_givens(x, y)
+        result_x = c * x + s * y
+        result_y = -s * x + c * y
+        
+        print(f"\n{name}: x={x}, y={y}")
+        print(f"  r = {np.sqrt(x*x+y*y):.12f}")
+        print(f"  c = {c:.12f}, s = {s:.12f}")
+        print(f"  [c s; -s c] @ [x;y] = [{result_x:.2e}, {result_y:.2e}]")
+    
+    # ------------------------------------------------------------------
+    # Test 5: Apply Givens Left/Right
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 5: applyGivensLeft / applyGivensRight")
+    print("=" * 70)
+    
+    A_test = np.array([[3.0, 4.0], [1.0, 2.0]], dtype=np.float64)
+    c, s = compute_givens(3.0, 1.0)
+    
+    print(f"\nInput matrix:\n{A_test}")
+    print(f"Givens rotation (rows 0,1): c={c:.12f}, s={s:.12f}")
+    
+    A_left = apply_givens_left(A_test, 0, 1, c, s)
+    print(f"\nAfter applyGivensLeft:\n{A_left}")
+    print(f"  A[1,0] = {A_left[1,0]:.2e} (should be ~0)")
+    
+    A_test = np.array([[3.0, 1.0], [4.0, 2.0]], dtype=np.float64)
+    c, s = compute_givens(3.0, 4.0)
+    
+    print(f"\nInput matrix:\n{A_test}")
+    print(f"Givens rotation (cols 0,1): c={c:.12f}, s={s:.12f}")
+    
+    A_right = apply_givens_right(A_test, 0, 1, c, s)
+    print(f"\nAfter applyGivensRight:\n{A_right}")
+    print(f"  A[0,1] = {A_right[0,1]:.2e} (should be ~0)")
+    
+    # ------------------------------------------------------------------
+    # Test 6: Full Bidiagonalization
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 6: householderBidiagonalization")
+    print("=" * 70)
+    
+    bidiag_tests = [
+        ("2x2 [[1,2],[3,4]]", np.array([[1.0, 2.0], [3.0, 4.0]])),
+        ("3x3 SPD [[5,3],[3,5]]", np.array([[5.0, 3.0], [3.0, 5.0]])),
+        ("3x3 diag [[10,0,0],[0,5,0],[0,0,2]]", 
+         np.array([[10.0, 0, 0], [0, 5.0, 0], [0, 0, 2.0]])),
+        ("3x3 full [[1,2,3],[4,5,6],[7,8,10]]", 
+         np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 10.0]])),
+        ("Tall 4x3", np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], dtype=np.float64)),
+    ]
+    
+    for name, A in bidiag_tests:
+        B, QL, QR = householder_bidiagonalization(A)
+        m, n = A.shape
+        p = min(m, n)
+        
+        print(f"\n{name}:")
+        print(f"  Original:\n{A}")
+        print(f"\n  Bidiagonal B:\n{B}")
+        print(f"  Diagonal: {[round(float(B[i,i]), 10) for i in range(p)]}")
+        print(f"  Superdiag: {[round(float(B[i,i+1]), 10) for i in range(min(p-1, n-1))]}")
+        
+        recon = QL @ B @ QR.T
+        err = np.linalg.norm(recon - A, 'fro')
+        print(f"  ||QL @ B @ QR^T - A||_F = {err:.2e}")
+    
+    # ------------------------------------------------------------------
+    # Test 7: Full SVD Reference Values
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 7: Full SVD Reference Values (scipy.linalg.svd)")
+    print("=" * 70)
+    
+    test_matrices = [
+        ("Simple 2x2", np.array([[1,2],[3,4]], dtype=np.float64)),
+        ("SPD 2x2", np.array([[5,3],[3,5]], dtype=np.float64)),
+        ("Full-rank 3x3", np.array([[1,2,3],[4,5,6],[7,8,10]], dtype=np.float64)),
+        ("Rank-deficient 3x3", np.array([[1,2,3],[4,5,6],[7,8,9]], dtype=np.float64)),
+        ("Diagonal 3x3", np.array([[10,0,0],[0,5,0],[0,0,2]], dtype=np.float64)),
+        ("Tall 4x3", np.array([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], dtype=np.float64)),
+        ("Wide 3x5", np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15]], dtype=np.float64)),
+        ("Symmetric tri 5x5", np.array([[2,-1,0,0,0],[-1,2,-1,0,0],[0,-1,2,-1,0],[0,0,-1,2,-1],[0,0,0,-1,2]], dtype=np.float64)),
+        ("Neg values 2x3", np.array([[0.5,-0.3,0.8],[-0.2,0.7,0.1]], dtype=np.float64)),
+        ("Near-singular 2x2", np.array([[1,0],[0,1e-6]], dtype=np.float64)),
+        ("Orthogonal 3x3", np.array([[np.cos(np.pi/4), -np.sin(np.pi/4), 0],
+                                      [np.sin(np.pi/4), np.cos(np.pi/4), 0],
+                                      [0, 0, 1]], dtype=np.float64)),
+        ("Identity 3x3", np.eye(3)),
+        ("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)),
+    ]
+    
+    for name, A in test_matrices:
+        U, s, Vt = svd(A, full_matrices=False)
+        
+        print(f"\n{name}: shape={A.shape}")
+        print(f"  Singular values: {[round(float(x), 12) for x in s]}")
+        print(f"  U:\n{np.array2string(U, precision=6, floatmode='maxprec_equal')}")
+        print(f"  Vt:\n{np.array2string(Vt, precision=6, floatmode='maxprec_equal')}")
+        recon_err = np.linalg.norm(A - U @ np.diag(s) @ Vt, 'fro')
+        print(f"  Reconstruction error: {recon_err:.2e}")
+    
+    # ------------------------------------------------------------------
+    # Test 8: Implicit QR Iteration on Bidiagonal
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("TEST 8: implicitQRIteration")
+    print("=" * 70)
+    
+    qr_tests = [
+        ("2x2 [[1,2],[3,4]]", np.array([[1.0, 2.0], [3.0, 4.0]])),
+        ("3x3 diag", np.array([[10.0, 0, 0], [0, 5.0, 0], [0, 0, 2.0]])),
+    ]
+    
+    for name, A in qr_tests:
+        B, QL, QR = householder_bidiagonalization(A)
+        Sigma, QR_final = implicit_qr_iteration(B.copy(), QR.copy())
+        
+        print(f"\n{name}:")
+        print(f"  Bidiagonal B:\n{B}")
+        print(f"  After QR iteration (Sigma):\n{Sigma}")
+        print(f"  Diagonal entries: {[round(float(Sigma[i,i]), 10) for i in range(min(Sigma.shape))]}")
+        
+        # Verify: QL @ Sigma @ QR_final^T ≈ A
+        recon = QL @ Sigma @ QR_final.T
+        err = np.linalg.norm(recon - A, 'fro')
+        print(f"  ||QL @ Sigma @ QR^T - A||_F = {err:.2e}")
+    
+    # ------------------------------------------------------------------
+    # JSON output for easy import into C++ tests
+    # ------------------------------------------------------------------
+    print("\n" + "=" * 70)
+    print("JSON OUTPUT (for easy C++ integration)")
+    print("=" * 70)
+    
+    json_data = {}
+    
+    # Householder test vectors
+    hh_tests = {}
+    for name, vec in test_vectors:
+        v, alpha = compute_householder(vec)
+        x = np.array(vec)
+        Hx = x - 2 * np.dot(v, x) * v
+        hh_tests[name] = {
+            "input": [float(xi) for xi in x],
+            "norm": float(np.linalg.norm(x)),
+            "alpha": float(alpha),
+            "v_normalized": [round(float(vi), 12) for vi in v],
+            "Hx": [round(float(xi), 12) for xi in Hx],
+        }
+    json_data["householder_vectors"] = hh_tests
+    
+    # Full SVD reference values
+    svd_tests = {}
+    for name, A in test_matrices:
+        U, s, Vt = svd(A, full_matrices=False)
+        svd_tests[name] = {
+            "shape": list(A.shape),
+            "singular_values": [round(float(x), 12) for x in s],
+            "U": [[round(float(U[i,j]), 8) for j in range(U.shape[1])] for i in range(U.shape[0])],
+            "Vt": [[round(float(Vt[i,j]), 8) for j in range(Vt.shape[1])] for i in range(Vt.shape[0])],
+        }
+    json_data["svd_reference"] = svd_tests
+    
+    print(json.dumps(json_data, indent=2))
+
+
+if __name__ == "__main__":
+    main()

unit-tests/timing-results/matrix-timing-tests.txt

@@ -0,0 +1,36 @@
+Running matrix-timing-tests with timing
+Randomness seeded to: 3567651885
+1.857 s: Addition
+1.857 s: Timing Tests
+1.788 s: Subtraction
+1.788 s: Timing Tests
+1.929 s: Multiplication
+1.929 s: Timing Tests
+1.268 s: Scalar Multiplication
+1.268 s: Timing Tests
+1.798 s: Element Multiply
+1.798 s: Timing Tests
+1.802 s: Element Divide
+1.803 s: Timing Tests
+1.553 s: Minor Matrix
+1.554 s: Timing Tests
+1.009 s: Determinant
+1.009 s: Timing Tests
+4.076 s: Matrix of Minors
+4.076 s: Timing Tests
+1.066 s: Invert
+1.066 s: Timing Tests
+1.246 s: Transpose
+1.246 s: Timing Tests
+2.284 s: Normalize
+2.284 s: Timing Tests
+0.606 s: GET ROW
+0.606 s: Timing Tests
+24.629 s: GET COLUMN
+24.630 s: Timing Tests
+3.064 s: QR Decomposition
+3.064 s: Timing Tests
+===============================================================================
+test cases: 1 | 1 passed
+assertions: - none -
+

unit-tests/vector-tests.cpp

@@ -0,0 +1,45 @@
+// 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 "Vector3D.hpp"
+#include "Matrix.hpp"
+
+// any other libraries
+#include <array>
+#include <cmath>
+#include <iostream>
+
+TEST_CASE("Vector Math", "Vector")
+{
+    V3D<float> v1{1, 2, 3};
+    V3D<float> v2{4, 5, 6};
+    V3D<float> v3{};
+
+    SECTION("Initialization")
+    {
+        // list initialization
+        REQUIRE(v1.x == 1);
+        REQUIRE(v1.y == 2);
+        REQUIRE(v1.z == 3);
+
+        // copy initialization
+        V3D<float> v4{v2};
+        REQUIRE(v4.x == 4);
+        REQUIRE(v4.y == 5);
+        REQUIRE(v4.z == 6);
+
+        // empty initialization
+        REQUIRE(v3.x == 0);
+        REQUIRE(v3.y == 0);
+        REQUIRE(v3.z == 0);
+
+        // matrix initialization
+        Matrix<1, 3> mat1{v1.ToArray()};
+        V3D<float> v5{mat1};
+        REQUIRE(v5.x == v1.x);
+        REQUIRE(v5.y == v1.y);
+        REQUIRE(v5.z == v1.z);
+    }
+}