Changes in src/helpers.cpp [5034e1:99593f]

File:

• src/helpers.cpp (modified) (5 diffs)

src/helpers.cpp

@@ -6 +6 @@
 
 #include "helpers.hpp"
+
+#include "memoryusageobserver.hpp"
 
 /********************************************** helpful functions *********************************/
@@ -47 +49 @@
   while (*b < lower_bound)
     *b += step;
-};
-
-/** Flips two doubles.
- * \param *x pointer to first double
- * \param *y pointer to second double
- */
-void flip(double *x, double *y)
-{
-  double tmp;
-  tmp = *x;
-  *x = *y;
-  *y = tmp;
 };
 
@@ -143 +133 @@
  * \return allocated NDIM*NDIM array with the symmetric matrix
  */
-double * ReturnFullMatrixforSymmetric(double *symm)
+double * ReturnFullMatrixforSymmetric(const double * const symm)
 {
   double *matrix = Malloc<double>(NDIM * NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
@@ -158 +148 @@
 };
 
+/** Calculate the inverse of a 3x3 matrix.
+ * \param *matrix NDIM_NDIM array
+ */
+double * InverseMatrix( const double * const A)
+{
+  double *B = Malloc<double>(NDIM * NDIM, "Vector::InverseMatrix: *B");
+  double detA = RDET3(A);
+  double detAReci;
+
+  for (int i=0;i<NDIM*NDIM;++i)
+    B[i] = 0.;
+  // calculate the inverse B
+  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
+    detAReci = 1./detA;
+    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
+    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
+    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
+    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
+    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
+    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
+    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
+    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
+    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
+  }
+  return B;
+};
+
+/** hard-coded determinant of a 3x3 matrix.
+ * \param a[9] matrix
+ * \return \f$det(a)\f$
+ */
+double RDET3(const double a[NDIM*NDIM])
+{
+  return ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3]);
+};
+
+/** hard-coded determinant of a 2x2 matrix.
+ * \param a[4] matrix
+ * \return \f$det(a)\f$
+ */
+double RDET2(const double a[4])
+{
+  return ((a[0])*(a[3])-(a[1])*(a[2]));
+};
+
+/** hard-coded determinant of a 2x2 matrix.
+ * \param a0 (0,0) entry of matrix
+ * \param a1 (0,1) entry of matrix
+ * \param a2 (1,0) entry of matrix
+ * \param a3 (1,1) entry of matrix
+ * \return \f$det(a)\f$
+ */
+double RDET2(const double a0, const double a1, const double a2, const double a3)
+{
+  return ((a0)*(a3)-(a1)*(a2));
+};
+
 /** Comparison function for GSL heapsort on distances in two molecules.
  * \param *a
@@ -202 +249 @@
  * Frees all memory registered by the memory observer and calls exit(225) afterwards.
  */
-static void performCriticalExit() {
+void performCriticalExit() {
   map<void*, size_t> pointers = MemoryUsageObserver::getInstance()->getPointersToAllocatedMemory();
   for (map<void*, size_t>::iterator runner = pointers.begin(); runner != pointers.end(); runner++) {