Create invertmatrix.h

Named matrix-inversion routines less cryptically. "ur0" meant that the
upper-right triangle was all zero, because the matrix was a
well-formed box-vector matrix, per docs in vec.h. Added unit tests and
Doxygen. There's not actually a maths module yet, but it's better
than it used to be.

Change-Id: I4fe53acef07c6710c3a6ee49fc344d312694cf2c

diff --git a/src/gromacs/ewald/pme.cpp b/src/gromacs/ewald/pme.cpp
index 7c9a29ead3f51499860fe81482f27d704f86832e..279671c7b77f7b0491e565a16c2b53b958232c8b 100644 (file)
--- a/src/gromacs/ewald/pme.cpp
+++ b/src/gromacs/ewald/pme.cpp
@@ -86,6 +86,7 @@
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/gmxlib/nrnb.h"
 #include "gromacs/math/gmxcomplex.h"
+#include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vectypes.h"
@@ -959,7 +960,7 @@ int gmx_pme_do(struct gmx_pme_t *pme,
         atc->f = f;
     }
 
-    m_inv_ur0(box, pme->recipbox);
+    gmx::invertBoxMatrix(box, pme->recipbox);
     bFirst = TRUE;
 
     /* For simplicity, we construct the splines for all particles if

diff --git a/src/gromacs/gmxana/gmx_vanhove.cpp b/src/gromacs/gmxana/gmx_vanhove.cpp
index 5b040492e7770a13d57eebb4aa7faec0f61b54fd..5b0249c4b22404236e9d9310d60c8e1ccf2f628b 100644 (file)
--- a/src/gromacs/gmxana/gmx_vanhove.cpp
+++ b/src/gromacs/gmxana/gmx_vanhove.cpp
@@ -49,6 +49,7 @@
 #include "gromacs/gmxana/gmx_ana.h"
 #include "gromacs/gmxana/gstat.h"
 #include "gromacs/math/functions.h"
+#include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/topology/index.h"
@@ -296,7 +297,7 @@ int gmx_vanhove(int argc, char *argv[])
             fprintf(stderr, "\rProcessing frame %d", f);
         }
         /* Scale all the configuration to the average box */
-        m_inv_ur0(sbox[f], corr);
+        gmx::invertBoxMatrix(sbox[f], corr);
         mmul_ur0(avbox, corr, corr);
         for (i = 0; i < isize; i++)
         {

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index 4cdd35a256087df07759bcf3ff0a1f2f1f035157..a87b28be64b7c17c5b5003f23bc67a70300b0be3 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -68,6 +68,7 @@
 #include "gromacs/gmxpreprocess/vsite_parm.h"
 #include "gromacs/imd/imd.h"
 #include "gromacs/math/functions.h"
+#include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/calc_verletbuf.h"
@@ -840,7 +841,7 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
             clear_rvec(invbox[j]);
             invbox[j][j] = 1;
         }
-        m_inv_ur0(invbox, invbox);
+        gmx::invertBoxMatrix(invbox, invbox);
     }
 
     /* Copy the reference coordinates to mtop */

diff --git a/src/gromacs/math/invertmatrix.cpp b/src/gromacs/math/invertmatrix.cpp
index a597f10aa20407d2f249c2633e8f947738e63e94..87100100d0af611977aa940f99057c83b43212b6 100644 (file)
--- a/src/gromacs/math/invertmatrix.cpp
+++ b/src/gromacs/math/invertmatrix.cpp
@@ -32,14 +32,26 @@
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
+/*! \internal \file
+ * \brief
+ * Routines to invert 3x3 matrices
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_math
+ */
 #include "gmxpre.h"
 
+#include "invertmatrix.h"
+
 #include <cmath>
 
-#include "gromacs/math/vec.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
 
-void m_inv_ur0(const matrix src, matrix dest)
+namespace gmx
+{
+
+void invertBoxMatrix(const matrix src, matrix dest)
 {
     double tmp = src[XX][XX]*src[YY][YY]*src[ZZ][ZZ];
     if (std::fabs(tmp) <= 100*GMX_REAL_MIN)
@@ -59,7 +71,7 @@ void m_inv_ur0(const matrix src, matrix dest)
     dest[YY][ZZ] = 0.0;
 }
 
-void m_inv(const matrix src, matrix dest)
+void invertMatrix(const matrix src, matrix dest)
 {
     const real smallreal = (real)1.0e-24;
     const real largereal = (real)1.0e24;
@@ -72,6 +84,7 @@ void m_inv(const matrix src, matrix dest)
     {
         gmx_fatal(FARGS, "Can not invert matrix, determinant = %e", determinant);
     }
+    GMX_ASSERT(dest != src, "Cannot do in-place inversion of matrix");
 
     dest[XX][XX] = c*(src[YY][YY]*src[ZZ][ZZ]-src[ZZ][YY]*src[YY][ZZ]);
     dest[XX][YY] = -c*(src[XX][YY]*src[ZZ][ZZ]-src[ZZ][YY]*src[XX][ZZ]);
@@ -83,3 +96,5 @@ void m_inv(const matrix src, matrix dest)
     dest[ZZ][YY] = -c*(src[XX][XX]*src[ZZ][YY]-src[ZZ][XX]*src[XX][YY]);
     dest[ZZ][ZZ] = c*(src[XX][XX]*src[YY][YY]-src[YY][XX]*src[XX][YY]);
 }
+
+} // namespace

diff --git a/src/gromacs/math/invertmatrix.h b/src/gromacs/math/invertmatrix.h
new file mode 100644 (file)
index 0000000..6bd31fe
--- /dev/null
+++ b/src/gromacs/math/invertmatrix.h
@@ -0,0 +1,70 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \libinternal \file
+ * \brief
+ * Declares routines to invert 3x3 matrices
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_math
+ * \inlibraryapi
+ */
+#ifndef GMX_MATH_INVERTMATRIX_H
+#define GMX_MATH_INVERTMATRIX_H
+
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/basedefinitions.h"
+
+namespace gmx
+{
+
+/*! \brief Invert a simulation-box matrix in \c src, return in \c dest
+ *
+ * This routine assumes that src is a simulation-box matrix, i.e. has
+ * zeroes in the upper-right triangle. A fatal error occurs if the
+ * product of the leading diagonal is too small. The inversion can be
+ * done "in place", i.e \c src and \c dest can be the same matrix.
+ */
+void invertBoxMatrix(const matrix src, matrix dest);
+
+/*! \brief Invert a general 3x3 matrix in \c src, return in \c dest
+ *
+ * A fatal error occurs if the determinant is too small. \c src and
+ * \c dest cannot be the same matrix.
+ */
+void invertMatrix(const matrix src, matrix dest);
+
+} // namespace gmx
+
+#endif

diff --git a/src/gromacs/math/tests/CMakeLists.txt b/src/gromacs/math/tests/CMakeLists.txt
index 44d9ea463fb454ce3750ee3b6dd252c6132d1324..781248c0238ee51ce2f4af912feff667d5a1c2b5 100644 (file)
--- a/src/gromacs/math/tests/CMakeLists.txt
+++ b/src/gromacs/math/tests/CMakeLists.txt
@@ -34,5 +34,6 @@
 
 gmx_add_unit_test(MathUnitTests math-test
                   functions.cpp
+                  invertmatrix.cpp
                   vectypes.cpp
                   )

diff --git a/src/gromacs/math/tests/invertmatrix.cpp b/src/gromacs/math/tests/invertmatrix.cpp
new file mode 100644 (file)
index 0000000..5dadb3c
--- /dev/null
+++ b/src/gromacs/math/tests/invertmatrix.cpp
@@ -0,0 +1,145 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests matrix inversion routines
+ *
+ * \todo Test error conditions when they throw exceptions
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_math
+ */
+#include "gmxpre.h"
+
+#include "gromacs/math/invertmatrix.h"
+
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/math/vec.h"
+
+#include "testutils/testasserts.h"
+
+namespace
+{
+
+using gmx::invertMatrix;
+using gmx::invertBoxMatrix;
+using gmx::test::defaultRealTolerance;
+
+TEST(InvertMatrixTest, IdentityIsImpotent)
+{
+    matrix in = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+    matrix out;
+
+    invertMatrix(in, out);
+
+    EXPECT_REAL_EQ_TOL(out[XX][XX], in[XX][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[XX][YY], in[XX][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[XX][ZZ], in[XX][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[YY][XX], in[YY][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[YY][YY], in[YY][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[YY][ZZ], in[YY][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[ZZ][XX], in[ZZ][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[ZZ][YY], in[ZZ][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[ZZ][ZZ], in[ZZ][ZZ], defaultRealTolerance());
+}
+
+TEST(InvertMatrixTest, ComputesInverse)
+{
+    matrix in = {{1, 2, 3}, {-1, 2.5, 1}, {10, -2, 1.2}};
+    matrix out;
+    matrix expected = {{-0.12019230769230768,
+                        0.20192307692307693,
+                        0.13221153846153844},
+                       {-0.26923076923076916,
+                        0.69230769230769229,
+                        0.096153846153846145},
+                       {0.55288461538461531,
+                        -0.52884615384615374,
+                        -0.10817307692307691}};
+
+    invertMatrix(in, out);
+
+    EXPECT_REAL_EQ_TOL(out[XX][XX], expected[XX][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[XX][YY], expected[XX][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[XX][ZZ], expected[XX][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[YY][XX], expected[YY][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[YY][YY], expected[YY][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[YY][ZZ], expected[YY][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[ZZ][XX], expected[ZZ][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[ZZ][YY], expected[ZZ][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(out[ZZ][ZZ], expected[ZZ][ZZ], defaultRealTolerance());
+}
+
+TEST(InvertBoxMatrixTest, IdentityIsImpotent)
+{
+    matrix in = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+
+    invertBoxMatrix(in, in);
+
+    EXPECT_REAL_EQ_TOL(in[XX][XX], in[XX][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[XX][YY], in[XX][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[XX][ZZ], in[XX][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[YY][XX], in[YY][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[YY][YY], in[YY][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[YY][ZZ], in[YY][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[ZZ][XX], in[ZZ][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[ZZ][YY], in[ZZ][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(in[ZZ][ZZ], in[ZZ][ZZ], defaultRealTolerance());
+}
+
+TEST(InvertBoxMatrixTest, ComputesInverseInPlace)
+{
+    matrix in       = {{1, 0, 0}, {-1, 2.5, 0}, {10, -2, 1.2}};
+    matrix expected = {{1, 0, 0},
+                       {0.4, 0.4, 0},
+                       {-23.0/3.0, 2.0/3.0, 5.0/6.0}};
+
+    invertBoxMatrix(in, in);
+
+    EXPECT_REAL_EQ_TOL(expected[XX][XX], in[XX][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[XX][YY], in[XX][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[XX][ZZ], in[XX][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[YY][XX], in[YY][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[YY][YY], in[YY][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[YY][ZZ], in[YY][ZZ], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[ZZ][XX], in[ZZ][XX], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[ZZ][YY], in[ZZ][YY], defaultRealTolerance());
+    EXPECT_REAL_EQ_TOL(expected[ZZ][ZZ], in[ZZ][ZZ], defaultRealTolerance());
+}
+
+} // namespace

diff --git a/src/gromacs/math/vec.h b/src/gromacs/math/vec.h
index 160e1bb1a0e4dea16f59822b0cd9f35c5e9d622b..651876fe9865ef75ec482490fe3e0a821c54cbae 100644 (file)
--- a/src/gromacs/math/vec.h
+++ b/src/gromacs/math/vec.h
@@ -91,8 +91,6 @@
    void m_add(matrix a,matrix b,matrix dest)        dest = a + b
    void m_sub(matrix a,matrix b,matrix dest)        dest = a - b
    void msmul(matrix m1,real r1,matrix dest)        dest = r1 * m1
-   void m_inv_ur0(matrix src,matrix dest)           dest = src^-1
-   void m_inv(matrix src,matrix dest)            !  dest = src^-1
    void mvmul(matrix a,rvec src,rvec dest)       !  dest = a . src
    void mvmul_ur0(matrix a,rvec src,rvec dest)      dest = a . src
    void tmvmul_ur0(matrix a,rvec src,rvec dest)     dest = a* . src
@@ -570,10 +568,6 @@ static inline void msmul(const matrix m1, real r1, matrix dest)
     dest[ZZ][ZZ] = r1*m1[ZZ][ZZ];
 }
 
-/* Routines defined in invertmatrix.cpp */
-void m_inv_ur0(const matrix src, matrix dest);
-void m_inv(const matrix src, matrix dest);
-
 static inline void mvmul(const matrix a, const rvec src, rvec dest)
 {
     dest[XX] = a[XX][XX]*src[XX]+a[XX][YY]*src[YY]+a[XX][ZZ]*src[ZZ];

diff --git a/src/gromacs/mdlib/coupling.cpp b/src/gromacs/mdlib/coupling.cpp
index e057bbf9e02167dd8bff85e5cc517e806fe8ddf2..c2a3cb2b46780efd2c6275a28930cac512f6ebec 100644 (file)
--- a/src/gromacs/mdlib/coupling.cpp
+++ b/src/gromacs/mdlib/coupling.cpp
@@ -43,6 +43,7 @@
 #include "gromacs/domdec/domdec_struct.h"
 #include "gromacs/gmxlib/nrnb.h"
 #include "gromacs/math/functions.h"
+#include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vecdump.h"
@@ -375,7 +376,7 @@ void parrinellorahman_pcoupl(FILE *fplog, gmx_int64_t step,
 
     real   maxl;
 
-    m_inv_ur0(box, invbox);
+    gmx::invertBoxMatrix(box, invbox);
 
     if (!bFirstStep)
     {

diff --git a/src/gromacs/mdlib/csettle.cpp b/src/gromacs/mdlib/csettle.cpp
index f0e1aaadb46dd155779af2916c5dcfe1a8322fe6..8a2cebb018f439abc1553019a0b1f32462766b24 100644 (file)
--- a/src/gromacs/mdlib/csettle.cpp
+++ b/src/gromacs/mdlib/csettle.cpp
@@ -40,6 +40,7 @@
 #include <stdio.h>
 
 #include "gromacs/math/functions.h"
+#include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/constr.h"
 #include "gromacs/pbcutil/ishift.h"
@@ -83,7 +84,7 @@ static void init_proj_matrix(settleparam_t *p,
     p->imH = invmH;
     /* We normalize the inverse masses with imO for the matrix inversion.
      * so we can keep using masses of almost zero for frozen particles,
-     * without running out of the float range in m_inv.
+     * without running out of the float range in invertMatrix.
      */
     imOn = 1;
     imHn = p->imH/p->imO;
@@ -99,7 +100,7 @@ static void init_proj_matrix(settleparam_t *p,
     mat[2][0] = mat[0][2];
     mat[2][1] = mat[1][2];
 
-    m_inv(mat, p->invmat);
+    gmx::invertMatrix(mat, p->invmat);
 
     msmul(p->invmat, 1/p->imO, p->invmat);
 

diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp
index 3d6e1f2c940d67034f3162f1bbb5ea292d67ab79..9e4208c0a8489f644ea761a0e263d0fb33f2c557 100644 (file)
--- a/src/gromacs/mdlib/update.cpp
+++ b/src/gromacs/mdlib/update.cpp
@@ -50,6 +50,7 @@
 #include "gromacs/listed-forces/disre.h"
 #include "gromacs/listed-forces/orires.h"
 #include "gromacs/math/functions.h"
+#include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vecdump.h"
@@ -1292,7 +1293,7 @@ static void deform(gmx_update_t upd,
             }
         }
     }
-    m_inv_ur0(box, invbox);
+    gmx::invertBoxMatrix(box, invbox);
     copy_mat(bnew, box);
     mmul_ur0(box, invbox, mu);
 

diff --git a/src/gromacs/mdlib/vcm.cpp b/src/gromacs/mdlib/vcm.cpp
index ca797c47816b8372972d29f96978825846d60881..529ecfe21cec8a8a87c3bddcb39f5ef23eea4ce5 100644 (file)
--- a/src/gromacs/mdlib/vcm.cpp
+++ b/src/gromacs/mdlib/vcm.cpp
@@ -41,6 +41,7 @@
 
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/math/functions.h"
+#include "gromacs/math/invertmatrix.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vecdump.h"
 #include "gromacs/mdtypes/md_enums.h"
@@ -278,7 +279,7 @@ static void get_minv(tensor A, tensor B)
             tmp[m][n] *= fac;
         }
     }
-    m_inv(tmp, B);
+    gmx::invertMatrix(tmp, B);
     for (m = 0; (m < DIM); m++)
     {
         for (n = 0; (n < DIM); n++)