Add trace functionality to 3x3 matrices

Replacement for vec.h trace

Change-Id: I5c452cb36bab277ffcb3f2f9645d3d025a1805a8

diff --git a/src/gromacs/math/matrix.h b/src/gromacs/math/matrix.h
index a725e69ef63e0aaefe370c4964896088e9e47fed..ed225ad5854c9a65f7825ccf9c131d3c874853fd 100644 (file)
--- a/src/gromacs/math/matrix.h
+++ b/src/gromacs/math/matrix.h
@@ -61,7 +61,11 @@ using BasicMatrix3x3 = MultiDimArray<std::array<ElementType, 3*3>, extents <3, 3
 /*! \brief Three-by-three real number matrix.
  * \note will replace the C-style real[3][3] "matrix"
  */
-using Matrix3x3 = BasicMatrix3x3<real>;
+using Matrix3x3          = BasicMatrix3x3<real>;
+//! Convenience alias for a matrix view
+using Matrix3x3Span      = Matrix3x3::view_type;
+//! Convenience alias for a const matrix view
+using Matrix3x3ConstSpan = Matrix3x3::const_view_type;
 
 //! Determinant of a 3x3 matrix
 template <class ElementType>
@@ -72,6 +76,11 @@ ElementType determinant(const BasicMatrix3x3<ElementType> matrix)
               +matrix(2, 0)*(matrix(0, 1)*matrix(1, 2)-matrix(1, 1)*matrix(0, 2)));
 }
 
+//! Calculates the trace of a 3x3 matrix view
+constexpr real trace(Matrix3x3ConstSpan matrixView)
+{
+    return matrixView(0, 0) + matrixView(1, 1) + matrixView(2, 2);
+}
 } // namespace gmx
 
 #endif

diff --git a/src/gromacs/math/tests/matrix.cpp b/src/gromacs/math/tests/matrix.cpp
index 37141209f1fdf057c0df7a73314e7010fe621d95..406022aa1df7e56058fb0defe9b5d0944ae8f745 100644 (file)
--- a/src/gromacs/math/tests/matrix.cpp
+++ b/src/gromacs/math/tests/matrix.cpp
@@ -168,6 +168,12 @@ TEST_F(MatrixTest, determinantOfDiagonalMatrix)
     EXPECT_EQ(determinant(mat), 24);
 }
 
+TEST_F(MatrixTest, traceWorks)
+{
+    const Matrix3x3 mat = {{1.5, 9, 9, 9, 2.0, 9, 9, 9, 0.25}};
+    EXPECT_EQ(trace(mat), 3.75);
+}
+
 } // namespace
 
 } // namespace test