From: Kevin Boyd <kevin.boyd@uconn.edu>
Date: Fri, 21 Jun 2019 01:25:42 +0000 (-0400)
Subject: Add function for calculating 3x3 matrix determinant
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=f034e66c5d50f3284f7ab6490174ea616579a528;p=alexxy%2Fgromacs.git

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 9d3a9dfcf6..a725e69ef6 100644
--- 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 d246d0e0ab..37141209f1 100644
--- 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