diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index e34f98ed8c25df7a77ba90bbd7cd05305b5eb115..1717b48c78deaedc7cabd9f60feb215e29a80ebe 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -21,16 +21,21 @@ class [[nodiscard]] TinyMatrix
using data_type = T;
private:
+ static_assert((M > 0), "TinyMatrix number of rows must be strictly positive");
+ static_assert((N > 0), "TinyMatrix number of columns must be strictly positive");
+
T m_values[M * N];
PUGS_FORCEINLINE
- constexpr size_t _index(size_t i, size_t j) const noexcept // LCOV_EXCL_LINE (due to forced inline)
+ constexpr size_t
+ _index(size_t i, size_t j) const noexcept // LCOV_EXCL_LINE (due to forced inline)
{
return i * N + j;
}
template <typename... Args>
- PUGS_FORCEINLINE constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
+ PUGS_FORCEINLINE constexpr void
+ _unpackVariadicInput(const T& t, Args&&... args) noexcept
{
m_values[M * N - 1 - sizeof...(args)] = t;
if constexpr (sizeof...(args) > 0) {
@@ -40,13 +45,15 @@ class [[nodiscard]] TinyMatrix
public:
PUGS_INLINE
- constexpr bool isSquare() const noexcept
+ constexpr bool
+ isSquare() const noexcept
{
return M == N;
}
PUGS_INLINE
- constexpr friend TinyMatrix<N, M, T> transpose(const TinyMatrix& A)
+ constexpr friend TinyMatrix<N, M, T>
+ transpose(const TinyMatrix& A)
{
TinyMatrix<N, M, T> tA;
for (size_t i = 0; i < M; ++i) {
@@ -58,31 +65,36 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr size_t dimension() const
+ constexpr size_t
+ dimension() const
{
return M * N;
}
PUGS_INLINE
- constexpr size_t numberOfValues() const
+ constexpr size_t
+ numberOfValues() const
{
return this->dimension();
}
PUGS_INLINE
- constexpr size_t numberOfRows() const
+ constexpr size_t
+ numberOfRows() const
{
return M;
}
PUGS_INLINE
- constexpr size_t numberOfColumns() const
+ constexpr size_t
+ numberOfColumns() const
{
return N;
}
PUGS_INLINE
- constexpr TinyMatrix operator-() const
+ constexpr TinyMatrix
+ operator-() const
{
TinyMatrix opposite;
for (size_t i = 0; i < M * N; ++i) {
@@ -92,20 +104,23 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr friend TinyMatrix operator*(const T& t, const TinyMatrix& A)
+ constexpr friend TinyMatrix
+ operator*(const T& t, const TinyMatrix& A)
{
TinyMatrix B = A;
return B *= t;
}
PUGS_INLINE
- constexpr friend TinyMatrix operator*(const T& t, TinyMatrix&& A)
+ constexpr friend TinyMatrix
+ operator*(const T& t, TinyMatrix&& A)
{
return std::move(A *= t);
}
PUGS_INLINE
- constexpr TinyMatrix& operator*=(const T& t)
+ constexpr TinyMatrix&
+ operator*=(const T& t)
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] *= t;
@@ -114,7 +129,8 @@ class [[nodiscard]] TinyMatrix
}
template <size_t P>
- PUGS_INLINE constexpr TinyMatrix<M, P, T> operator*(const TinyMatrix<N, P, T>& B) const
+ PUGS_INLINE constexpr TinyMatrix<M, P, T>
+ operator*(const TinyMatrix<N, P, T>& B) const
{
const TinyMatrix& A = *this;
TinyMatrix<M, P, T> AB;
@@ -131,7 +147,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyVector<M, T> operator*(const TinyVector<N, T>& x) const
+ constexpr TinyVector<M, T>
+ operator*(const TinyVector<N, T>& x) const
{
const TinyMatrix& A = *this;
TinyVector<M, T> Ax;
@@ -146,7 +163,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
+ constexpr friend std::ostream&
+ operator<<(std::ostream& os, const TinyMatrix& A)
{
os << '[';
for (size_t i = 0; i < M; ++i) {
@@ -165,7 +183,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr bool operator==(const TinyMatrix& A) const
+ constexpr bool
+ operator==(const TinyMatrix& A) const
{
for (size_t i = 0; i < M * N; ++i) {
if (m_values[i] != A.m_values[i])
@@ -175,13 +194,15 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr bool operator!=(const TinyMatrix& A) const
+ constexpr bool
+ operator!=(const TinyMatrix& A) const
{
return not this->operator==(A);
}
PUGS_INLINE
- constexpr TinyMatrix operator+(const TinyMatrix& A) const
+ constexpr TinyMatrix
+ operator+(const TinyMatrix& A) const
{
TinyMatrix sum;
for (size_t i = 0; i < M * N; ++i) {
@@ -191,14 +212,16 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix operator+(TinyMatrix&& A) const
+ constexpr TinyMatrix
+ operator+(TinyMatrix&& A) const
{
A += *this;
return std::move(A);
}
PUGS_INLINE
- constexpr TinyMatrix operator-(const TinyMatrix& A) const
+ constexpr TinyMatrix
+ operator-(const TinyMatrix& A) const
{
TinyMatrix difference;
for (size_t i = 0; i < M * N; ++i) {
@@ -208,7 +231,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix operator-(TinyMatrix&& A) const
+ constexpr TinyMatrix
+ operator-(TinyMatrix&& A) const
{
for (size_t i = 0; i < M * N; ++i) {
A.m_values[i] = m_values[i] - A.m_values[i];
@@ -217,7 +241,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix& operator+=(const TinyMatrix& A)
+ constexpr TinyMatrix&
+ operator+=(const TinyMatrix& A)
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
@@ -226,7 +251,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr void operator+=(const volatile TinyMatrix& A) volatile
+ constexpr void
+ operator+=(const volatile TinyMatrix& A) volatile
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] += A.m_values[i];
@@ -234,7 +260,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix& operator-=(const TinyMatrix& A)
+ constexpr TinyMatrix&
+ operator-=(const TinyMatrix& A)
{
for (size_t i = 0; i < M * N; ++i) {
m_values[i] -= A.m_values[i];
@@ -243,21 +270,24 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr T& operator()(size_t i, size_t j) noexcept(NO_ASSERT)
+ constexpr T&
+ operator()(size_t i, size_t j) noexcept(NO_ASSERT)
{
Assert((i < M) and (j < N));
return m_values[_index(i, j)];
}
PUGS_INLINE
- constexpr const T& operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
+ constexpr const T&
+ operator()(size_t i, size_t j) const noexcept(NO_ASSERT)
{
Assert((i < M) and (j < N));
return m_values[_index(i, j)];
}
PUGS_INLINE
- constexpr TinyMatrix& operator=(ZeroType) noexcept
+ constexpr TinyMatrix&
+ operator=(ZeroType) noexcept
{
static_assert(std::is_arithmetic<T>(), "Cannot assign 'zero' value for non-arithmetic types");
for (size_t i = 0; i < M * N; ++i) {
@@ -267,7 +297,8 @@ class [[nodiscard]] TinyMatrix
}
PUGS_INLINE
- constexpr TinyMatrix& operator=(IdentityType) noexcept
+ constexpr TinyMatrix&
+ operator=(IdentityType) noexcept
{
static_assert(std::is_arithmetic<T>(), "Cannot assign 'identity' value for non-arithmetic types");
for (size_t i = 0; i < M; ++i) {
@@ -329,7 +360,7 @@ class [[nodiscard]] TinyMatrix
constexpr TinyMatrix(const TinyMatrix&) noexcept = default;
PUGS_INLINE
- TinyMatrix(TinyMatrix && A) noexcept = default;
+ TinyMatrix(TinyMatrix&& A) noexcept = default;
PUGS_INLINE
~TinyMatrix() = default;
@@ -450,6 +481,19 @@ getMinor(const TinyMatrix<M, N, T>& A, size_t I, size_t J)
return m;
}
+template <size_t N, typename T>
+PUGS_INLINE T
+trace(const TinyMatrix<N, N, T>& A)
+{
+ static_assert(std::is_arithmetic<T>::value, "trace is not defined for non-arithmetic types");
+
+ T t = A(0, 0);
+ for (size_t i = 1; i < N; ++i) {
+ t += A(i, i);
+ }
+ return t;
+}
+
template <size_t N, typename T>
PUGS_INLINE constexpr TinyMatrix<N, N, T> inverse(const TinyMatrix<N, N, T>& A);
diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp
index 2f4c09d4ad9ce1dbb5e3d856d8737d27344459f2..f6f09cf55dd2ad1ffef55c9275b989c203598f7e 100644
--- a/tests/test_TinyMatrix.cpp
+++ b/tests/test_TinyMatrix.cpp
@@ -194,6 +194,17 @@ TEST_CASE("TinyMatrix", "[algebra]")
REQUIRE(det(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0);
}
+ SECTION("checking for trace calculations")
+ {
+ REQUIRE(trace(TinyMatrix<1, 1, int>(6)) == 6);
+ REQUIRE(trace(TinyMatrix<2, 2, int>(5, 1, -3, 6)) == 5 + 6);
+ REQUIRE(trace(TinyMatrix<3, 3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1 + 2 + 3);
+ REQUIRE(trace(TinyMatrix<3, 3, int>(6, 5, 3, 8, 34, 6, 35, 6, 7)) == 6 + 34 + 7);
+ REQUIRE(trace(TinyMatrix<4>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) ==
+ Catch::Approx(1 + 4 + 2 + 17.5).epsilon(1E-14));
+ REQUIRE(trace(TinyMatrix<4>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 1 + 0 + 1 + 2);
+ }
+
SECTION("checking for inverse calculations")
{
{