diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp index 96ba7680fb125726f5c8b015f5adfe117a9ae807..1bb3a3a4488819923b81faea28c5d9ff9f649205 100644 --- a/src/algebra/TinyMatrix.hpp +++ b/src/algebra/TinyMatrix.hpp @@ -409,6 +409,17 @@ TinyMatrix<1,T> inverse(const TinyMatrix<1,T>& A) return std::move(A_1); } +template <size_t N, typename T> +KOKKOS_INLINE_FUNCTION +T cofactor(const TinyMatrix<N,T>& A, + const size_t& i, + const size_t& j) +{ + static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non arithmetic types"); + const T sign = ((i+j)%2) ? 1 : -1; + return sign * det(getMinor(A, i, j)); +} + template <typename T> KOKKOS_INLINE_FUNCTION TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A) @@ -419,11 +430,12 @@ TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A) const T determinent = det(A); const T inv_determinent = 1./determinent; - TinyMatrix<2,T> cofactor_T(A(1,1), -A(1,0), - -A(0,1), A(0,0)); - return std::move(cofactor_T *= inv_determinent); + TinyMatrix<2,T> A_cofactors_T(A(1,1), -A(1,0), + -A(0,1), A(0,0)); + return std::move(A_cofactors_T *= inv_determinent); } + template <typename T> KOKKOS_INLINE_FUNCTION TinyMatrix<3,T> inverse(const TinyMatrix<3,T>& A) @@ -432,13 +444,12 @@ TinyMatrix<3,T> inverse(const TinyMatrix<3,T>& A) static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only"); const T determinent = det(A); - const T inv_determinent = 1./determinent; - TinyMatrix<3,T> cofactor_T(det(getMinor(A,0,0)), det(getMinor(A,1,0)), det(getMinor(A,2,0)), - det(getMinor(A,0,1)), det(getMinor(A,1,1)), det(getMinor(A,2,1)), - det(getMinor(A,0,2)), det(getMinor(A,1,2)), det(getMinor(A,2,2))); + TinyMatrix<3,T> A_cofactors_T(cofactor(A,0,0), cofactor(A,1,0), cofactor(A,2,0), + cofactor(A,0,1), cofactor(A,1,1), cofactor(A,2,1), + cofactor(A,0,2), cofactor(A,1,2), cofactor(A,2,2)); - return std::move(cofactor_T *= inv_determinent); + return std::move(A_cofactors_T *= 1./determinent); } #endif // TINYMATRIX_HPP