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

@@ -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;
     }
   }