Add function for calculating 3x3 matrix determinant

C++ matrix compatible version of the det() function in math/vec.h

Change-Id: Iaa2749dc8059fb6c3c112c4798f8e8cc8eae993c

diff --git a/src/gromacs/math/matrix.h b/src/gromacs/math/matrix.h
index 9d3a9dfcf640bb8194292a801348a7798d36297f..a725e69ef63e0aaefe370c4964896088e9e47fed 100644 (file)
--- a/src/gromacs/math/matrix.h
+++ b/src/gromacs/math/matrix.h
@@ -34,7 +34,7 @@
  */
 /*! \libinternal
  * \file
- * \brief Declares special case of 3x3 matrix frequently used.
+ * \brief Declares special case of 3x3 matrix frequently used, and associated functions.
  *
  * \author Christian Blau <cblau@gwdg.de>
  * \ingroup module_math
@@ -63,6 +63,15 @@ using BasicMatrix3x3 = MultiDimArray<std::array<ElementType, 3*3>, extents <3, 3
  */
 using Matrix3x3 = BasicMatrix3x3<real>;
 
+//! Determinant of a 3x3 matrix
+template <class ElementType>
+ElementType determinant(const BasicMatrix3x3<ElementType> matrix)
+{
+    return (  matrix(0, 0)*(matrix(1, 1)*matrix(2, 2)-matrix(2, 1)*matrix(1, 2))
+              -matrix(1, 0)*(matrix(0, 1)*matrix(2, 2)-matrix(2, 1)*matrix(0, 2))
+              +matrix(2, 0)*(matrix(0, 1)*matrix(1, 2)-matrix(1, 1)*matrix(0, 2)));
+}
+
 } // namespace gmx
 
 #endif

diff --git a/src/gromacs/math/tests/matrix.cpp b/src/gromacs/math/tests/matrix.cpp
index d246d0e0ab6c884d5c55472e89ee97e6409217f2..37141209f1fdf057c0df7a73314e7010fe621d95 100644 (file)
--- a/src/gromacs/math/tests/matrix.cpp
+++ b/src/gromacs/math/tests/matrix.cpp
@@ -142,6 +142,32 @@ TEST_F(MatrixTest, staticMultiDimArrayExtent)
     EXPECT_EQ(matrix_.extent(1), 3);
 }
 
+TEST_F(MatrixTest, determinantWorksForInt)
+{
+    const BasicMatrix3x3<int> mat = {{6, 4, 2, 1, -2, 8, 1, 5, 7}};
+    EXPECT_EQ(determinant(mat), -306);
+}
+
+TEST_F(MatrixTest, determinantWorksForFloat)
+{
+    const BasicMatrix3x3<float> mat = {{1.0, 2.0, 3.0,
+                                        0.0, 1.0, 4.0,
+                                        5.0, 6.0, 0.0}};
+    EXPECT_EQ(determinant(mat), 1);
+}
+
+TEST_F(MatrixTest, noninvertableDeterminantIsZero)
+{
+    const BasicMatrix3x3<int> mat = {{1, 0, 0, 0, 1, 0, 0, 0, 0}};
+    EXPECT_EQ(determinant(mat), 0);
+}
+
+TEST_F(MatrixTest, determinantOfDiagonalMatrix)
+{
+    const BasicMatrix3x3<int> mat = {{2, 0, 0, 0, 3, 0, 0, 0, 4}};
+    EXPECT_EQ(determinant(mat), 24);
+}
+
 } // namespace
 
 } // namespace test