Skip to content
Snippets Groups Projects
Commit 14615ef6 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Merge branch 'issue/tinymatrix_det' into 'develop'

Fix TinyMatrix determinant for N>3

See merge request !53
parents 3feead36 c07e9b26
No related branches found
No related tags found
1 merge request!53Fix TinyMatrix determinant for N>3
......@@ -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;
......@@ -343,6 +347,8 @@ det(const TinyMatrix<N, T>& A)
}
}
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];
......@@ -351,6 +357,9 @@ det(const TinyMatrix<N, T>& A)
M(K, l) -= factor * M(I, l);
}
}
} else {
return 0;
}
}
for (size_t i = 0; i < N; ++i) {
......
......@@ -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")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment