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

Began determinent code

(must check perfs for pivoting approach even for 3x3 matrices)
parent 9f0c321e
No related branches found
No related tags found
No related merge requests found
......@@ -82,6 +82,10 @@ public:
KOKKOS_INLINE_FUNCTION
friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
{
if constexpr(N==1) {
os << A(0,0);
}
else {
os << '[';
for (size_t i=0; i<N; ++i) {
os << '(' << A(i,0);
......@@ -91,6 +95,7 @@ public:
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
......@@ -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);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment