Compare commits

1 Commits

Author SHA1 Message Date
c0e07eac56 Made my own equally wrong QR factorization
Some checks failed
Merge-Checker / build_and_test (pull_request) Failing after 17s
2025-06-02 20:39:44 -04:00
8 changed files with 277 additions and 311 deletions

View File

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

View File

@@ -2,11 +2,8 @@
This matrix math library is focused on embedded development and avoids any 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. It uses templates to pre-allocate matrices on the stack.
# Building There are still several operations that are works in progress such as:
1. Initialize the repositiory with the command: - Add a function to calculate eigenvalues/vectors
```bash - Add a function to compute RREF
cmake -S . -B build -G Ninja - Add a function for SVD decomposition
``` - Add a function for LQ decomposition
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`

View File

@@ -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 #ifdef MATRIX_H_ // since the .cpp file has to be included by the .hpp file this
// will evaluate to true // will evaluate to true
#include "Matrix.hpp" #include "Matrix.hpp"
@@ -13,6 +6,12 @@
#include <cmath> #include <cmath>
#include <cstdlib> #include <cstdlib>
#include <cstring> #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> template <uint8_t rows, uint8_t columns>
Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) { Matrix<rows, columns>::Matrix(const std::array<float, rows * columns> &array) {
@@ -26,14 +25,6 @@ Matrix<rows, columns>::Matrix(Args... args) {
static_cast<uint16_t>(columns)}; static_cast<uint16_t>(columns)};
std::initializer_list<float> initList{static_cast<float>(args)...}; 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 // choose whichever buffer size is smaller for the copy length
uint32_t minSize = uint32_t minSize =
std::min(arraySize, static_cast<uint16_t>(initList.size())); std::min(arraySize, static_cast<uint16_t>(initList.size()));
@@ -41,13 +32,11 @@ Matrix<rows, columns>::Matrix(Args... args) {
} }
template <uint8_t rows, uint8_t columns> template <uint8_t rows, uint8_t columns>
Matrix<rows, columns> Matrix<rows, columns>::Identity() { void Matrix<rows, columns>::Identity() {
Matrix<rows, columns> identityMatrix{0}; this->Fill(0);
uint32_t minDimension = std::min(rows, columns); for (uint8_t idx{0}; idx < rows; idx++) {
for (uint8_t idx{0}; idx < minDimension; idx++) { this->matrix[idx * columns + idx] = 1;
identityMatrix[idx][idx] = 1;
} }
return identityMatrix;
} }
template <uint8_t rows, uint8_t columns> template <uint8_t rows, uint8_t columns>
@@ -544,12 +533,7 @@ void Matrix<rows, columns>::QRDecomposition(Matrix<rows, columns> &Q,
Matrix<1, rows> Q_column_k_T = Q_column_k.Transpose(); Matrix<1, rows> Q_column_k_T = Q_column_k.Transpose();
u = u - Q_column_k * (Q_column_k_T * a_col); u = u - Q_column_k * (Q_column_k_T * a_col);
} }
float norm = u.EuclideanNorm(); u = u / u.EuclideanNorm();
if (norm > 1e-4) {
u = u / norm;
} else {
u.Fill(0);
}
Q.SetSubMatrix(0, column, u); Q.SetSubMatrix(0, column, u);
// ----------------------- // -----------------------
@@ -568,19 +552,16 @@ void Matrix<rows, columns>::EigenQR(Matrix<rows, rows> &eigenVectors,
uint32_t maxIterations, uint32_t maxIterations,
float tolerance) const { float tolerance) const {
static_assert(rows > 1, "Matrix size must be > 1 for QR iteration"); 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> Ak = *this; // Copy original matrix
Matrix<rows, rows> QQ{Matrix<rows, rows>::Identity()}; Matrix<rows, rows> QQ{};
Matrix<rows, rows> shift{0}; QQ.Identity();
for (uint32_t iter = 0; iter < maxIterations; ++iter) { for (uint32_t iter = 0; iter < maxIterations; ++iter) {
Matrix<rows, rows> Q, R; Matrix<rows, rows> Q, R;
Ak.QRDecomposition(Q, R);
// // QR shift lets us "attack" the first diagonal to speed up the algorithm Ak = R * Q;
// shift = Matrix<rows, rows>::Identity() * Ak[rows - 1][rows - 1];
(Ak - shift).QRDecomposition(Q, R);
Ak = R * Q + shift;
QQ = QQ * Q; QQ = QQ * Q;
// Check convergence: off-diagonal norm // Check convergence: off-diagonal norm

View File

@@ -1,4 +1,5 @@
#pragma once #ifndef MATRIX_H_
#define MATRIX_H_
#include <array> #include <array>
#include <cstdint> #include <cstdint>
@@ -18,6 +19,11 @@ public:
*/ */
Matrix() = default; 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 * @brief Initialize a matrix with an array
*/ */
@@ -34,9 +40,9 @@ public:
template <typename... Args> Matrix(Args... args); template <typename... Args> Matrix(Args... args);
/** /**
* @brief Create an identity matrix * @brief set the matrix diagonals to 1 and all other values to 0
*/ */
static Matrix<rows, columns> Identity(); void Identity();
/** /**
* @brief Set all elements in this to value * @brief Set all elements in this to value
@@ -252,6 +258,6 @@ private:
void setMatrixToArray(const std::array<float, rows * columns> &array); void setMatrixToArray(const std::array<float, rows * columns> &array);
}; };
#ifndef MATRIX_H_
#include "Matrix.cpp" #include "Matrix.cpp"
#endif // MATRIX_H_ #endif // MATRIX_H_

View File

@@ -2,11 +2,12 @@
#define QUATERNION_H_ #define QUATERNION_H_
#include "Matrix.hpp" #include "Matrix.hpp"
class Quaternion : public Matrix<1, 4> { class Quaternion : public Matrix<1, 4>
{
public: public:
Quaternion() : Matrix<1, 4>() {} Quaternion() : Matrix<1, 4>() {}
Quaternion(float w, float v1, float v2, float v3) Quaternion(float fillValue) : Matrix<1, 4>(fillValue) {}
: Matrix<1, 4>(w, v1, v2, v3) {} 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 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 Matrix<1, 4> &matrix) : Matrix<1, 4>(matrix) {}
Quaternion(const std::array<float, 4> &array) : Matrix<1, 4>(array) {} Quaternion(const std::array<float, 4> &array) : Matrix<1, 4>(array) {}

View File

@@ -10,61 +10,41 @@
#include <cmath> #include <cmath>
#include <iostream> #include <iostream>
// Helper functions
template <uint8_t rows, uint8_t columns>
float matrixSum(const Matrix<rows, columns> &matrix) {
float sum = 0;
for (uint32_t i = 0; i < rows * columns; i++) {
float number = matrix.ToArray()[i];
sum += number * number;
}
return std::sqrt(sum);
}
template <uint8_t rows, uint8_t columns>
void printLabeledMatrix(const std::string &label,
const Matrix<rows, columns> &matrix) {
std::string strBuf = "";
matrix.ToString(strBuf);
std::cout << label << ":\n" << strBuf << std::endl;
}
TEST_CASE("Initialization", "Matrix") {
SECTION("Array Initialization") {
std::array<float, 4> arr2{5, 6, 7, 8};
Matrix<2, 2> mat2{arr2};
// array initialization
REQUIRE(mat2.Get(0, 0) == 5);
REQUIRE(mat2.Get(0, 1) == 6);
REQUIRE(mat2.Get(1, 0) == 7);
REQUIRE(mat2.Get(1, 1) == 8);
}
SECTION("Argument Pack Initialization") {
Matrix<2, 2> mat1{1, 2, 3, 4};
// template pack initialization
REQUIRE(mat1.Get(0, 0) == 1);
REQUIRE(mat1.Get(0, 1) == 2);
REQUIRE(mat1.Get(1, 0) == 3);
REQUIRE(mat1.Get(1, 1) == 4);
}
SECTION("Single Argument Pack Initialization") {
Matrix<2, 2> mat1{2};
// template pack initialization
REQUIRE(mat1.Get(0, 0) == 2);
REQUIRE(mat1.Get(0, 1) == 2);
REQUIRE(mat1.Get(1, 0) == 2);
REQUIRE(mat1.Get(1, 1) == 2);
}
}
TEST_CASE("Elementary Matrix Operations", "Matrix") { TEST_CASE("Elementary Matrix Operations", "Matrix") {
std::array<float, 4> arr2{5, 6, 7, 8}; std::array<float, 4> arr2{5, 6, 7, 8};
Matrix<2, 2> mat1{1, 2, 3, 4}; Matrix<2, 2> mat1{1, 2, 3, 4};
Matrix<2, 2> mat2{arr2}; Matrix<2, 2> mat2{arr2};
Matrix<2, 2> mat3{}; Matrix<2, 2> mat3{};
SECTION("Initialization") {
// array initialization
REQUIRE(mat1.Get(0, 0) == 1);
REQUIRE(mat1.Get(0, 1) == 2);
REQUIRE(mat1.Get(1, 0) == 3);
REQUIRE(mat1.Get(1, 1) == 4);
// empty initialization
REQUIRE(mat3.Get(0, 0) == 0);
REQUIRE(mat3.Get(0, 1) == 0);
REQUIRE(mat3.Get(1, 0) == 0);
REQUIRE(mat3.Get(1, 1) == 0);
// template pack initialization
REQUIRE(mat2.Get(0, 0) == 5);
REQUIRE(mat2.Get(0, 1) == 6);
REQUIRE(mat2.Get(1, 0) == 7);
REQUIRE(mat2.Get(1, 1) == 8);
// large matrix
Matrix<255, 255> mat6{};
mat6.Fill(4);
for (uint8_t row{0}; row < 255; row++) {
for (uint8_t column{0}; column < 255; column++) {
REQUIRE(mat6.Get(row, column) == 4);
}
}
}
SECTION("Fill") { SECTION("Fill") {
mat1.Fill(0); mat1.Fill(0);
REQUIRE(mat1.Get(0, 0) == 0); REQUIRE(mat1.Get(0, 0) == 0);
@@ -86,6 +66,10 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
} }
SECTION("Addition") { SECTION("Addition") {
std::string strBuf1 = "";
mat1.ToString(strBuf1);
std::cout << "Matrix 1:\n" << strBuf1 << std::endl;
mat1.Add(mat2, mat3); mat1.Add(mat2, mat3);
REQUIRE(mat3.Get(0, 0) == 6); REQUIRE(mat3.Get(0, 0) == 6);
@@ -352,86 +336,46 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
Matrix<3, 3> mat4 = startMatrix; Matrix<3, 3> mat4 = startMatrix;
Matrix<2, 2> mat5{10, 11, 12, 13}; 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, 0) == 10);
REQUIRE(mat4.Get(0, 1) == 11); REQUIRE(mat4.Get(0, 1) == 11);
REQUIRE(mat4.Get(1, 0) == 12); REQUIRE(mat4.Get(1, 0) == 12);
REQUIRE(mat4.Get(1, 1) == 13); REQUIRE(mat4.Get(1, 1) == 13);
mat4 = startMatrix; mat4 = startMatrix;
mat4.SetSubMatrix(1, 1, mat5); mat4.SetSubMatrix<2, 2, 1, 1>(mat5);
REQUIRE(mat4.Get(1, 1) == 10); REQUIRE(mat4.Get(1, 1) == 10);
REQUIRE(mat4.Get(1, 2) == 11); REQUIRE(mat4.Get(1, 2) == 11);
REQUIRE(mat4.Get(2, 1) == 12); REQUIRE(mat4.Get(2, 1) == 12);
REQUIRE(mat4.Get(2, 2) == 13); REQUIRE(mat4.Get(2, 2) == 13);
Matrix<3, 1> mat6{10, 11, 12}; 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(0, 0) == 10);
REQUIRE(mat4.Get(1, 0) == 11); REQUIRE(mat4.Get(1, 0) == 11);
REQUIRE(mat4.Get(2, 0) == 12); REQUIRE(mat4.Get(2, 0) == 12);
Matrix<1, 3> mat7{10, 11, 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, 0) == 10);
REQUIRE(mat4.Get(0, 1) == 11); REQUIRE(mat4.Get(0, 1) == 11);
REQUIRE(mat4.Get(0, 2) == 12); REQUIRE(mat4.Get(0, 2) == 12);
} }
} }
TEST_CASE("Identity Matrix", "Matrix") { template <uint8_t rows, uint8_t columns>
SECTION("Square Matrix") { float matrixSum(const Matrix<rows, columns> &matrix) {
Matrix<5, 5> matrix = Matrix<5, 5>::Identity(); float sum = 0;
uint32_t oneColumnIndex{0}; for (uint32_t i = 0; i < rows * columns; i++) {
for (uint32_t row = 0; row < 5; row++) { float number = matrix.ToArray()[i];
for (uint32_t column = 0; column < 5; column++) { sum += number * number;
float value = matrix[row][column];
if (oneColumnIndex == column) {
REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
} else {
REQUIRE_THAT(value, Catch::Matchers::WithinRel(0.0f, 1e-6f));
}
}
oneColumnIndex++;
}
}
SECTION("Wide Matrix") {
Matrix<2, 5> matrix = Matrix<2, 5>::Identity();
uint32_t oneColumnIndex{0};
for (uint32_t row = 0; row < 2; row++) {
for (uint32_t column = 0; column < 5; column++) {
float value = matrix[row][column];
if (oneColumnIndex == column && row < 3) {
REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
} else {
REQUIRE_THAT(value, Catch::Matchers::WithinRel(0.0f, 1e-6f));
}
}
oneColumnIndex++;
}
}
SECTION("Tall Matrix") {
Matrix<5, 2> matrix = Matrix<5, 2>::Identity();
uint32_t oneColumnIndex{0};
for (uint32_t row = 0; row < 5; row++) {
for (uint32_t column = 0; column < 2; column++) {
float value = matrix[row][column];
if (oneColumnIndex == column) {
REQUIRE_THAT(value, Catch::Matchers::WithinRel(1.0f, 1e-6f));
} else {
REQUIRE_THAT(value, Catch::Matchers::WithinRel(0.0f, 1e-6f));
}
}
oneColumnIndex++;
}
} }
return std::sqrt(sum);
} }
// TODO: Add test for scalar division // TODO: Add test for scalar division
TEST_CASE("Euclidean Norm", "Matrix") {
TEST_CASE("Normalization", "Matrix") {
SECTION("2x2 Normalize") { SECTION("2x2 Normalize") {
Matrix<2, 2> mat1{1, 2, 3, 4}; Matrix<2, 2> mat1{1, 2, 3, 4};
@@ -529,6 +473,13 @@ TEST_CASE("QR Decompositions", "Matrix") {
Matrix<3, 3> Q{}, R{}; Matrix<3, 3> Q{}, R{};
A.QRDecomposition(Q, R); A.QRDecomposition(Q, R);
std::string strBuf1 = "";
Q.ToString(strBuf1);
std::cout << "Q:\n" << strBuf1 << std::endl;
strBuf1 = "";
R.ToString(strBuf1);
std::cout << "R:\n" << strBuf1 << std::endl;
// Check that Q * R ≈ A // Check that Q * R ≈ A
Matrix<3, 3> QR{}; Matrix<3, 3> QR{};
QR = Q * R; QR = Q * R;
@@ -539,13 +490,11 @@ TEST_CASE("QR Decompositions", "Matrix") {
} }
// Check that Qᵀ * Q ≈ I // Check that Qᵀ * Q ≈ I
// Since the rank of this matrix is 2, only the top left 2x2 sub-matrix will
// equal I.
Matrix<3, 3> Qt = Q.Transpose(); Matrix<3, 3> Qt = Q.Transpose();
Matrix<3, 3> QtQ{}; Matrix<3, 3> QtQ{};
QtQ = Qt * Q; QtQ = Qt * Q;
for (int i = 0; i < 2; ++i) { for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 2; ++j) { for (int j = 0; j < 3; ++j) {
if (i == j) if (i == j)
REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f)); REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f));
else else
@@ -559,6 +508,28 @@ TEST_CASE("QR Decompositions", "Matrix") {
REQUIRE(std::fabs(R[i][j]) < 1e-4f); REQUIRE(std::fabs(R[i][j]) < 1e-4f);
} }
} }
// check that all Q values are correct
REQUIRE_THAT(Q[0][0], Catch::Matchers::WithinRel(0.1231f, 1e-4f));
REQUIRE_THAT(Q[0][1], Catch::Matchers::WithinRel(0.904534f, 1e-4f));
REQUIRE_THAT(Q[0][2], Catch::Matchers::WithinRel(0.0f, 1e-4f));
REQUIRE_THAT(Q[1][0], Catch::Matchers::WithinRel(0.49237f, 1e-4f));
REQUIRE_THAT(Q[1][1], Catch::Matchers::WithinRel(0.301511f, 1e-4f));
REQUIRE_THAT(Q[1][2], Catch::Matchers::WithinRel(0.0f, 1e-4f));
REQUIRE_THAT(Q[2][0], Catch::Matchers::WithinRel(0.86164f, 1e-4f));
REQUIRE_THAT(Q[2][1], Catch::Matchers::WithinRel(-0.30151f, 1e-4f));
REQUIRE_THAT(Q[2][2], Catch::Matchers::WithinRel(0.0f, 1e-4f));
// check that all R values are correct
REQUIRE_THAT(R[0][0], Catch::Matchers::WithinRel(8.124038f, 1e-4f));
REQUIRE_THAT(R[0][1], Catch::Matchers::WithinRel(9.60114f, 1e-4f));
REQUIRE_THAT(R[0][2], Catch::Matchers::WithinRel(11.07823f, 1e-4f));
REQUIRE_THAT(R[1][0], Catch::Matchers::WithinRel(0.0f, 1e-4f));
REQUIRE_THAT(R[1][1], Catch::Matchers::WithinRel(0.90453f, 1e-4f));
REQUIRE_THAT(R[1][2], Catch::Matchers::WithinRel(1.80907f, 1e-4f));
REQUIRE_THAT(R[2][0], Catch::Matchers::WithinRel(0.0f, 1e-4f));
REQUIRE_THAT(R[2][1], Catch::Matchers::WithinRel(0.0f, 1e-4f));
REQUIRE_THAT(R[2][2], Catch::Matchers::WithinRel(1.0f, 1e-4f));
} }
SECTION("4x2 QRDecomposition") { SECTION("4x2 QRDecomposition") {
@@ -600,41 +571,42 @@ TEST_CASE("QR Decompositions", "Matrix") {
} }
} }
TEST_CASE("Eigenvalues and Vectors", "Matrix") { // TEST_CASE("Eigenvalues and Vectors", "Matrix") {
SECTION("2x2 Eigen") { // SECTION("2x2 Eigen") {
Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f}; // Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
Matrix<2, 2> vectors{}; // Matrix<2, 2> vectors{};
Matrix<2, 1> values{}; // Matrix<2, 1> values{};
A.EigenQR(vectors, values, 1000000, 1e-20f); // A.EigenQR(vectors, values, 1000000, 1e-20f);
REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.41597f, 1e-4f)); // REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.41597f, 1e-4f));
REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.90938f, 1e-4f)); // REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.90938f, 1e-4f));
REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(5.372282f, 1e-4f)); // REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(5.372282f, 1e-4f));
REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(-0.372281f, 1e-4f)); // REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(-0.372281f,
} // 1e-4f));
// }
SECTION("3x3 Rank Defficient Eigen") { // SECTION("3x3 Eigen") {
SKIP("Skipping this because QR decomposition isn't ready for it"); // // this symmetrix tridiagonal matrix is well behaved for testing
// this symmetrix tridiagonal matrix is well behaved for testing // Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
Matrix<3, 3> vectors{}; // Matrix<3, 3> vectors{};
Matrix<3, 1> values{}; // Matrix<3, 1> values{};
A.EigenQR(vectors, values, 1000000, 1e-8f); // A.EigenQR(vectors, values, 1000000, 1e-8f);
std::string strBuf1 = ""; // std::string strBuf1 = "";
vectors.ToString(strBuf1); // vectors.ToString(strBuf1);
std::cout << "Vectors:\n" << strBuf1 << std::endl; // std::cout << "Vectors:\n" << strBuf1 << std::endl;
strBuf1 = ""; // strBuf1 = "";
values.ToString(strBuf1); // values.ToString(strBuf1);
std::cout << "Values:\n" << strBuf1 << std::endl; // std::cout << "Values:\n" << strBuf1 << std::endl;
REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.23197f, 1e-4f)); // REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.23197f, 1e-4f));
REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.525322f, 1e-4f)); // REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.525322f,
REQUIRE_THAT(vectors[2][0], Catch::Matchers::WithinRel(0.81867f, 1e-4f)); // 1e-4f)); REQUIRE_THAT(vectors[2][0], Catch::Matchers::WithinRel(0.81867f,
REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(-1.11684f, 1e-4f)); // 1e-4f)); REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(-1.11684f,
REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(0.0f, 1e-4f)); // 1e-4f)); REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(0.0f,
REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(16.1168f, 1e-4f)); // 1e-4f)); REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(16.1168f,
} // 1e-4f));
} // }
// }

View File

@@ -8,7 +8,6 @@
// any other libraries // any other libraries
#include <array> #include <array>
#include <cmath> #include <cmath>
#include <cstdint>
// basically re-run all of the matrix tests with huge matrices and time the // basically re-run all of the matrix tests with huge matrices and time the
// results. // results.
@@ -30,13 +29,13 @@ TEST_CASE("Timing Tests", "Matrix") {
Matrix<4, 4> mat5{}; Matrix<4, 4> mat5{};
SECTION("Addition") { SECTION("Addition") {
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat3 = mat1 + mat2; mat3 = mat1 + mat2;
} }
} }
SECTION("Subtraction") { SECTION("Subtraction") {
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat3 = mat1 - mat2; mat3 = mat1 - mat2;
} }
} }
@@ -48,19 +47,19 @@ TEST_CASE("Timing Tests", "Matrix") {
} }
SECTION("Scalar Multiplication") { SECTION("Scalar Multiplication") {
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat3 = mat1 * 3; mat3 = mat1 * 3;
} }
} }
SECTION("Element Multiply") { SECTION("Element Multiply") {
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat1.ElementMultiply(mat2, mat3); mat1.ElementMultiply(mat2, mat3);
} }
} }
SECTION("Element Divide") { SECTION("Element Divide") {
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat1.ElementDivide(mat2, mat3); mat1.ElementDivide(mat2, mat3);
} }
} }
@@ -69,59 +68,52 @@ TEST_CASE("Timing Tests", "Matrix") {
// what about matrices of 0,0 or 1,1? // what about matrices of 0,0 or 1,1?
// minor matrix for 2x2 matrix // minor matrix for 2x2 matrix
Matrix<49, 49> minorMat1{}; Matrix<49, 49> minorMat1{};
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat1.MinorMatrix(minorMat1, 0, 0); mat1.MinorMatrix(minorMat1, 0, 0);
} }
} }
SECTION("Determinant") { SECTION("Determinant") {
for (uint32_t i{0}; i < 1000000; i++) { for (uint32_t i{0}; i < 100000; i++) {
float det1 = mat4.Det(); float det1 = mat4.Det();
} }
} }
SECTION("Matrix of Minors") { SECTION("Matrix of Minors") {
for (uint32_t i{0}; i < 1000000; i++) { for (uint32_t i{0}; i < 100000; i++) {
mat4.MatrixOfMinors(mat5); mat4.MatrixOfMinors(mat5);
} }
} }
SECTION("Invert") { SECTION("Invert") {
for (uint32_t i{0}; i < 1000000; i++) { for (uint32_t i{0}; i < 100000; i++) {
mat5 = mat4.Invert(); mat5 = mat4.Invert();
} }
}; };
SECTION("Transpose") { SECTION("Transpose") {
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat3 = mat1.Transpose(); mat3 = mat1.Transpose();
} }
} }
SECTION("Normalize") { SECTION("Normalize") {
for (uint32_t i{0}; i < 100000; i++) { for (uint32_t i{0}; i < 10000; i++) {
mat3 = mat1 / mat1.EuclideanNorm(); mat3 = mat1 / mat1.EuclideanNorm();
} }
} }
SECTION("GET ROW") { SECTION("GET ROW") {
Matrix<1, 50> mat1Rows{}; Matrix<1, 50> mat1Rows{};
for (uint32_t i{0}; i < 100000000; i++) { for (uint32_t i{0}; i < 1000000; i++) {
mat1.GetRow(0, mat1Rows); mat1.GetRow(0, mat1Rows);
} }
} }
SECTION("GET COLUMN") { SECTION("GET COLUMN") {
Matrix<50, 1> mat1Columns{}; Matrix<50, 1> mat1Columns{};
for (uint32_t i{0}; i < 100000000; i++) { for (uint32_t i{0}; i < 1000000; i++) {
mat1.GetColumn(0, mat1Columns); mat1.GetColumn(0, mat1Columns);
} }
} }
SECTION("QR Decomposition") {
Matrix<50, 50> Q, R{};
for (uint32_t i{0}; i < 500; i++) {
mat1.QRDecomposition(Q, R);
}
}
} }

View File

@@ -1,36 +1,56 @@
Running matrix-timing-tests with timing Randomness seeded to: 2444679151
Randomness seeded to: 3567651885 0.180 s: Addition
1.857 s: Addition 0.180 s: Timing Tests
1.857 s: Timing Tests 0.177 s: Subtraction
1.788 s: Subtraction 0.177 s: Timing Tests
1.788 s: Timing Tests 1.868 s: Multiplication
1.929 s: Multiplication 1.868 s: Timing Tests
1.929 s: Timing Tests 0.127 s: Scalar Multiplication
1.268 s: Scalar Multiplication 0.127 s: Timing Tests
1.268 s: Timing Tests 0.173 s: Element Multiply
1.798 s: Element Multiply 0.173 s: Timing Tests
1.798 s: Timing Tests 0.178 s: Element Divide
1.802 s: Element Divide 0.178 s: Timing Tests
1.803 s: Timing Tests 0.172 s: Minor Matrix
1.553 s: Minor Matrix 0.172 s: Timing Tests
1.554 s: Timing Tests 0.103 s: Determinant
1.009 s: Determinant 0.103 s: Timing Tests
1.009 s: Timing Tests 0.411 s: Matrix of Minors
4.076 s: Matrix of Minors 0.411 s: Timing Tests
4.076 s: Timing Tests 0.109 s: Invert
1.066 s: Invert 0.109 s: Timing Tests
1.066 s: Timing Tests 0.122 s: Transpose
1.246 s: Transpose 0.122 s: Timing Tests
1.246 s: Timing Tests 0.190 s: Normalize
2.284 s: Normalize 0.190 s: Timing Tests
2.284 s: Timing Tests 0.006 s: GET ROW
0.606 s: GET ROW 0.006 s: Timing Tests
0.606 s: Timing Tests 0.235 s: GET COLUMN
24.629 s: GET COLUMN 0.235 s: Timing Tests
24.630 s: Timing Tests
3.064 s: QR Decomposition
3.064 s: Timing Tests
=============================================================================== ===============================================================================
test cases: 1 | 1 passed test cases: 1 | 1 passed
assertions: - none - assertions: - none -
Command being timed: "build/unit-tests/matrix-timing-tests -d yes"
User time (seconds): 4.05
System time (seconds): 0.00
Percent of CPU this job got: 100%
Elapsed (wall clock) time (h:mm:ss or m:ss): 0:04.05
Average shared text size (kbytes): 0
Average unshared data size (kbytes): 0
Average stack size (kbytes): 0
Average total size (kbytes): 0
Maximum resident set size (kbytes): 3200
Average resident set size (kbytes): 0
Major (requiring I/O) page faults: 184
Minor (reclaiming a frame) page faults: 171
Voluntary context switches: 1
Involuntary context switches: 26
Swaps: 0
File system inputs: 12
File system outputs: 1
Socket messages sent: 0
Socket messages received: 0
Signals delivered: 0
Page size (bytes): 4096
Exit status: 0