Added cofactor calculation

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),
+  TinyMatrix<2,T> A_cofactors_T(A(1,1), -A(1,0),
                                 -A(0,1), A(0,0));
-  return std::move(cofactor_T *= inv_determinent);
+  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