diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp index 94bf1ad7d020913999f01e88af84e39e6dd5f97c..c7447be5faff25326816a4e9f5417887ac72b9b6 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 d68d94569eb35c1bbd943ed35772f3d2c736db8a..e00e02c508baedac31c3810c180f452f86e44738 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);