From b0bced2782587822a04565c8b7dc4fb9da544030 Mon Sep 17 00:00:00 2001
From: Kevin Boyd <kevin.boyd@uconn.edu>
Date: Mon, 1 Jul 2019 23:31:02 -0400
Subject: [PATCH] Add trace functionality to 3x3 matrices

Replacement for vec.h trace

Change-Id: I5c452cb36bab277ffcb3f2f9645d3d025a1805a8
---
 src/gromacs/math/matrix.h         | 11 ++++++++++-
 src/gromacs/math/tests/matrix.cpp |  6 ++++++
 2 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/src/gromacs/math/matrix.h b/src/gromacs/math/matrix.h
index a725e69ef6..ed225ad585 100644
--- 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 37141209f1..406022aa1d 100644
--- 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
-- 
2.22.0