Wrote dimension dependent matrix inversion method and corresponding unit test....

Wrote dimension dependent matrix inversion method and corresponding unit test. Started making Gizmo dimension independent.

.gitignore

@@ -76,6 +76,7 @@ tests/testAdiabaticIndex
 tests/testRiemannExact
 tests/testRiemannTRRS
 tests/testRiemannHLLC
+tests/testMatrixInversion
 
 theory/latex/swift.pdf
 theory/kernel/kernels.pdf

src/dimension.h

@@ -33,6 +33,8 @@
 #include "inline.h"
 #include "vector.h"
 
+#include <math.h>
+
 /* First define some constants */
 #if defined(HYDRO_DIMENSION_3D)
 
@@ -114,6 +116,92 @@ __attribute__((always_inline)) INLINE static float pow_dimension_plus_one(
 #endif
 }
 
+/**
+ * @brief Inverts the given dimension by dimension matrix (in place)
+ *
+ * @param A A 3x3 matrix of which we want to invert the top left dxd part
+ */
+__attribute__((always_inline)) INLINE static void
+invert_dimension_by_dimension_matrix(float A[3][3]) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+  float detA, Ainv[3][3];
+
+  detA = A[0][0] * A[1][1] * A[2][2] + A[0][1] * A[1][2] * A[2][0] +
+         A[0][2] * A[1][0] * A[2][1] - A[0][2] * A[1][1] * A[2][0] -
+         A[0][1] * A[1][0] * A[2][2] - A[0][0] * A[1][2] * A[2][1];
+
+  if (detA && !isnan(detA)) {
+    Ainv[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) / detA;
+    Ainv[0][1] = (A[0][2] * A[2][1] - A[0][1] * A[2][2]) / detA;
+    Ainv[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) / detA;
+    Ainv[1][0] = (A[1][2] * A[2][0] - A[1][0] * A[2][2]) / detA;
+    Ainv[1][1] = (A[0][0] * A[2][2] - A[0][2] * A[2][0]) / detA;
+    Ainv[1][2] = (A[0][2] * A[1][0] - A[0][0] * A[1][2]) / detA;
+    Ainv[2][0] = (A[1][0] * A[2][1] - A[1][1] * A[2][0]) / detA;
+    Ainv[2][1] = (A[0][1] * A[2][0] - A[0][0] * A[2][1]) / detA;
+    Ainv[2][2] = (A[0][0] * A[1][1] - A[0][1] * A[1][0]) / detA;
+  } else {
+    Ainv[0][0] = 0.0f;
+    Ainv[0][1] = 0.0f;
+    Ainv[0][2] = 0.0f;
+    Ainv[1][0] = 0.0f;
+    Ainv[1][1] = 0.0f;
+    Ainv[1][2] = 0.0f;
+    Ainv[2][0] = 0.0f;
+    Ainv[2][1] = 0.0f;
+    Ainv[2][2] = 0.0f;
+  }
+
+  A[0][0] = Ainv[0][0];
+  A[0][1] = Ainv[0][1];
+  A[0][2] = Ainv[0][2];
+  A[1][0] = Ainv[1][0];
+  A[1][1] = Ainv[1][1];
+  A[1][2] = Ainv[1][2];
+  A[2][0] = Ainv[2][0];
+  A[2][1] = Ainv[2][1];
+  A[2][2] = Ainv[2][2];
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+  float detA, Ainv[2][2];
+
+  detA = A[0][0] * A[1][1] - A[0][1] * A[1][0];
+
+  if (detA && !isnan(detA)) {
+    Ainv[0][0] = A[1][1] / detA;
+    Ainv[0][1] = -A[0][1] / detA;
+    Ainv[1][0] = -A[1][0] / detA;
+    Ainv[1][1] = A[0][0] / detA;
+  } else {
+    Ainv[0][0] = 0.0f;
+    Ainv[0][1] = 0.0f;
+    Ainv[1][0] = 0.0f;
+    Ainv[1][1] = 0.0f;
+  }
+
+  A[0][0] = Ainv[0][0];
+  A[0][1] = Ainv[0][1];
+  A[1][0] = Ainv[1][0];
+  A[1][1] = Ainv[1][1];
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+  if (A[0][0] && !isnan(A[0][0])) {
+    A[0][0] = 1.0f / A[0][0];
+  } else {
+    A[0][0] = 0.0f;
+  }
+
+#else
+
+  error("The dimension is not defined !");
+
+#endif
+}
+
 /* ------------------------------------------------------------------------- */
 #ifdef WITH_VECTORIZATION
 

src/hydro/Gizmo/hydro.h

@@ -117,9 +117,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float h = p->h;
   const float ih = 1.0f / h;
   const float ih2 = ih * ih;
+  const float ihdim = pow_dimension(ih);
 
-  float detE, volume;
-  float E[3][3];
+  float volume;
   float m, momentum[3], energy;
 
 #ifndef THERMAL_ENERGY
@@ -129,48 +129,20 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* Final operation on the geometry. */
   /* we multiply with the smoothing kernel normalization ih3 and calculate the
    * volume */
-  volume = ih * ih2 * (p->geometry.volume + kernel_root);
+  volume = ihdim * (p->geometry.volume + kernel_root);
   p->geometry.volume = volume = 1. / volume;
   /* we multiply with the smoothing kernel normalization */
-  p->geometry.matrix_E[0][0] = E[0][0] = ih * ih2 * p->geometry.matrix_E[0][0];
-  p->geometry.matrix_E[0][1] = E[0][1] = ih * ih2 * p->geometry.matrix_E[0][1];
-  p->geometry.matrix_E[0][2] = E[0][2] = ih * ih2 * p->geometry.matrix_E[0][2];
-  p->geometry.matrix_E[1][0] = E[1][0] = ih * ih2 * p->geometry.matrix_E[1][0];
-  p->geometry.matrix_E[1][1] = E[1][1] = ih * ih2 * p->geometry.matrix_E[1][1];
-  p->geometry.matrix_E[1][2] = E[1][2] = ih * ih2 * p->geometry.matrix_E[1][2];
-  p->geometry.matrix_E[2][0] = E[2][0] = ih * ih2 * p->geometry.matrix_E[2][0];
-  p->geometry.matrix_E[2][1] = E[2][1] = ih * ih2 * p->geometry.matrix_E[2][1];
-  p->geometry.matrix_E[2][2] = E[2][2] = ih * ih2 * p->geometry.matrix_E[2][2];
-
-  /* invert the E-matrix */
-  /* code shamelessly stolen from the public version of GIZMO */
-  /* But since we should never invert a matrix, this code has to be replaced */
-  detE = E[0][0] * E[1][1] * E[2][2] + E[0][1] * E[1][2] * E[2][0] +
-         E[0][2] * E[1][0] * E[2][1] - E[0][2] * E[1][1] * E[2][0] -
-         E[0][1] * E[1][0] * E[2][2] - E[0][0] * E[1][2] * E[2][1];
-  /* check for zero determinant */
-  if ((detE != 0) && !isnan(detE)) {
-    p->geometry.matrix_E[0][0] = (E[1][1] * E[2][2] - E[1][2] * E[2][1]) / detE;
-    p->geometry.matrix_E[0][1] = (E[0][2] * E[2][1] - E[0][1] * E[2][2]) / detE;
-    p->geometry.matrix_E[0][2] = (E[0][1] * E[1][2] - E[0][2] * E[1][1]) / detE;
-    p->geometry.matrix_E[1][0] = (E[1][2] * E[2][0] - E[1][0] * E[2][2]) / detE;
-    p->geometry.matrix_E[1][1] = (E[0][0] * E[2][2] - E[0][2] * E[2][0]) / detE;
-    p->geometry.matrix_E[1][2] = (E[0][2] * E[1][0] - E[0][0] * E[1][2]) / detE;
-    p->geometry.matrix_E[2][0] = (E[1][0] * E[2][1] - E[1][1] * E[2][0]) / detE;
-    p->geometry.matrix_E[2][1] = (E[0][1] * E[2][0] - E[0][0] * E[2][1]) / detE;
-    p->geometry.matrix_E[2][2] = (E[0][0] * E[1][1] - E[0][1] * E[1][0]) / detE;
-  } else {
-    /* if the E-matrix is not well behaved, we cannot use it */
-    p->geometry.matrix_E[0][0] = 0.0f;
-    p->geometry.matrix_E[0][1] = 0.0f;
-    p->geometry.matrix_E[0][2] = 0.0f;
-    p->geometry.matrix_E[1][0] = 0.0f;
-    p->geometry.matrix_E[1][1] = 0.0f;
-    p->geometry.matrix_E[1][2] = 0.0f;
-    p->geometry.matrix_E[2][0] = 0.0f;
-    p->geometry.matrix_E[2][1] = 0.0f;
-    p->geometry.matrix_E[2][2] = 0.0f;
-  }
+  p->geometry.matrix_E[0][0] = ihdim * p->geometry.matrix_E[0][0];
+  p->geometry.matrix_E[0][1] = ihdim * p->geometry.matrix_E[0][1];
+  p->geometry.matrix_E[0][2] = ihdim * p->geometry.matrix_E[0][2];
+  p->geometry.matrix_E[1][0] = ihdim * p->geometry.matrix_E[1][0];
+  p->geometry.matrix_E[1][1] = ihdim * p->geometry.matrix_E[1][1];
+  p->geometry.matrix_E[1][2] = ihdim * p->geometry.matrix_E[1][2];
+  p->geometry.matrix_E[2][0] = ihdim * p->geometry.matrix_E[2][0];
+  p->geometry.matrix_E[2][1] = ihdim * p->geometry.matrix_E[2][1];
+  p->geometry.matrix_E[2][2] = ihdim * p->geometry.matrix_E[2][2];
+
+  invert_dimension_by_dimension_matrix(p->geometry.matrix_E);
 
   hydro_gradients_prepare_force_loop(p, ih2, volume);
 

src/kernel_hydro.h

@@ -252,7 +252,7 @@ __attribute__((always_inline)) INLINE static void kernel_deval(
 
   //#if kernel_ivals == 1
   ///* Only one branch in this case */
-  //const float *const coeffs = &kernel_coeffs[0];
+  // const float *const coeffs = &kernel_coeffs[0];
   //#else
   /* Pick the correct branch of the kernel */
   const int temp = (int)(x * kernel_ivals_f);
@@ -290,7 +290,7 @@ __attribute__((always_inline)) INLINE static void kernel_eval(
 
   //#if kernel_ivals == 1
   ///* Only one branch in this case */
-  //const float *const coeffs = &kernel_coeffs[0];
+  // const float *const coeffs = &kernel_coeffs[0];
   //#else
   /* Pick the correct branch of the kernel */
   const int temp = (int)(x * kernel_ivals_f);

tests/Makefile.am

@@ -24,14 +24,15 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS
 TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
         testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh  \
         testParser.sh testSPHStep test125cells.sh testKernelGrav testFFT \
-        testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC
+        testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
+        testMatrixInversion
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
 		 testSPHStep testPair test27cells test125cells testParser \
                  testKernel testKernelGrav testFFT testInteractions testMaths testSymmetry \
                  testAdiabaticIndex testRiemannExact testRiemannTRRS \
-                 testRiemannHLLC
+                 testRiemannHLLC testMatrixInversion
 
 # Sources for the individual programs
 testGreetings_SOURCES = testGreetings.c
@@ -72,6 +73,8 @@ testRiemannTRRS_SOURCES = testRiemannTRRS.c
 
 testRiemannHLLC_SOURCES = testRiemannHLLC.c
 
+testMatrixInversion_SOURCES = testMatrixInversion.c
+
 # Files necessary for distribution
 EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
 	     test27cells.sh test27cellsPerturbed.sh tolerance.dat testParser.sh \

tests/testMatrixInversion.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com).
+ *
+ * This program 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include "const.h"
+#include "dimension.h"
+#include "error.h"
+#include "tools.h"
+
+void setup_matrix(float A[3][3]) {
+  A[0][0] = random_uniform(-1.0, 1.0);
+  A[0][1] = random_uniform(-1.0, 1.0);
+  A[0][2] = random_uniform(-1.0, 1.0);
+  A[1][0] = random_uniform(-1.0, 1.0);
+  A[1][1] = random_uniform(-1.0, 1.0);
+  A[1][2] = random_uniform(-1.0, 1.0);
+  A[2][0] = random_uniform(-1.0, 1.0);
+  A[2][1] = random_uniform(-1.0, 1.0);
+  A[2][2] = random_uniform(-1.0, 1.0);
+}
+
+int is_unit_matrix(float A[3][3]) {
+  int check = 1;
+
+  check &= (fabsf(A[0][0] - 1.0f) < 1.e-6f);
+
+#if defined(HYDRO_DIMENSION_2D) && defined(HYDRO_DIMENSION_3D)
+  check &= (fabsf(A[0][1]) < 1.e-6f);
+  check &= (fabsf(A[1][0]) < 1.e-6f);
+  check &= (fabsf(A[1][1] - 1.0f) < 1.e-6f);
+#if defined(HYDRO_DIMENSION_3D)
+  check &= (fabsf(A[0][2]) < 1.e-6f);
+  check &= (fabsf(A[1][2]) < 1.e-6f);
+  check &= (fabsf(A[2][0]) < 1.e-6f);
+  check &= (fabsf(A[2][1]) < 1.e-6f);
+  check &= (fabsf(A[2][2] - 1.0f) < 1.e-6f);
+#endif  // 3D
+#endif  // 2D and 3D
+
+  return check;
+}
+
+void print_matrix(float A[3][3], const char* s) {
+  message("Matrix %s:", s);
+#if defined(HYDRO_DIMENSION_1D)
+  message("[%.3e]", A[0][0]);
+#elif defined(HYDRO_DIMENSION_2D)
+  message("[%.3e, %.3e]", A[0][0], A[0][1]);
+  message("[%.3e, %.3e]", A[1][0], A[1][1]);
+#elif defined(HYDRO_DIMENSION_3D)
+  message("[%.8e, %.8e, %.8e]", A[0][0], A[0][1], A[0][2]);
+  message("[%.8e, %.8e, %.8e]", A[1][0], A[1][1], A[1][2]);
+  message("[%.8e, %.8e, %.8e]", A[2][0], A[2][1], A[2][2]);
+#endif
+}
+
+void multiply_matrices(float A[3][3], float B[3][3], float C[3][3]) {
+#if defined(HYDRO_DIMENSION_1D)
+  C[0][0] = A[0][0] * B[0][0];
+#elif defined(HYDRO_DIMENSION_2D)
+  for (int i = 0; i < 2; ++i) {
+    for (int j = 0; j < 2; ++j) {
+      C[i][j] = 0.0f;
+      for (int k = 0; k < 2; ++k) {
+        C[i][j] += A[i][k] * B[k][j];
+      }
+    }
+  }
+#elif defined(HYDRO_DIMENSION_3D)
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      C[i][j] = 0.0f;
+      for (int k = 0; k < 3; ++k) {
+        C[i][j] += A[i][k] * B[k][j];
+      }
+    }
+  }
+#endif
+}
+
+int main() {
+
+  float A[3][3], B[3][3], C[3][3];
+  setup_matrix(A);
+
+  memcpy(B, A, 9 * sizeof(float));
+
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      if (A[i][j] != B[i][j]) {
+        error("Matrices not equal after copy!");
+      }
+    }
+  }
+
+  invert_dimension_by_dimension_matrix(A);
+
+  multiply_matrices(A, B, C);
+
+  if (!is_unit_matrix(C)) {
+    print_matrix(A, "A");
+    print_matrix(B, "B");
+    print_matrix(C, "C");
+    error("Inverted matrix is wrong!");
+  }
+
+  return 0;
+}