diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp index 3ef5e180710d90e0b62e0cf8f8d222cb3a84ca14..a4499baa4677f09bc7edd3c1e3c59d41c9b3557d 100644 --- a/src/algebra/TinyMatrix.hpp +++ b/src/algebra/TinyMatrix.hpp @@ -50,14 +50,16 @@ class 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); } @@ -73,7 +75,8 @@ class TinyMatrix } PUGS_INLINE - constexpr TinyMatrix operator*(const TinyMatrix& B) const + constexpr TinyMatrix + operator*(const TinyMatrix& B) const { const TinyMatrix& A = *this; TinyMatrix AB; @@ -90,7 +93,8 @@ class TinyMatrix } PUGS_INLINE - constexpr TinyVector<N, T> operator*(const TinyVector<N, T>& x) const + constexpr TinyVector<N, T> + operator*(const TinyVector<N, T>& x) const { const TinyMatrix& A = *this; TinyVector<N, T> Ax; @@ -342,14 +346,19 @@ det(const TinyMatrix<N, T>& A) determinent *= -1; } } - const size_t I = index[i]; - const T inv_Mii = 1. / M(I, i); - for (size_t k = i + 1; k < N; ++k) { - const size_t K = index[k]; - const T factor = M(K, i) * inv_Mii; - for (size_t l = i + 1; l < N; ++l) { - M(K, l) -= factor * M(I, l); + const size_t I = index[i]; + const T Mii = M(I, i); + if (Mii != 0) { + const T inv_Mii = 1. / M(I, i); + for (size_t k = i + 1; k < N; ++k) { + const size_t K = index[k]; + const T factor = M(K, i) * inv_Mii; + for (size_t l = i + 1; l < N; ++l) { + M(K, l) -= factor * M(I, l); + } } + } else { + return 0; } } diff --git a/tests/test_TinyMatrix.cpp b/tests/test_TinyMatrix.cpp index e0d7ef39da53ae03358b307e448136d69216d1cf..0d3510dc361a9fd43b1a17559b514c2c3748bc08 100644 --- a/tests/test_TinyMatrix.cpp +++ b/tests/test_TinyMatrix.cpp @@ -161,9 +161,12 @@ TEST_CASE("TinyMatrix", "[algebra]") { REQUIRE(det(TinyMatrix<1, int>(6)) == 6); REQUIRE(det(TinyMatrix<2, int>(3, 1, -3, 6)) == 21); + REQUIRE(det(TinyMatrix<3, int>(1, 1, 1, 1, 2, 1, 2, 1, 3)) == 1); REQUIRE(det(B) == -1444); REQUIRE(det(TinyMatrix<4, double>(1, 2.3, 7, -6.2, 3, 4, 9, 1, 4.1, 5, 2, -3, 2, 27, 3, 17.5)) == Approx(6661.455).epsilon(1E-14)); + + REQUIRE(det(TinyMatrix<4, double>(1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 2, 0, 0, 2, 2)) == 0); } SECTION("checking for inverse calculations")