From 4f57cfc742c5d69b831b3ead1075ccbeb6e7f0b8 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 27 Apr 2018 01:18:08 +0200
Subject: [PATCH] Began determinent code
(must check perfs for pivoting approach even for 3x3 matrices)
---
src/algebra/TinyMatrix.hpp | 62 +++++++++++++++++++++++++++++++++-----
src/main.cpp | 12 ++++++++
2 files changed, 67 insertions(+), 7 deletions(-)
diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 94bf1ad7d..c7447be5f 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -82,15 +82,20 @@ public:
KOKKOS_INLINE_FUNCTION
friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
{
- os << '[';
- for (size_t i=0; i<N; ++i) {
- os << '(' << A(i,0);
- for (size_t j=1; j<N; ++j) {
- os << ',' << A(i,j);
+ if constexpr(N==1) {
+ os << A(0,0);
+ }
+ else {
+ os << '[';
+ for (size_t i=0; i<N; ++i) {
+ os << '(' << A(i,0);
+ for (size_t j=1; j<N; ++j) {
+ os << ',' << A(i,j);
+ }
+ os << ')';
}
- os << ')';
+ os << ']';
}
- os << ']';
return os;
}
@@ -246,4 +251,47 @@ TinyMatrix<N,T> tensorProduct(const TinyVector<N,T>& x,
return std::move(A);
}
+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;
+}
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION
+T det(const TinyMatrix<1,T>& A)
+{
+ static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+ TinyMatrix<1,T> M=A;
+
+ return A(0,0);
+}
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION
+T det(const TinyMatrix<2,T>& A)
+{
+ static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+ TinyMatrix<2,T> M=A;
+
+ return A(0,0)*A(1,1)-A(1,0)*A(0,1);
+}
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION
+T det(const TinyMatrix<3,T>& A)
+{
+ static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
+ 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));
+}
+
#endif // TINYMATRIX_HPP
diff --git a/src/main.cpp b/src/main.cpp
index d68d94569..e00e02c50 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -120,6 +120,18 @@ int main(int argc, char *argv[])
// }
+
+ TinyMatrix<2,double> A( 2,1,
+ -1,3);
+ std::cout << "det" << A << "=" << det(A) << '\n';
+
+ TinyMatrix<3,double> B(6,1,1,
+ 4,-2,5,
+ 2,8,7);
+ std::cout << "det" << B << "=" << det(B) << '\n';
+ std::exit(0);
+
+
if (filename != "") {
std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
GmshReader gmsh_reader(filename);
--
GitLab