From c07e9b265213a847670c24a094de1735d1147186 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 5 Oct 2020 18:39:33 +0200
Subject: [PATCH] Fix TinyMatrix determinant for N>3
The case of zero determinant was not treated correctly when dealing
with N*N matrices with N>4
For these matrices, a Gauss-pivot is used but it was incorrect in the
singular case
---
src/algebra/TinyMatrix.hpp | 31 ++++++++++++++++++++-----------
tests/test_TinyMatrix.cpp | 3 +++
2 files changed, 23 insertions(+), 11 deletions(-)
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 3ef5e1807..a4499baa4 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 e0d7ef39d..0d3510dc3 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")
--
GitLab