Made my own equally wrong QR factorization

.vscode/settings.json

@@ -76,8 +76,5 @@
     "clangd.enable": true,
     "C_Cpp.dimInactiveRegions": false,
     "editor.defaultFormatter": "xaver.clang-format",
-    "clangd.inactiveRegions.useBackgroundHighlight": true,
-    "clangd.arguments": [
-        "--compile-commands-dir=${workspaceFolder}/build"
-    ],
+    "clangd.inactiveRegions.useBackgroundHighlight": true
 }

src/Matrix.cpp

@@ -533,12 +533,7 @@ void Matrix<rows, columns>::QRDecomposition(Matrix<rows, columns> &Q,
       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);
-    }
+    u = u / u.EuclideanNorm();
     Q.SetSubMatrix(0, column, u);
 
     // -----------------------

unit-tests/matrix-tests.cpp

@@ -336,27 +336,27 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
     Matrix<3, 3> mat4 = startMatrix;
 
     Matrix<2, 2> mat5{10, 11, 12, 13};
-    mat4.SetSubMatrix(0, 0, mat5);
+    mat4.SetSubMatrix<2, 2, 0, 0>(mat5);
     REQUIRE(mat4.Get(0, 0) == 10);
     REQUIRE(mat4.Get(0, 1) == 11);
     REQUIRE(mat4.Get(1, 0) == 12);
     REQUIRE(mat4.Get(1, 1) == 13);
 
     mat4 = startMatrix;
-    mat4.SetSubMatrix(1, 1, mat5);
+    mat4.SetSubMatrix<2, 2, 1, 1>(mat5);
     REQUIRE(mat4.Get(1, 1) == 10);
     REQUIRE(mat4.Get(1, 2) == 11);
     REQUIRE(mat4.Get(2, 1) == 12);
     REQUIRE(mat4.Get(2, 2) == 13);
 
     Matrix<3, 1> mat6{10, 11, 12};
-    mat4.SetSubMatrix(0, 0, mat6);
+    mat4.SetSubMatrix<3, 1, 0, 0>(mat6);
     REQUIRE(mat4.Get(0, 0) == 10);
     REQUIRE(mat4.Get(1, 0) == 11);
     REQUIRE(mat4.Get(2, 0) == 12);
 
     Matrix<1, 3> mat7{10, 11, 12};
-    mat4.SetSubMatrix(0, 0, mat7);
+    mat4.SetSubMatrix<1, 3, 0, 0>(mat7);
     REQUIRE(mat4.Get(0, 0) == 10);
     REQUIRE(mat4.Get(0, 1) == 11);
     REQUIRE(mat4.Get(0, 2) == 12);
@@ -375,7 +375,7 @@ float matrixSum(const Matrix<rows, columns> &matrix) {
 
 // TODO: Add test for scalar division
 
-TEST_CASE("Euclidean Norm", "Matrix") {
+TEST_CASE("Normalization", "Matrix") {
 
   SECTION("2x2 Normalize") {
     Matrix<2, 2> mat1{1, 2, 3, 4};
@@ -423,48 +423,48 @@ TEST_CASE("Euclidean Norm", "Matrix") {
 }
 
 TEST_CASE("QR Decompositions", "Matrix") {
-  // SECTION("2x2 QRDecomposition") {
-  //   Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
-  //   Matrix<2, 2> Q{}, R{};
-  //   A.QRDecomposition(Q, R);
+  SECTION("2x2 QRDecomposition") {
+    Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
+    Matrix<2, 2> Q{}, R{};
+    A.QRDecomposition(Q, R);
 
-  //   // Check that Q * R ≈ A
-  //   Matrix<2, 2> QR{};
-  //   Q.Mult(R, QR);
-  //   for (int i = 0; i < 2; ++i) {
-  //     for (int j = 0; j < 2; ++j) {
-  //       REQUIRE_THAT(QR[i][j], Catch::Matchers::WithinRel(A[i][j], 1e-4f));
-  //     }
-  //   }
+    // Check that Q * R ≈ A
+    Matrix<2, 2> QR{};
+    Q.Mult(R, QR);
+    for (int i = 0; i < 2; ++i) {
+      for (int j = 0; j < 2; ++j) {
+        REQUIRE_THAT(QR[i][j], Catch::Matchers::WithinRel(A[i][j], 1e-4f));
+      }
+    }
 
-  //   // Check that Q is orthonormal: Qᵀ * Q ≈ I
-  //   Matrix<2, 2> Qt = Q.Transpose();
-  //   Matrix<2, 2> QtQ{};
-  //   Qt.Mult(Q, QtQ);
-  //   for (int i = 0; i < 2; ++i) {
-  //     for (int j = 0; j < 2; ++j) {
-  //       if (i == j)
-  //         REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f));
-  //       else
-  //         REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
-  //     }
-  //   }
+    // Check that Q is orthonormal: Qᵀ * Q ≈ I
+    Matrix<2, 2> Qt = Q.Transpose();
+    Matrix<2, 2> QtQ{};
+    Qt.Mult(Q, QtQ);
+    for (int i = 0; i < 2; ++i) {
+      for (int j = 0; j < 2; ++j) {
+        if (i == j)
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f));
+        else
+          REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
+      }
+    }
 
-  //   // Optional: R should be upper triangular
-  //   REQUIRE(std::fabs(R[1][0]) < 1e-4f);
+    // Optional: R should be upper triangular
+    REQUIRE(std::fabs(R[1][0]) < 1e-4f);
 
-  //   // check that all Q values are correct
-  //   REQUIRE_THAT(Q[0][0], Catch::Matchers::WithinRel(0.3162f, 1e-4f));
-  //   REQUIRE_THAT(Q[0][1], Catch::Matchers::WithinRel(0.94868f, 1e-4f));
-  //   REQUIRE_THAT(Q[1][0], Catch::Matchers::WithinRel(0.94868f, 1e-4f));
-  //   REQUIRE_THAT(Q[1][1], Catch::Matchers::WithinRel(-0.3162f, 1e-4f));
+    // check that all Q values are correct
+    REQUIRE_THAT(Q[0][0], Catch::Matchers::WithinRel(0.3162f, 1e-4f));
+    REQUIRE_THAT(Q[0][1], Catch::Matchers::WithinRel(0.94868f, 1e-4f));
+    REQUIRE_THAT(Q[1][0], Catch::Matchers::WithinRel(0.94868f, 1e-4f));
+    REQUIRE_THAT(Q[1][1], Catch::Matchers::WithinRel(-0.3162f, 1e-4f));
 
-  //   // check that all R values are correct
-  //   REQUIRE_THAT(R[0][0], Catch::Matchers::WithinRel(3.16228f, 1e-4f));
-  //   REQUIRE_THAT(R[0][1], Catch::Matchers::WithinRel(4.42719f, 1e-4f));
-  //   REQUIRE_THAT(R[1][0], Catch::Matchers::WithinRel(0.0f, 1e-4f));
-  //   REQUIRE_THAT(R[1][1], Catch::Matchers::WithinRel(0.63246f, 1e-4f));
-  // }
+    // check that all R values are correct
+    REQUIRE_THAT(R[0][0], Catch::Matchers::WithinRel(3.16228f, 1e-4f));
+    REQUIRE_THAT(R[0][1], Catch::Matchers::WithinRel(4.42719f, 1e-4f));
+    REQUIRE_THAT(R[1][0], Catch::Matchers::WithinRel(0.0f, 1e-4f));
+    REQUIRE_THAT(R[1][1], Catch::Matchers::WithinRel(0.63246f, 1e-4f));
+  }
 
   SECTION("3x3 QRDecomposition") {
     // this symmetrix tridiagonal matrix is well behaved for testing