From 89c7b8cfe822a184572d97e6e2e509088e98e1cf Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 27 Apr 2018 15:27:23 +0200
Subject: [PATCH] computes determinent using Gauss pivoting for matrices of Rdd
for d>3
---
src/algebra/TinyMatrix.hpp | 45 +++++++++++++++++++++++++++++++++-----
1 file changed, 39 insertions(+), 6 deletions(-)
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index c7447be5f..05cbd6813 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -255,10 +255,43 @@ template <size_t N, typename T>
KOKKOS_INLINE_FUNCTION
T det(const TinyMatrix<N,T>& A)
{
- TinyMatrix<N,T> M=A;
-#warning must code a Gauss pivoting approach and compare perfs with Camers for 3x3 matrices
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
- return 0;
+ TinyMatrix<N,T> M = A;
+
+ TinyVector<N, int> index;
+ for (int i=0; i<N; ++i) index[i]=i;
+
+ T determinent = 1;
+ for (int i=0; i<N; ++i) {
+ for (int j=i; j<N; ++j) {
+ int l = j;
+ const int J = index[j];
+ for (int k=j+1; k<N; ++k) {
+ if (std::abs(M(index[k],i)) > std::abs(M(J,i))) {
+ l=k;
+ }
+ }
+ if (l != j) {
+ std::swap(index[l], index[j]);
+ determinent *= -1;
+ }
+ }
+ const int I = index[i];
+ if (M(I,i)==0) return 0;
+ double inv_Mii = 1./M(I,i);
+ for (int k=i+1; k<N; ++k) {
+ const int K = index[k];
+ const T factor = M(K,i)*inv_Mii;
+ for (int l=i+1; l<N; ++l) {
+ M(K,l) -= factor*M(I,l);
+ }
+ }
+ }
+
+ for (int i=0; i<N; ++i) {
+ determinent *= M(index[i],i);
+ }
+ return determinent;
}
template <typename T>
@@ -289,9 +322,9 @@ T det(const TinyMatrix<3,T>& A)
TinyMatrix<3,T> M=A;
return
- A(0,0)*(A(1,1)*A(2,2)-A(1,2)*A(2,1))
- -A(1,0)*(A(0,1)*A(2,2)-A(0,2)*A(2,1))
- +A(2,0)*(A(1,2)*A(0,1)-A(1,1)*A(0,2));
+ A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
+ -A(1,0)*(A(0,1)*A(2,2)-A(2,1)*A(0,2))
+ +A(2,0)*(A(0,1)*A(1,2)-A(1,1)*A(0,2));
}
#endif // TINYMATRIX_HPP
--
GitLab