Working on getting the QR decomposition to compile

.gitea/workflows/Merge-Checker.yaml

@@ -1,102 +0,0 @@
-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/
-.cache/
+venv/

.vscode/launch.json

@@ -8,14 +8,14 @@
             "name": "Debug Matrix Unit Tests",
             "type": "cppdbg",
             "request": "launch",
-            "program": "${workspaceFolder}/build/unit-tests/matrix-tests",
+            "program": "${workspaceFolder}/build/unit-tests/matrix-tests", 
-            "args": [],
+            "args": [], 
             "stopAtEntry": false,
             "cwd": "${workspaceFolder}",
             "environment": [],
             "externalConsole": false,
-            "MIMode": "gdb",
+            "MIMode": "gdb", 
-            "miDebuggerPath": "/usr/bin/gdb", // Adjust to your debugger path
+            "miDebuggerPath": "/usr/bin/gdb",  // Adjust to your debugger path
             "setupCommands": [
                 {
                     "description": "Enable pretty-printing for gdb",
@@ -23,29 +23,20 @@
                     "ignoreFailures": true
                 }
             ],
-            "preLaunchTask": "build_tests", // Task to compile unit tests
+            "preLaunchTask": "build_tests",  // Task to compile unit tests
             "internalConsoleOptions": "openOnSessionStart"
         },
         {
-            "name": "Debug Quaternion Unit Tests",
+            "name": "Run Matrix Unit Tests",
-            "type": "cppdbg",
+            "type": "cpp",
             "request": "launch",
-            "program": "${workspaceFolder}/build/unit-tests/quaternion-tests",
+            "program": "${workspaceFolder}/build/unit-tests/matrix-tests",
-            "args": [],
+            "args": [], 
             "stopAtEntry": false,
             "cwd": "${workspaceFolder}",
             "environment": [],
             "externalConsole": false,
-            "MIMode": "gdb",
+            "preLaunchTask": "build_tests",  // Compile unit tests before running
-            "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,5 +1,8 @@
 {
     "C_Cpp.intelliSenseEngine": "default",
+    "clangd.arguments": [
+        "--include-directory=build/unit-tests"
+    ],
     "C_Cpp.default.intelliSenseMode": "linux-gcc-x64",
     "files.associations": {
         "*.h": "cpp",
@@ -68,16 +71,8 @@
         "typeinfo": "cpp",
         "variant": "cpp",
         "shared_mutex": "cpp",
-        "charconv": "cpp",
+        "complex": "cpp"
-        "format": "cpp",
-        "csignal": "cpp",
-        "span": "cpp"
     },
-    "clangd.enable": true,
+    "clangd.enable": false,
-    "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,14 +4,12 @@
         {
             "label": "build_tests",
             "type": "shell",
-            "command": "cd build && ninja",
+            "command": "cd build && ninja matrix-tests",
             "group": {
                 "kind": "build",
                 "isDefault": true
             },
-            "problemMatcher": [
+            "problemMatcher": ["$gcc"],
-                "$gcc"
-            ],
             "detail": "Generated task to build unit test executable"
         }
     ]

CMakeLists.txt

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

Matrix.cpp

@@ -1,10 +1,3 @@
-// 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"
@@ -12,45 +5,18 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdlib>
-#include <cstring>
+#include <type_traits>
+
+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++) {
@@ -61,6 +27,19 @@ 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) {
@@ -112,18 +91,21 @@ 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<columns, 1> other_column;
+  Matrix<rows, 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 < other_columns; column_idx++) {
+    for (uint8_t column_idx{0}; column_idx < 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.Transpose());
+          Matrix<rows, columns>::dotProduct(this_row, other_column_t);
     }
   }
 
@@ -143,13 +125,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>::Invert() const {
+Matrix<rows, columns> &
+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()};
@@ -178,8 +160,8 @@ Matrix<rows, columns> Matrix<rows, columns>::Invert() const {
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<columns, rows> Matrix<rows, columns>::Transpose() const {
+Matrix<columns, rows> &
-  Matrix<columns, rows> result{};
+Matrix<rows, columns>::Transpose(Matrix<columns, rows> &result) const {
   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);
@@ -191,10 +173,9 @@ Matrix<columns, rows> Matrix<rows, columns>::Transpose() const {
 
 // explicitly define the determinant for a 2x2 matrix because it is definitely
 // the fastest way to calculate a 2x2 matrix determinant
-// template <>
+template <> float Matrix<0, 0>::Det() const { return 1e+6; }
-// inline float Matrix<0, 0>::Det() const { return 1e+6; }
+template <> float Matrix<1, 1>::Det() const { return this->matrix[0]; }
-template <> inline float Matrix<1, 1>::Det() const { return this->matrix[0]; }
+template <> float Matrix<2, 2>::Det() const {
-template <> inline float Matrix<2, 2>::Det() const {
   return this->matrix[0] * this->matrix[3] - this->matrix[1] * this->matrix[2];
 }
 
@@ -290,13 +271,8 @@ void Matrix<rows, columns>::ToString(std::string &stringBuffer) const {
 }
 
 template <uint8_t rows, uint8_t columns>
-const float *Matrix<rows, columns>::ToArray() const {
+std::array<float, columns> &Matrix<rows, columns>::
-  return this->matrix.data();
+operator[](uint8_t row_index) {
-}
-
-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;
@@ -308,36 +284,38 @@ Matrix<rows, columns>::operator[](uint8_t row_index) {
 }
 
 template <uint8_t rows, uint8_t columns>
-Matrix<rows, columns> &
+Matrix<rows, columns> &Matrix<rows, columns>::
-Matrix<rows, columns>::operator=(const Matrix<rows, columns> &other) {
+operator=(const Matrix<rows, columns> &other) {
-  memcpy(this->matrix.begin(), other.matrix.begin(),
+  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
-         rows * columns * sizeof(float));
+    for (uint8_t column_idx{0}; column_idx < columns; column_idx++) {
-
+      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>::
-Matrix<rows, columns>::operator+(const Matrix<rows, columns> &other) const {
+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>::
-Matrix<rows, columns>::operator-(const Matrix<rows, columns> &other) const {
+operator-(const Matrix<rows, columns> &other) const {
   Matrix<rows, columns> buffer{};
   this->Sub(other, buffer);
   return buffer;
 }
 
 template <uint8_t rows, uint8_t columns>
-template <uint8_t other_columns>
+Matrix<rows, columns> Matrix<rows, columns>::
-Matrix<rows, other_columns> Matrix<rows, columns>::operator*(
+operator*(const Matrix<rows, columns> &other) const {
-    const Matrix<columns, other_columns> &other) const {
+  Matrix<rows, columns> buffer{};
-  Matrix<rows, other_columns> buffer{};
   this->Mult(other, buffer);
   return buffer;
 }
@@ -349,25 +327,9 @@ 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++) {
@@ -379,7 +341,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++) {
@@ -391,11 +353,7 @@ 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) {
-  for (uint8_t row_idx{0}; row_idx < rows; row_idx++) {
+  this->matrix.fill(value);
-    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>
@@ -451,8 +409,8 @@ Matrix<rows, columns>::adjugate(Matrix<rows, columns> &result) const {
 }
 
 template <uint8_t rows, uint8_t columns>
-float Matrix<rows, columns>::EuclideanNorm() const {
+Matrix<rows, columns> &
-
+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++) {
@@ -461,147 +419,90 @@ float Matrix<rows, columns>::EuclideanNorm() const {
     }
   }
 
-  return sqrt(sum);
+  if (sum == 0) {
+    // 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>
-template <uint8_t sub_rows, uint8_t sub_columns, uint8_t row_offset,
+Matrix<rows, rows> Matrix<rows, columns>::Eye() {
-          uint8_t column_offset>
+  Matrix<rows, rows> i_matrix;
-Matrix<sub_rows, sub_columns> Matrix<rows, columns>::SubMatrix() const {
+  i_matrix.Fill(0);
-  // static assert that sub_rows + row_offset <= rows
+  for (uint8_t i{0}; i < rows; i++) {
-  // static assert that sub_columns + column_offset <= columns
+    i_matrix[i][i] = 1;
-  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 buffer;
+  return i_matrix;
 }
 
 template <uint8_t rows, uint8_t columns>
-template <uint8_t sub_rows, uint8_t sub_columns>
+void Matrix<rows, columns>::QR_Decomposition(Matrix<rows, columns> &Q,
-void Matrix<rows, columns>::SetSubMatrix(
+                                             Matrix<rows, columns> &R) const {
-    uint8_t rowOffset, uint8_t columnOffset,
+  Q = Matrix<rows, columns>::Eye(); // Q starts as the identity matrix
-    const Matrix<sub_rows, sub_columns> &sub_matrix) {
+  R = *this; // R starts as a copy of this matrix (For this algorithm we'll call
-  int16_t adjustedSubRows = sub_rows;
+             // this matrix A)
-  int16_t adjustedSubColumns = sub_columns;
-  int16_t adjustedRowOffset = rowOffset;
-  int16_t adjustedColumnOffset = columnOffset;
 
-  // a bunch of safety checks to make sure we don't overflow the matrix
+  for (uint8_t row{0}; row < rows; row++) {
-  if (sub_rows > rows) {
+    // compute the householder vector
-    adjustedSubRows = rows;
+    const uint8_t houseHoldVectorSize{rows - row};
-  }
+    const uint8_t subMatrixSize{columns - row};
-  if (sub_columns > columns) {
+    Matrix<houseHoldVectorSize, 1> x{};
-    adjustedSubColumns = columns;
+    this->SubMatrix(row, row, x);
-  }
 
-  if (adjustedSubRows + adjustedRowOffset >= rows) {
+    Matrix<houseHoldVectorSize, 1> e1{};
-    adjustedRowOffset =
+    e1.Fill(0);
-        std::max(0, static_cast<int16_t>(rows) - adjustedSubRows);
+    if (x[0][0] >= 0) {
-  }
+      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 {
-      u.Fill(0);
+      e1[0][0] = -x.Norm();
     }
-    Q.SetSubMatrix(0, column, u);
 
-    // -----------------------
+    Matrix<houseHoldVectorSize, 1> v = x + e1;
-    // ----- CALCULATE R -----
+    v = v * (1 / v.Norm()); // normalize V
-    // -----------------------
+
-    for (uint8_t k = 0; k <= column; k++) {
+    // ************************************
-      Q.GetColumn(k, e);
+    // Apply the reflection to the R matrix
-      R[k][column] = (a_col.Transpose() * e).Get(0, 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);
+
+    // ************************************
+    // Apply the reflection to the Q matrix
+    // ************************************
+    // 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(
+                                        Q_subMatrix.Mult(v, Qv), v, Qv_outer);
+    Q.CopySubMatrixInto(0, row, Q_subMatrix);
   }
 }
 
-template <uint8_t rows, uint8_t columns>
-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
-  Matrix<rows, rows> QQ{Matrix<rows, rows>::Identity()};
-  Matrix<rows, rows> shift{0};
-
-  for (uint32_t iter = 0; iter < maxIterations; ++iter) {
-    Matrix<rows, rows> Q, R;
-
-    // // QR shift lets us "attack" the first diagonal to speed up the algorithm
-    // shift = Matrix<rows, rows>::Identity() * Ak[rows - 1][rows - 1];
-    (Ak - shift).QRDecomposition(Q, R);
-    Ak = R * Q + shift;
-    QQ = QQ * Q;
-
-    // Check convergence: off-diagonal norm
-    float offDiagSum = 0.0f;
-    for (uint32_t row = 1; row < rows; row++) {
-      for (uint32_t column = 0; column < row; column++) {
-        offDiagSum += fabs(Ak[row][column]);
-      }
-    }
-
-    if (offDiagSum < tolerance) {
-      break;
-    }
-  }
-
-  // Diagonal elements are the eigenvalues
-  for (uint8_t i = 0; i < rows; i++) {
-    eigenValues[i][0] = Ak[i][i];
-  }
-  eigenVectors = QQ;
-}
-
 #endif // MATRIX_H_

Matrix.hpp

@@ -1,9 +1,8 @@
-#pragma once
+#ifndef MATRIX_H_
+#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
@@ -12,13 +11,16 @@
 
 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
    */
@@ -30,17 +32,9 @@ public:
   Matrix(const Matrix<rows, columns> &other);
 
   /**
-   * @brief Initialize a matrix directly with scalar values
+   * @brief Initialize a matrix directly with any number of arguments
-   * Uses SFINAE to only accept arithmetic types (int, float, double, etc.)
    */
-  template <typename... Args,
+  template <typename... Args> Matrix(Args... 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
@@ -102,8 +96,8 @@ public:
                                        Matrix<rows, columns> &result) const;
 
   Matrix<rows - 1, columns - 1> &
-  MinorMatrix(Matrix<rows - 1, columns - 1> &result, uint8_t row_idx,
+      MinorMatrix(Matrix<rows - 1, columns - 1> &result, uint8_t row_idx,
-              uint8_t column_idx) const;
+                  uint8_t column_idx) const;
 
   /**
    * @return Get the determinant of the matrix
@@ -118,20 +112,79 @@ 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() const;
+  Matrix<rows, columns> &Invert(Matrix<rows, columns> &result) const;
 
   /**
    * @brief Transpose this matrix
    * @param result A buffer to store the result into
    */
-  Matrix<columns, rows> Transpose() const;
+  Matrix<columns, rows> &Transpose(Matrix<columns, rows> &result) const;
 
   /**
-   * @brief Returns the euclidean magnitude of the matrix. Also known as the L2
+   * @brief reduce the matrix so the sum of its elements equal 1
-   * norm
    * @param result a buffer to store the result into
    */
-  float EuclideanNorm() const;
+  Matrix<rows, columns> &Normalize(Matrix<rows, columns> &result) 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
@@ -158,16 +211,8 @@ 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
@@ -176,6 +221,10 @@ 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
@@ -194,68 +243,29 @@ public:
 
   Matrix<rows, columns> operator-(const Matrix<rows, columns> &other) const;
 
-  template <uint8_t other_columns>
+  Matrix<rows, columns> operator*(const Matrix<rows, columns> &other) const;
-  Matrix<rows, other_columns>
-  operator*(const Matrix<columns, other_columns> &other) const;
 
   Matrix<rows, columns> operator*(float scalar) const;
 
-  Matrix<rows, columns> operator/(float scalar) const;
+private:
-
-  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_

README.md

@@ -1,12 +1 @@
-# Introduction
+A Simple matrix math library focused on embedded development which avoids and heap memory allocation unless you explicitly ask for it.
-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

@@ -0,0 +1,81 @@
+#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

@@ -1,20 +0,0 @@
-{
-    "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

@@ -0,0 +1,159 @@
+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

@@ -1,76 +0,0 @@
-# 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/Quaternion.cpp

@@ -1,120 +0,0 @@
-#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

@@ -1,90 +0,0 @@
-#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

@@ -1,540 +0,0 @@
-// 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];
-    }
-  }
-}
-
-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;
-}
-
-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;
-  }
-}
-
-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;
-  }
-}
-
-// ============================================================================
-// 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 ----
-  // Reduce W to upper bidiagonal form using Householder reflections.
-  // For a p×q matrix (p ≤ q after transpose), we do p steps:
-  //   Step k: zero out subdiagonal in column k, then superdiagonal in row k
-
-  float hhVec[5]; // Householder vector storage
-
-  for (uint8_t k = 0; k < p; k++) {
-    // --- Left Householder on column k, rows k..min(m,n)-1 ---
-    {
-      uint8_t len = (m > n) ? m - k : n - k;
-      if (len <= 1)
-        continue;
-
-      // Extract the column segment
-      float x[5];
-      for (uint8_t i = 0; i < len; i++) {
-        x[i] = W[k + i][k];
-      }
-
-      // Compute Householder reflection: H·x = [α, 0, ..., 0]ᵀ
-      float norm = 0;
-      for (uint8_t i = 0; i < len; i++)
-        norm += x[i] * x[i];
-      norm = sqrtf(norm);
-
-      if (norm < 1e-30f)
-        continue;
-
-      float alpha = (x[0] >= 0) ? -norm : norm;
-      float v0 = x[0] - alpha;
-      float vv = v0 * v0 + norm * norm - alpha * x[0];
-      if (vv < 1e-30f)
-        continue;
-
-      float scale = 1.0f / sqrtf(vv);
-      // Store Householder vector (first element is implicit 1, rest in hhVec)
-      hhVec[0] =
-          v0 * scale; // this is the first element of the reflected vector
-      for (uint8_t i = 1; i < len; i++) {
-        hhVec[i] = x[i] * scale;
-      }
-
-      // Apply H from left to W: W = H·W
-      // For each column j, w[k+i][j] -= 2·v_i·(vᵀ·w_col) / (vᵀv)
-      float vvNorm = 1.0f + hhVec[0] * hhVec[0];
-      for (uint8_t i = 1; i < len; i++) {
-        vvNorm += hhVec[i] * hhVec[i];
-      }
-
-      for (uint8_t j = k; j < n; j++) {
-        float dot = 0;
-        for (uint8_t i = 0; i < len; i++) {
-          dot += hhVec[i] * W[k + i][j];
-        }
-        dot *= 2.0f / vvNorm;
-        for (uint8_t i = 0; i < len; i++) {
-          W[k + i][j] -= dot * hhVec[i];
-        }
-      }
-
-      // Apply H from right to QL: QL = QL · H
-      for (uint8_t j = k; j < m; j++) {
-        float dot = 0;
-        for (uint8_t i = 0; i < len; i++) {
-          dot += hhVec[i] * QL[j][k + i];
-        }
-        dot *= 2.0f / vvNorm;
-        for (uint8_t i = 0; i < len; i++) {
-          QL[j][k + i] -= dot * hhVec[i];
-        }
-      }
-    }
-
-    // --- Right Householder on row k, columns k+1..min(m,n)-1 ---
-    {
-      uint8_t len = (p > 1) ? p - 1 - k : 0;
-      if (len <= 0)
-        continue;
-
-      // Extract the row segment
-      float x[5];
-      for (uint8_t i = 0; i < len; i++) {
-        x[i] = W[k][k + 1 + i];
-      }
-
-      // Compute Householder reflection
-      float norm = 0;
-      for (uint8_t i = 0; i < len; i++)
-        norm += x[i] * x[i];
-      norm = sqrtf(norm);
-
-      if (norm < 1e-30f)
-        continue;
-
-      float alpha = (x[0] >= 0) ? -norm : norm;
-      float v0 = x[0] - alpha;
-      float vv = v0 * v0 + norm * norm - alpha * x[0];
-      if (vv < 1e-30f)
-        continue;
-
-      float scale = 1.0f / sqrtf(vv);
-      hhVec[0] = v0 * scale;
-      for (uint8_t i = 1; i < len; i++) {
-        hhVec[i] = x[i] * scale;
-      }
-
-      // Compute vᵀv
-      float vvNorm = 1.0f + hhVec[0] * hhVec[0];
-      for (uint8_t i = 1; i < len; i++) {
-        vvNorm += hhVec[i] * hhVec[i];
-      }
-
-      // Apply H from right to W: W = W·H
-      for (uint8_t i = 0; i < m; i++) {
-        float dot = 0;
-        for (uint8_t j = 0; j < len; j++) {
-          dot += hhVec[j] * W[i][k + 1 + j];
-        }
-        dot *= 2.0f / vvNorm;
-        for (uint8_t j = 0; j < len; j++) {
-          W[i][k + 1 + j] -= dot * hhVec[j];
-        }
-      }
-
-      // Apply H from right to QR: QR = QR · H
-      for (uint8_t i = 0; i < n; i++) {
-        float dot = 0;
-        for (uint8_t j = 0; j < len; j++) {
-          dot += hhVec[j] * QR[i][k + 1 + j];
-        }
-        dot *= 2.0f / vvNorm;
-        for (uint8_t j = 0; j < len; j++) {
-          QR[i][k + 1 + j] -= dot * hhVec[j];
-        }
-      }
-    }
-  }
-
-  // ---- Phase 2: Implicit QR Iteration on Bidiagonal Matrix ----
-  // W now contains the upper bidiagonal matrix B.
-  // We apply implicit QR steps to diagonalize it.
-
-  uint32_t maxIter = 1000;
-  float tol = 1e-8f;
-
-  for (uint32_t iter = 0; iter < maxIter; iter++) {
-    // Deflate: zero out negligible subdiagonal elements
-    for (uint8_t i = p - 1; i > 0; i--) {
-      float test = fabsf(W[i][i - 1]);
-      float scale = fabsf(W[i - 1][i - 1]) + fabsf(W[i][i]);
-      if (test < tol * (scale + 1e-30f)) {
-        W[i][i - 1] = 0;
-      }
-    }
-
-    // Find the smallest unreduced block [start..end]
-    uint8_t start = 0, end = p - 1;
-    for (uint8_t i = 0; i < p - 1; i++) {
-      if (fabsf(W[i + 1][i]) >
-          tol * (fabsf(W[i][i]) + fabsf(W[i + 1][i + 1]) + 1e-30f)) {
-        start = i + 1;
-      }
-    }
-    for (int8_t i = (int8_t)p - 2; i >= 0; i--) {
-      if (fabsf(W[i + 1][i]) >
-          tol * (fabsf(W[i][i]) + fabsf(W[i + 1][i + 1]) + 1e-30f)) {
-        end = (uint8_t)i;
-        break;
-      }
-    }
-
-    // Check convergence of the block
-    if (start >= end) {
-      continue;
-    }
-
-    bool blockConverged = true;
-    for (uint8_t i = start; i <= end; i++) {
-      if (i > start && fabsf(W[i][i - 1]) > tol * (fabsf(W[i - 1][i - 1]) +
-                                                   fabsf(W[i][i]) + 1e-30f)) {
-        blockConverged = false;
-        break;
-      }
-      if (i < end &&
-          fabsf(W[i][i + 1]) >
-              tol * (fabsf(W[i][i]) + fabsf(W[i + 1][i + 1]) + 1e-30f)) {
-        blockConverged = false;
-        break;
-      }
-    }
-    if (blockConverged)
-      continue;
-
-    // Wilkinson shift from bottom 2×2 corner
-    float a = W[end - 1][end - 1];
-    float b = W[end - 1][end];
-    float c = W[end][end - 1];
-    float d = W[end][end];
-
-    float trace = a + d;
-    float det = a * d - b * c;
-    float disc = trace * trace - 4.0f * det;
-    float shift;
-    if (disc >= 0) {
-      float sqrtDisc = sqrtf(disc);
-      float e1 = (trace + sqrtDisc) / 2.0f;
-      float e2 = (trace - sqrtDisc) / 2.0f;
-      shift = (fabsf(e1 - d) < fabsf(e2 - d)) ? e1 : e2;
-    } else {
-      shift = d;
-    }
-
-    // --- Implicit QR step using Givens rotations ---
-    // First, apply Givens rotation from the left to zero out (W[start][start-1]
-    // - shift) For the bidiagonal structure, we process from top to bottom.
-
-    float x = W[start][start] - shift;
-    float y = (start > 0) ? W[start][start - 1] : 0.0f;
-
-    for (uint8_t i = start; i <= end; i++) {
-      float r = sqrtf(x * x + y * y);
-      if (r < 1e-30f) {
-        x = W[i][i];
-        y = (i < end) ? W[i + 1][i] : 0.0f;
-        continue;
-      }
-
-      float cs = x / r;
-      float sn = y / r;
-
-      // Apply Givens from left to rows i, i+1 of W (columns i..p-1)
-      for (uint8_t j = i; j < p; j++) {
-        float t1 = W[i][j];
-        float t2 = W[i + 1][j];
-        W[i][j] = cs * t1 + sn * t2;
-        W[i + 1][j] = -sn * t1 + cs * t2;
-      }
-
-      // Apply Givens from right to columns i, i+1 of W (rows 0..i)
-      if (i > start) {
-        for (uint8_t j = 0; j <= i; j++) {
-          float t1 = W[j][i];
-          float t2 = W[j][i + 1];
-          W[j][i] = cs * t1 + sn * t2;
-          W[j][i + 1] = -sn * t1 + cs * t2;
-        }
-      }
-
-      // Accumulate right transformations into QR
-      for (uint8_t j = 0; j < n; j++) {
-        float t1 = QR[j][i];
-        float t2 = QR[j][i + 1];
-        QR[j][i] = cs * t1 + sn * t2;
-        QR[j][i + 1] = -sn * t1 + cs * t2;
-      }
-
-      // Prepare next Givens rotation
-      x = W[i + 1][i];
-      y = (i + 1 < end) ? W[i + 1][i + 1] : 0.0f;
-    }
-  }
-
-  // ---- Phase 3: Extract Results ----
-  // Singular values are the absolute values of diagonal elements of W
-  for (uint8_t i = 0; i < p; 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]) {
-        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: Compute Final U and Vt ----
-  // If transposeNeeded (wide matrix), we computed SVD(Aᵀ) = Ũ·Σ·Ṽᵀ
-  //   Then U = Ṽ (= QR[:,0:p]) and Vt = Ũᵀ (= QL[:,0:p]ᵀ)
-  // Otherwise, SVD(A) = QL[:,0:p] · Σ · (QR[:,0:p])ᵀ
-  //   So U = QL[:,0:p] and Vt = QR[:,0:p]ᵀ
-
-  for (uint8_t i = 0; i < m; i++) {
-    for (uint8_t j = 0; j < n; j++) {
-      if (j < p) {
-        if (transposeNeeded) {
-          // U = QR[:, 0:p]ᵀ → U[i][j] = QR[j][i]
-          U[i][j] = QR[j][i];
-        } else {
-          // U = QL[:, 0:p]
-          U[i][j] = QL[i][j];
-        }
-      } else {
-        U[i][j] = 0;
-      }
-    }
-  }
-
-  for (uint8_t i = 0; i < n; i++) {
-    for (uint8_t j = 0; j < m; j++) {
-      if (i < p && j < m) {
-        if (transposeNeeded) {
-          // Vt = QL[:, 0:p]ᵀ → Vt[i][j] = QL[j][i]
-          Vt[i][j] = QL[j][i];
-        } else {
-          // Vt = QR[:, 0:p]ᵀ → Vt[i][j] = QR[j][i]
-          Vt[i][j] = QR[j][i];
-        }
-      } else {
-        Vt[i][j] = 0;
-      }
-    }
-  }
-}
-
-#endif

src/SVD.hpp

@@ -1,131 +0,0 @@
-#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 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

@@ -1,162 +0,0 @@
-#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

@@ -1,58 +0,0 @@
-#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,55 +1,21 @@
-# Quaternion tests
+cmake_minimum_required (VERSION 3.11)
-add_executable(quaternion-tests quaternion-tests.cpp)
  
-target_link_libraries(quaternion-tests 
+project ("test_driver")
-    PRIVATE 
+ 
-    quaternion
+include(FetchContent)
-    Catch2::Catch2WithMain
+ 
+FetchContent_Declare(
+    Catch2
+    GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+    GIT_TAG v3.0.1 # or a later release
 )
-
+ 
-# matrix tests
+FetchContent_MakeAvailable(Catch2)
+ 
 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

@@ -0,0 +1,14 @@
+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-timing-tests.cpp

@@ -1,128 +0,0 @@
-// 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

@@ -1,103 +0,0 @@
-// 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

@@ -0,0 +1,7 @@
+# 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-build-blocks-tests.cpp

@@ -1,785 +0,0 @@
-// include the unit test framework first
-#include <catch2/catch_test_macros.hpp>
-#include <catch2/matchers/catch_matchers_floating_point.hpp>
-
-// include the module you're going to test next
-#include "Matrix.hpp"
-#include "SVD.hpp"
-
-// any other libraries
-#include <array>
-#include <cmath>
-#include <iostream>
-
-// ============================================================================
-// Helper: Frobenius norm of a 5×5 matrix
-// ============================================================================
-static float frobeniusNorm5(const Matrix<5, 5> &M) {
-  float sum = 0.0f;
-  for (uint8_t i = 0; i < 5; i++) {
-    for (uint8_t j = 0; j < 5; j++) {
-      float v = M.Get(i, j);
-      sum += v * v;
-    }
-  }
-  return sqrtf(sum);
-}
-
-// ============================================================================
-// Helper: Check if a matrix is orthogonal (Mᵀ·M ≈ I)
-// ============================================================================
-static bool isOrthogonal5(const Matrix<5, 5> &M, float tol = 1e-6f) {
-  Matrix<5, 5> Mt = M.Transpose();
-  Matrix<5, 5> MtM{0};
-  Mt.Mult(M, MtM);
-
-  for (uint8_t i = 0; i < 5; i++) {
-    for (uint8_t j = 0; j < 5; j++) {
-      float expected = (i == j) ? 1.0f : 0.0f;
-      if (fabsf(MtM.Get(i, j) - expected) > tol) {
-        return false;
-      }
-    }
-  }
-  return true;
-}
-
-// ============================================================================
-// TEST 1: ComputeHouseholder
-// ============================================================================
-TEST_CASE("SVD Building Block: ComputeHouseholder", "[Matrix][SVD]") {
-  // Test case: [3, 4] should give alpha = -5 (norm), v normalized
-  // Reference: scipy.linalg.householder([3, 4]) → v ≈ [0.894427191,
-  // 0.447213596], α = -5
-  {
-    float x[] = {3.0f, 4.0f};
-    float v[5] = {0};
-    float alpha = 0;
-
-    float norm = SVD::ComputeHouseholder(x, 2, v, alpha);
-
-    // Norm should be 5.0
-    REQUIRE_THAT(norm, Catch::Matchers::WithinRel(5.0f, 1e-6f));
-
-    // Alpha should be -5 (negative norm)
-    REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-5.0f, 1e-6f));
-
-    // v should be normalized: ||v|| ≈ 1
-    float vNorm = sqrtf(v[0] * v[0] + v[1] * v[1]);
-    REQUIRE_THAT(vNorm, Catch::Matchers::WithinRel(1.0f, 1e-6f));
-
-    // Verify H·x = [alpha, 0]: (I - 2vvᵀ)·x should give [-5, 0]
-    float hx0 = x[0] - 2.0f * v[0] * (v[0] * x[0] + v[1] * x[1]);
-    float hx1 = x[1] - 2.0f * v[1] * (v[0] * x[0] + v[1] * x[1]);
-    REQUIRE_THAT(hx0, Catch::Matchers::WithinRel(alpha, 1e-6f));
-    REQUIRE_THAT(hx1, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [1, 3]
-  // Reference: norm = √10 ≈ 3.16228, alpha = -√10
-  {
-    float x[] = {1.0f, 3.0f};
-    float v[5] = {0};
-    float alpha = 0;
-
-    float norm = SVD::ComputeHouseholder(x, 2, v, alpha);
-
-    REQUIRE_THAT(norm, Catch::Matchers::WithinRel(sqrtf(10.0f), 1e-6f));
-    REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-sqrtf(10.0f), 1e-6f));
-
-    // Verify H·x = [alpha, 0]
-    float dot = v[0] * x[0] + v[1] * x[1];
-    float hx0 = x[0] - 2.0f * v[0] * dot;
-    float hx1 = x[1] - 2.0f * v[1] * dot;
-    REQUIRE_THAT(hx0, Catch::Matchers::WithinRel(alpha, 1e-6f));
-    REQUIRE_THAT(hx1, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [1, 2, 3] (3D)
-  // Reference: norm = √14 ≈ 3.74166
-  {
-    float x[] = {1.0f, 2.0f, 3.0f};
-    float v[5] = {0};
-    float alpha = 0;
-
-    float norm = SVD::ComputeHouseholder(x, 3, v, alpha);
-
-    REQUIRE_THAT(norm, Catch::Matchers::WithinRel(sqrtf(14.0f), 1e-6f));
-    REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-sqrtf(14.0f), 1e-6f));
-
-    // Verify v is normalized
-    float vNorm = sqrtf(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
-    REQUIRE_THAT(vNorm, Catch::Matchers::WithinRel(1.0f, 1e-6f));
-
-    // Verify H·x = [alpha, 0, 0]
-    float dot = v[0] * x[0] + v[1] * x[1] + v[2] * x[2];
-    for (uint8_t i = 0; i < 3; i++) {
-      float hx_i = x[i] - 2.0f * v[i] * dot;
-      if (i == 0) {
-        REQUIRE_THAT(hx_i, Catch::Matchers::WithinRel(alpha, 1e-6f));
-      } else {
-        REQUIRE_THAT(hx_i, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-      }
-    }
-  }
-
-  // Test case: [0, 0, 1] (already has leading zeros)
-  {
-    float x[] = {0.0f, 0.0f, 1.0f};
-    float v[5] = {0};
-    float alpha = 0;
-
-    float norm = SVD::ComputeHouseholder(x, 3, v, alpha);
-
-    REQUIRE_THAT(norm, Catch::Matchers::WithinRel(1.0f, 1e-6f));
-    REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-1.0f, 1e-6f));
-
-    // Verify H·x = [-1, 0, 0]
-    float dot = v[0] * x[0] + v[1] * x[1] + v[2] * x[2];
-    float hx0 = x[0] - 2.0f * v[0] * dot;
-    float hx1 = x[1] - 2.0f * v[1] * dot;
-    float hx2 = x[2] - 2.0f * v[2] * dot;
-    REQUIRE_THAT(hx0, Catch::Matchers::WithinRel(alpha, 1e-6f));
-    REQUIRE_THAT(hx1, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-    REQUIRE_THAT(hx2, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [5, -3, 2, 1] (4D)
-  {
-    float x[] = {5.0f, -3.0f, 2.0f, 1.0f};
-    float v[5] = {0};
-    float alpha = 0;
-
-    float norm = SVD::ComputeHouseholder(x, 4, v, alpha);
-
-    REQUIRE_THAT(norm, Catch::Matchers::WithinRel(sqrtf(39.0f), 1e-6f));
-    REQUIRE_THAT(alpha, Catch::Matchers::WithinRel(-sqrtf(39.0f), 1e-6f));
-
-    // Verify H·x = [alpha, 0, 0, 0]
-    float dot = v[0] * x[0] + v[1] * x[1] + v[2] * x[2] + v[3] * x[3];
-    for (uint8_t i = 0; i < 4; i++) {
-      float hx_i = x[i] - 2.0f * v[i] * dot;
-      if (i == 0) {
-        REQUIRE_THAT(hx_i, Catch::Matchers::WithinRel(alpha, 1e-6f));
-      } else {
-        REQUIRE_THAT(hx_i, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-      }
-    }
-  }
-
-  // Test case: zero vector
-  {
-    float x[] = {0.0f, 0.0f};
-    float v[5] = {0};
-    float alpha = 0;
-
-    float norm = SVD::ComputeHouseholder(x, 2, v, alpha);
-
-    REQUIRE_THAT(norm, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-    REQUIRE(alpha == 0.0f);
-  }
-}
-
-// ============================================================================
-// TEST 2: ApplyHouseholderLeft
-// ============================================================================
-TEST_CASE("SVD Building Block: ApplyHouseholderLeft", "[Matrix][SVD]") {
-  // Test: Apply Householder to zero out column 0, rows 1:2 of a 3×3 matrix
-  // Input: [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
-  // After applying HH on col 0 (rows 1:2): A[2,0] should be ~0
-  {
-    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0,    0, 4.0f, 5.0f, 6.0f, 0,
-                   0,    7.0f, 8.0f, 9.0f, 0, 0,    0,    0,    0,
-                   0,    0,    0,    0,    0, 0,    0};
-
-    // Compute Householder for column 0, rows 1:2 → vector [4, 7]
-    float x[] = {4.0f, 7.0f};
-    float v[5] = {0};
-    float alpha = 0;
-    SVD::ComputeHouseholder(x, 2, v, alpha);
-
-    // Apply from left
-    SVD::ApplyHouseholderLeft(W, v, 1, 2);
-
-    // A[2,0] should be ~0
-    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-
-    // Verify orthogonality of the transformation: W = H·W_original
-    Matrix<5, 5> W_orig{1.0f, 2.0f, 3.0f, 0,    0, 4.0f, 5.0f, 6.0f, 0,
-                        0,    7.0f, 8.0f, 9.0f, 0, 0,    0,    0,    0,
-                        0,    0,    0,    0,    0, 0,    0};
-
-    // Compute H_left explicitly: I - 2*v*vᵀ (on rows 1:2)
-    Matrix<5, 5> H_left{0};
-    for (uint8_t i = 0; i < 5; i++) {
-      H_left[i][i] = 1.0f;
-    }
-    // Apply -2*v*vᵀ to the sub-block
-    float vv = v[0] * v[0] + v[1] * v[1];
-    for (uint8_t i = 1; i <= 2; i++) {
-      for (uint8_t j = 1; j <= 2; j++) {
-        H_left[i][j] -= 2.0f * v[i - 1] * v[j - 1] / vv;
-      }
-    }
-
-    // Verify: W ≈ H_left · W_orig
-    Matrix<5, 5> HLeftW{0};
-    H_left.Mult(W_orig, HLeftW);
-    float err = frobeniusNorm5(W - HLeftW);
-    REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
-
-    // Verify H_left is orthogonal
-    REQUIRE(isOrthogonal5(H_left));
-  }
-
-  // Test: Apply to a larger block (4 rows)
-  {
-    Matrix<5, 5> W{1.0f, 2.0f, 0,    0, 0, 3.0f, 4.0f, 0,    0,
-                   0,    5.0f, 6.0f, 0, 0, 0,    7.0f, 8.0f, 0,
-                   0,    0,    0,    0, 0, 0,    0};
-
-    // Householder on [3, 5, 7] (rows 1:3)
-    float x[] = {3.0f, 5.0f, 7.0f};
-    float v[5] = {0};
-    float alpha = 0;
-    SVD::ComputeHouseholder(x, 3, v, alpha);
-
-    SVD::ApplyHouseholderLeft(W, v, 1, 3);
-
-    // A[2,0] and A[3,0] should be ~0
-    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-    REQUIRE_THAT(W.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-}
-
-// ============================================================================
-// TEST 3: ApplyHouseholderRight
-// ============================================================================
-TEST_CASE("SVD Building Block: ApplyHouseholderRight", "[Matrix][SVD]") {
-  // Test: Apply Householder to zero out row 0, cols 1:2 of a 3×3 matrix
-  // Input: [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
-  // After applying HH on row 0 (cols 1:2): A[0,2] should be ~0
-  {
-    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0,    0, 4.0f, 5.0f, 6.0f, 0,
-                   0,    7.0f, 8.0f, 9.0f, 0, 0,    0,    0,    0,
-                   0,    0,    0,    0,    0, 0,    0};
-
-    // Householder for row 0, cols 1:2 → vector [2, 3]
-    float x[] = {2.0f, 3.0f};
-    float v[5] = {0};
-    float alpha = 0;
-    SVD::ComputeHouseholder(x, 2, v, alpha);
-
-    // Apply from right
-    SVD::ApplyHouseholderRight(W, v, 1, 2);
-
-    // A[0,2] should be ~0
-    REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-
-    // Verify W ≈ W_orig · H_right
-    Matrix<5, 5> W_orig{1.0f, 2.0f, 3.0f, 0,    0, 4.0f, 5.0f, 6.0f, 0,
-                        0,    7.0f, 8.0f, 9.0f, 0, 0,    0,    0,    0,
-                        0,    0,    0,    0,    0, 0,    0};
-
-    // Compute H_right = I - 2*v*vᵀ (on cols 1:2)
-    Matrix<5, 5> H_right{0};
-    for (uint8_t i = 0; i < 5; i++) {
-      H_right[i][i] = 1.0f;
-    }
-    float vv = v[0] * v[0] + v[1] * v[1];
-    for (uint8_t i = 1; i <= 2; i++) {
-      for (uint8_t j = 1; j <= 2; j++) {
-        H_right[i][j] -= 2.0f * v[i - 1] * v[j - 1] / vv;
-      }
-    }
-
-    Matrix<5, 5> WOrigH{0};
-    W_orig.Mult(H_right, WOrigH);
-    float err = frobeniusNorm5(W - WOrigH);
-    REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-4f));
-
-    // Verify H_right is orthogonal
-    REQUIRE(isOrthogonal5(H_right));
-  }
-
-  // Test: Apply to wider block (4 cols)
-  {
-    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 4.0f, 0, 5.0f, 6.0f, 7.0f, 8.0f,
-                   0,    0,    0,    0,    0, 0,    0,    0,    0,
-                   0,    0,    0,    0,    0, 0,    0};
-
-    // Householder on [2, 3, 4] (cols 1:3)
-    float x[] = {2.0f, 3.0f, 4.0f};
-    float v[5] = {0};
-    float alpha = 0;
-    SVD::ComputeHouseholder(x, 3, v, alpha);
-
-    SVD::ApplyHouseholderRight(W, v, 1, 3);
-
-    // A[0,2] and A[0,3] should be ~0
-    REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-    REQUIRE_THAT(W.Get(0, 3), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-}
-
-// ============================================================================
-// TEST 4: ComputeGivens
-// ============================================================================
-TEST_CASE("SVD Building Block: ComputeGivens", "[Matrix][SVD]") {
-  // Test case: [3, 4] → c = 0.6, s = 0.8 (3-4-5 triangle)
-  {
-    float c, s;
-    SVD::ComputeGivens(3.0f, 4.0f, c, s);
-
-    REQUIRE_THAT(c, Catch::Matchers::WithinRel(0.6f, 1e-6f));
-    REQUIRE_THAT(s, Catch::Matchers::WithinRel(0.8f, 1e-6f));
-
-    // Verify: [c s; -s c] · [3; 4] = [5; 0]
-    float r = c * 3.0f + s * 4.0f;
-    float z = -s * 3.0f + c * 4.0f;
-    REQUIRE_THAT(r, Catch::Matchers::WithinRel(5.0f, 1e-6f));
-    REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-
-    // Verify c² + s² = 1
-    REQUIRE_THAT(c * c + s * s, Catch::Matchers::WithinRel(1.0f, 1e-6f));
-  }
-
-  // Test case: [1, 0] → c = 1, s = 0
-  {
-    float c, s;
-    SVD::ComputeGivens(1.0f, 0.0f, c, s);
-    REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-6f));
-    REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [0, 5] → c = 0, s = 1
-  {
-    float c, s;
-    SVD::ComputeGivens(0.0f, 5.0f, c, s);
-    REQUIRE_THAT(c, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-    REQUIRE_THAT(s, Catch::Matchers::WithinRel(1.0f, 1e-6f));
-
-    // Verify: [c s; -s c] · [0; 5] = [5; 0]
-    float r = c * 0.0f + s * 5.0f;
-    float z = -s * 0.0f + c * 5.0f;
-    REQUIRE_THAT(r, Catch::Matchers::WithinRel(5.0f, 1e-6f));
-    REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [-3, -4] → c = -0.6, s = -0.8
-  {
-    float c, s;
-    SVD::ComputeGivens(-3.0f, -4.0f, c, s);
-    REQUIRE_THAT(c, Catch::Matchers::WithinRel(-0.6f, 1e-6f));
-    REQUIRE_THAT(s, Catch::Matchers::WithinRel(-0.8f, 1e-6f));
-
-    // Verify: [c s; -s c] · [-3; -4] = [5; 0]
-    float r = c * (-3.0f) + s * (-4.0f);
-    float z = -s * (-3.0f) + c * (-4.0f);
-    REQUIRE_THAT(r, Catch::Matchers::WithinRel(5.0f, 1e-6f));
-    REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [1, -1] → c = 1/√2, s = -1/√2 (45°)
-  {
-    float c, s;
-    SVD::ComputeGivens(1.0f, -1.0f, c, s);
-    float invSqrt2 = 1.0f / sqrtf(2.0f);
-    REQUIRE_THAT(c, Catch::Matchers::WithinRel(invSqrt2, 1e-6f));
-    REQUIRE_THAT(s, Catch::Matchers::WithinRel(-invSqrt2, 1e-6f));
-
-    // Verify: [c s; -s c] · [1; -1] = [√2; 0]
-    float r = c * 1.0f + s * (-1.0f);
-    float z = -s * 1.0f + c * (-1.0f);
-    REQUIRE_THAT(r, Catch::Matchers::WithinRel(sqrtf(2.0f), 1e-6f));
-    REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [0, 0] → c = 1, s = 0 (identity)
-  {
-    float c, s;
-    SVD::ComputeGivens(0.0f, 0.0f, c, s);
-    REQUIRE_THAT(c, Catch::Matchers::WithinRel(1.0f, 1e-6f));
-    REQUIRE_THAT(s, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-
-  // Test case: [7, 24] → c = 7/25, s = 24/25 (7-24-25 triangle)
-  {
-    float c, s;
-    SVD::ComputeGivens(7.0f, 24.0f, c, s);
-    REQUIRE_THAT(c, Catch::Matchers::WithinRel(7.0f / 25.0f, 1e-6f));
-    REQUIRE_THAT(s, Catch::Matchers::WithinRel(24.0f / 25.0f, 1e-6f));
-
-    float r = c * 7.0f + s * 24.0f;
-    float z = -s * 7.0f + c * 24.0f;
-    REQUIRE_THAT(r, Catch::Matchers::WithinRel(25.0f, 1e-6f));
-    REQUIRE_THAT(z, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-}
-
-// ============================================================================
-// TEST 5: ApplyGivensLeft
-// ============================================================================
-TEST_CASE("SVD Building Block: ApplyGivensLeft", "[Matrix][SVD]") {
-  // Test: Apply Givens to zero out W[1,0] of a 2×2 matrix
-  // Input: [[3, 4], [1, 2]]
-  // Givens on rows 0,1 with x=W[0,0]=3, y=W[1,0]=1
-  {
-    Matrix<5, 5> W{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0,
-                   0,    0,    0, 0, 0, 0,    0,    0, 0, 0, 0, 0};
-
-    float c, s;
-    SVD::ComputeGivens(3.0f, 1.0f, c, s);
-
-    SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4);
-
-    // W[1,0] should be ~0
-    REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-
-    // Verify W ≈ G · W_orig
-    Matrix<5, 5> W_orig{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0,
-                        0,    0,    0, 0, 0, 0,    0,    0, 0, 0, 0, 0};
-
-    // Givens rotation matrix (5×5)
-    Matrix<5, 5> G{0};
-    for (uint8_t i = 0; i < 5; i++) {
-      G[i][i] = 1.0f;
-    }
-    G[0][0] = c;
-    G[0][1] = s;
-    G[1][0] = -s;
-    G[1][1] = c;
-
-    Matrix<5, 5> GW{0};
-    G.Mult(W_orig, GW);
-    float err = frobeniusNorm5(W - GW);
-    REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-
-    // Verify G is orthogonal
-    REQUIRE(isOrthogonal5(G));
-  }
-
-  // Test: Apply to larger range of columns
-  {
-    Matrix<5, 5> W{3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 1.0f, 2.0f, 3.0f, 4.0f,
-                   5.0f, 0,    0,    0,    0,    0,    0,    0,    0,
-                   0,    0,    0,    0,    0,    0,    0};
-
-    float c, s;
-    SVD::ComputeGivens(3.0f, 1.0f, c, s);
-
-    SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4);
-
-    REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-}
-
-// ============================================================================
-// TEST 6: ApplyGivensRight
-// ============================================================================
-TEST_CASE("SVD Building Block: ApplyGivensRight", "[Matrix][SVD]") {
-  // Test: Apply Givens to zero out W[0,1] of a 2×2 matrix
-  // Input: [[3, 4], [1, 2]]
-  // Givens on cols 0,1 with x=W[0,0]=3, y=W[0,1]=4
-  {
-    Matrix<5, 5> W{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0,
-                   0,    0,    0, 0, 0, 0,    0,    0, 0, 0, 0, 0};
-
-    float c, s;
-    SVD::ComputeGivens(3.0f, 4.0f, c, s);
-
-    SVD::ApplyGivensRight(W, 0, 1, c, s, 0, 4);
-
-    // W[0,1] should be ~0
-    REQUIRE_THAT(W.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-
-    // Verify W ≈ W_orig · G
-    Matrix<5, 5> W_orig{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 0, 0, 0,
-                        0,    0,    0, 0, 0, 0,    0,    0, 0, 0, 0, 0};
-
-    // Givens rotation matrix (5×5)
-    Matrix<5, 5> G{0};
-    for (uint8_t i = 0; i < 5; i++) {
-      G[i][i] = 1.0f;
-    }
-    G[0][0] = c;
-    G[0][1] = -s;
-    G[1][0] = s;
-    G[1][1] = c;
-
-    Matrix<5, 5> WG{0};
-    W_orig.Mult(G, WG);
-    float err = frobeniusNorm5(W - WG);
-    REQUIRE_THAT(err, Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-
-    // Verify G is orthogonal
-    REQUIRE(isOrthogonal5(G));
-  }
-
-  // Test: Apply to larger range of rows
-  {
-    Matrix<5, 5> W{3.0f, 4.0f, 0, 0, 0, 1.0f, 2.0f, 0, 0, 0, 5.0f, 6.0f, 0,
-                   0,    0,    0, 0, 0, 0,    0,    0, 0, 0, 0,    0};
-
-    float c, s;
-    SVD::ComputeGivens(3.0f, 4.0f, c, s);
-
-    SVD::ApplyGivensRight(W, 0, 1, c, s, 0, 2);
-
-    REQUIRE_THAT(W.Get(0, 1), Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-}
-
-// ============================================================================
-// TEST 7: Full Bidiagonalization (composing Householder steps)
-// ============================================================================
-TEST_CASE("SVD Building Block: Householder Bidiagonalization",
-          "[Matrix][SVD]") {
-  // Test: Bidiagonalize a 3×3 matrix and verify reconstruction
-  // Input: [[1, 2, 3], [4, 5, 6], [7, 8, 10]]
-  {
-    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0,     0, 4.0f, 5.0f, 6.0f, 0,
-                   0,    7.0f, 8.0f, 10.0f, 0, 0,    0,    0,    0,
-                   0,    0,    0,    0,     0, 0,    0};
-
-    // Step 1: Left HH on column 0, rows 1:2 → zero out W[2,0]
-    {
-      float x[] = {4.0f, 7.0f};
-      float v[5] = {0};
-      float alpha = 0;
-      SVD::ComputeHouseholder(x, 2, v, alpha);
-      SVD::ApplyHouseholderLeft(W, v, 1, 2);
-    }
-
-    // Step 2: Right HH on row 0, cols 1:2 → zero out W[0,2]
-    {
-      float x[] = {W.Get(0, 1), W.Get(0, 2)};
-      float v[5] = {0};
-      float alpha = 0;
-      SVD::ComputeHouseholder(x, 2, v, alpha);
-      SVD::ApplyHouseholderRight(W, v, 1, 2);
-    }
-
-    // Step 3: Left HH on column 1, rows 2:2 → nothing to do (single element)
-
-    // Verify bidiagonal structure: for 3x3, zero elements are A[2][0] (below
-    // subdiag in col 0) and A[0][2] (above superdiag in row 0) A[2][1] is the
-    // subdiagonal element of col 1 — valid in bidiagonal form
-    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
-    REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-4f));
-  }
-
-  // Test: Bidiagonalize a 4×3 matrix
-  {
-    Matrix<5, 5> W{1.0f, 2.0f, 3.0f, 0,    0, 4.0f, 5.0f,  6.0f,  0,
-                   0,    7.0f, 8.0f, 9.0f, 0, 0,    10.0f, 11.0f, 12.0f,
-                   0,    0,    0,    0,    0, 0,    0};
-
-    // Step 1: Left HH on col 0, rows 1:3 → zero out W[2,0], W[3,0]
-    {
-      float x[] = {4.0f, 7.0f, 10.0f};
-      float v[5] = {0};
-      float alpha = 0;
-      SVD::ComputeHouseholder(x, 3, v, alpha);
-      SVD::ApplyHouseholderLeft(W, v, 1, 3);
-    }
-
-    // Step 2: Right HH on row 0, cols 1:2 → zero out W[0,2]
-    {
-      float x[] = {W.Get(0, 1), W.Get(0, 2)};
-      float v[5] = {0};
-      float alpha = 0;
-      SVD::ComputeHouseholder(x, 2, v, alpha);
-      SVD::ApplyHouseholderRight(W, v, 1, 2);
-    }
-
-    // Step 3: Left HH on col 1, rows 2:3 → zero out W[3,1]
-    {
-      float x[] = {W.Get(2, 1), W.Get(3, 1)};
-      float v[5] = {0};
-      float alpha = 0;
-      SVD::ComputeHouseholder(x, 2, v, alpha);
-      SVD::ApplyHouseholderLeft(W, v, 2, 3);
-    }
-
-    // Verify bidiagonal structure
-    REQUIRE_THAT(W.Get(2, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-    REQUIRE_THAT(W.Get(3, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-    REQUIRE_THAT(W.Get(3, 1), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-  }
-
-  // Test: Diagonal matrix (no transformations needed)
-  {
-    Matrix<5, 5> W{10.0f, 0, 0, 0, 0, 0, 5.0f, 0, 0, 0, 0, 0, 2.0f,
-                   0,     0, 0, 0, 0, 0, 0,    0, 0, 0, 0, 0};
-
-    // Householder on zero vector should be identity
-    float x[] = {0.0f, 0.0f};
-    float v[5] = {0};
-    float alpha = 0;
-    SVD::ComputeHouseholder(x, 2, v, alpha);
-
-    // Applying identity should not change anything
-    Matrix<5, 5> W_copy{10.0f, 0, 0, 0, 0, 0, 5.0f, 0, 0, 0, 0, 0, 2.0f,
-                        0,     0, 0, 0, 0, 0, 0,    0, 0, 0, 0, 0};
-    SVD::ApplyHouseholderLeft(W_copy, v, 1, 2);
-
-    REQUIRE_THAT(frobeniusNorm5(W - W_copy),
-                 Catch::Matchers::WithinAbs(0.0f, 1e-6f));
-  }
-}
-
-// ============================================================================
-// //
-// ============================================================================
-// TEST 8: Givens QR step on bidiagonal matrix
-// ===========================================================================
-TEST_CASE("SVD Building Block: Givens QR Step on Bidiagonal", "[Matrix][SVD]") {
-  // Test: Apply left Givens to zero subdiagonal of a bidiagonal matrix,
-  // then apply right Givens with restricted row range to restore bidiagonal
-  // form.
-  //
-  // Input: 3x3 bidiagonal [[1, 2, 0], [3, -4, 5], [0, 6, -7]]
-  // Step 1: Left Givens on rows 0,1 with x=W[0][0]=1, y=W[1][0]=3 -> zero
-  // W[1][0] Step 2: Right Givens on cols 1,2 with x=W[0][1], y=W[0][2] -> zero
-  // W[0][2]
-  //         Only applied to row 0 (to not reintroduce subdiagonal non-zeros)
-  {
-    Matrix<5, 5> W{1.0f, 2.0f, 0.0f, 0,     0, 3.0f, -4.0f, 5.0f, 0,
-                   0,    0.0f, 6.0f, -7.0f, 0, 0,    0,     0,    0,
-                   0,    0,    0,    0,     0, 0,    0};
-
-    float c, s;
-    SVD::ComputeGivens(W.Get(0, 0), W.Get(1, 0), c, s);
-
-    // Apply from left to zero subdiagonal at W[1][0]
-    SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4);
-
-    REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-
-    // After left Givens, W[0][2] may have become non-zero (fill-in from row 0)
-    // Apply right Givens to cols 1,2 with x=W[0][1], y=W[0][2] -> zero W[0][2]
-    // Only apply to rows 0 (to preserve bidiagonal structure below row 0)
-    float c2, s2;
-    SVD::ComputeGivens(W.Get(0, 1), W.Get(0, 2), c2, s2);
-    SVD::ApplyGivensRight(W, 1, 2, c2, s2, 0, 0);
-
-    // Should be bidiagonal: W[1][0] ~ 0 (from left Givens), W[0][2] ~ 0 (from
-    // right Givens)
-    REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-    REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-  }
-
-  // Test: Verify that a full QR step (left + right Givens) preserves the
-  // bidiagonal structure when applied correctly with proper row ranges.
-  {
-    Matrix<5, 5> W{2.0f, 3.0f, 0, 0, 0, -1.0f, 4.0f, 5.0f, 0, 0, 0, 6.0f, -7.0f,
-                   0,    0,    0, 0, 0, 0,     0,    0,    0, 0, 0, 0};
-
-    // Left Givens on col 0 (rows 0,1)
-    float c, s;
-    SVD::ComputeGivens(W.Get(0, 0), W.Get(1, 0), c, s);
-    SVD::ApplyGivensLeft(W, 0, 1, c, s, 0, 4);
-
-    REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-
-    // Right Givens on row 0 (cols 1,2) - only affect row 0
-    float c2, s2;
-    SVD::ComputeGivens(W.Get(0, 1), W.Get(0, 2), c2, s2);
-    SVD::ApplyGivensRight(W, 1, 2, c2, s2, 0, 0);
-
-    // Bidiagonal structure preserved
-    REQUIRE_THAT(W.Get(1, 0), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-    REQUIRE_THAT(W.Get(0, 2), Catch::Matchers::WithinAbs(0.0f, 1e-5f));
-  }
-}
-
-// ============================================================================TEST
-// 9: Orthogonality preservation of Householder transformations
-// ============================================================================
-TEST_CASE("SVD Building Block: Householder preserves orthogonality",
-          "[Matrix][SVD]") {
-  // Starting with an orthogonal matrix, applying Householder should preserve it
-  {
-    // Identity matrix is orthogonal
-    Matrix<5, 5> M{0};
-    for (uint8_t i = 0; i < 5; i++) {
-      M[i][i] = 1.0f;
-    }
-
-    // Householder on first 3 elements of column 0
-    float x[] = {1.0f, 0.0f, 0.0f};
-    float v[5] = {0};
-    float alpha = 0;
-    SVD::ComputeHouseholder(x, 3, v, alpha);
-
-    // Apply from left
-    Matrix<5, 5> M_left = M;
-    SVD::ApplyHouseholderLeft(M_left, v, 0, 2);
-
-    // M_left should still be orthogonal
-    REQUIRE(isOrthogonal5(M_left));
-
-    // Apply from right
-    Matrix<5, 5> M_right = M;
-    SVD::ApplyHouseholderRight(M_right, v, 0, 2);
-
-    REQUIRE(isOrthogonal5(M_right));
-  }
-
-  // Random orthogonal matrix (rotation)
-  {
-    float c = sqrtf(0.5f);
-    float s = sqrtf(0.5f);
-    Matrix<5, 5> M{0};
-    M[0][0] = c;
-    M[0][1] = -s;
-    M[1][0] = s;
-    M[1][1] = c;
-    for (uint8_t i = 2; i < 5; i++) {
-      M[i][i] = 1.0f;
-    }
-
-    REQUIRE(isOrthogonal5(M));
-
-    // Apply Householder on rows 0,1
-    float x[] = {c, s};
-    float v[5] = {0};
-    float alpha = 0;
-    SVD::ComputeHouseholder(x, 2, v, alpha);
-
-    Matrix<5, 5> M_test = M;
-    SVD::ApplyHouseholderLeft(M_test, v, 0, 1);
-
-    REQUIRE(isOrthogonal5(M_test));
-  }
-}
-
-// ============================================================================
-// TEST 10: Orthogonality preservation of Givens transformations
-// ============================================================================
-TEST_CASE("SVD Building Block: Givens preserves orthogonality",
-          "[Matrix][SVD]") {
-  // Starting with an orthogonal matrix, applying Givens should preserve it
-  {
-    Matrix<5, 5> M{0};
-    for (uint8_t i = 0; i < 5; i++) {
-      M[i][i] = 1.0f;
-    }
-
-    float c, s;
-    SVD::ComputeGivens(3.0f, 4.0f, c, s);
-
-    // Apply from left
-    Matrix<5, 5> M_left = M;
-    SVD::ApplyGivensLeft(M_left, 0, 1, c, s, 0, 4);
-
-    REQUIRE(isOrthogonal5(M_left));
-
-    // Apply from right
-    Matrix<5, 5> M_right = M;
-    SVD::ApplyGivensRight(M_right, 0, 1, c, s, 0, 4);
-
-    REQUIRE(isOrthogonal5(M_right));
-  }
-}

unit-tests/svd-integration-test.cpp

@@ -1,269 +0,0 @@
-#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::WithinRel(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::WithinRel(0.0f, 1e-2f));
-  REQUIRE_THAT(vtErr, Catch::Matchers::WithinRel(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::WithinRel(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::WithinRel(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: U (3x5) * diag(sigma) (5x3) = 3x3, then * Vt (3x5) =
-  // 3x5
-  Matrix<3, 5> recon{0};
-  Matrix<3, 5> Usig{0};
-  for (int i = 0; i < 3; i++)
-    for (int j = 0; j < 5; j++)
-      Usig[i][j] = U.Get(i, j) * sigma.Get(j, 0);
-
-  // For wide matrix: A = U * Sigma * Vt where U is 3x5, Sigma is 5x5
-  // (diagonal), Vt is 5x5 But our implementation returns sigma as 3x1 and Vt as
-  // 3x5 So we need: recon = Usig (3x5) * Vt (3x5)^T ... no that doesn't work
-  // either The SVD for wide matrices is: A = U * Sigma * Vt where:
-  //   U is m×m (3×3), Sigma is m×n (3×5), Vt is n×n (5×5)
-  // But our API returns U as m×n (3×5), sigma as n×1 (3×1), Vt as n×n (3×5)
-  // So: recon = U (3x5) * diag(sigma) (5x5) * Vt (5x5)^T ...
-  // Actually, looking at the implementation, for wide matrices we swap roles.
-  // Let me just check reconstruction using the actual dimensions returned.
-
-  // For wide matrix: A (3x5) = U (3x5) * diag(sigma) (5x5 padded) * Vt (5x5)
-  // But our API returns Vt as 3x5, not 5x5
-  // The implementation stores: U = QR[:,0:p]^T (3x5), sigma (3x1), Vt =
-  // QL[:,0:p]^T (3x5) Reconstruction: A[i][j] = sum_k U[i][k]*sigma[k]*Vt[j][k]
-
-  float err2 = 0.0f;
-  for (int i = 0; i < 3; i++) {
-    for (int j = 0; j < 5; j++) {
-      float recon_val = 0.0f;
-      for (int k = 0; k < 3; k++) {
-        recon_val += U.Get(i, k) * sigma.Get(k, 0) * Vt.Get(j, k);
-      }
-      float diff = recon_val - A.Get(i, j);
-      err2 += diff * diff;
-    }
-  }
-  err2 = sqrtf(err2);
-  REQUIRE_THAT(err2, Catch::Matchers::WithinRel(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::WithinRel(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::WithinRel(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

@@ -1,482 +0,0 @@
-#!/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

@@ -1,36 +0,0 @@
-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

@@ -1,45 +0,0 @@
-// 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);
-    }
-}