Rewrote matMul function so it doesn't declare memory inside.

Makefile

@@ -1,8 +1,8 @@
 ## Compilers
 ARCH =
-H5CXX = h5c++
+H5CXX = h5c++ -std=gnu++11
-CXX = clang++
+CXX = g++
-FC = flang
+FC = gfortran
 
 ## Compiler flags
 H5CXXFLAGS = -O0 -g
@@ -55,7 +55,7 @@ $(OBJ_DIR)/%_h5.o: $(TST_DIR)/%_h5.cpp $(INC_DIR)/* | $(OBJ_DIR)
 
 ## Fortran modules
 $(OBJ_DIR)/%_mod.o: $(SRC_DIR)/%_mod.f90 | $(OBJ_DIR)
-	$(FC) $(FFLAGS) $(ARCH) -module $(OBJ_DIR)/ -c -o $@ $<
+	$(FC) $(FFLAGS) $(ARCH) -J$(OBJ_DIR)/ -c -o $@ $<
 
 ## Fortran objects
 $(OBJ_DIR)/%.o: $(TST_DIR)/%.f90 | $(OBJ_DIR)

include/Helpers.hpp

@@ -3,6 +3,7 @@
 #include <iostream>
 #include <cmath>
 #include <string>
+#include <cstring>
 using namespace std;
 
 template<typename T>
@@ -76,6 +77,18 @@ T *matMul(T *A, T *B, unsigned int M) {
     return C;
 }
 
+template<typename T>
+void matMul2(T *A, T *B, T *C, unsigned int M) {
+    memset(C,0,M*M*sizeof(T));
+    for (unsigned int i = 0; i < M; i++) {
+        for (unsigned int j = 0; j < M; j++) {
+            for (unsigned int k = 0; k < M; k++) {
+                C[i*M+j] += A[i*M+k] * B[k*M+j];
+            }
+        }
+    }
+}
+
 template<typename T1, typename T2>
 T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
     T1 *C = new T1[M*M];

src/SM_MaponiA3.cpp

@@ -68,16 +68,13 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     }
   }
 
-  // Keep the memory location of the passed array 'Slater_inv' before
+  double *next = new double[Dim*Dim] {0};
-  // 'Slater_inv' gets reassigned by 'matMul(...)' in the next line,
+  double *last;
-  // by creating a new pointer 'copy' that points to whereever
+  last =  Slater_inv;
-  // 'Slater_inv' points to now.
-  double *copy = Slater_inv;
-
   // Construct A-inverse from A0-inverse and the ylk
-  for (l = 0; l < N_updates; l++) {       // l    = 0, 1
+  for (l = 0; l < N_updates; l++) {
-    k = l + 1;                            // k    = 1, 2
+    k = l + 1;
-    component = Updates_index[p[k] - 1];  // comp = 2, 4
+    component = Updates_index[p[k] - 1];
     beta = 1 + ylk[l][p[k]][component];
     for (i = 0; i < Dim; i++) {
       for (j = 0; j < Dim; j++) {
@@ -85,24 +82,21 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
                       * ylk[l][p[k]][i + 1] / beta;
       }
     }
-    Slater_inv = matMul(Al, Slater_inv, Dim);
-  }
 
-  // Assign the new values of 'Slater_inv' to the old values
+    matMul2(Al, last, next, Dim);
-  // in 'copy[][]'
+    double *tmp = next;
-  for (i = 0; i < Dim; i++) {
+    next = last;
-    for (j = 0; j < Dim; j++) {
+    last = tmp;
-      copy[i * Dim + j] = Slater_inv[i * Dim + j];
-    }
   }
-
+  memcpy(Slater_inv, last, Dim*Dim*sizeof(double));
+ 
   for (l = 0; l < N_updates; l++) {
     for (k = 0; k < N_updates + 1; k++) {
       delete[] ylk[l][k];
     }
     delete[] ylk[l];
   }
-  delete[] Al;
+  delete[] Al, next;
   delete[] p, breakdown;
 }
 

tests/QMCChem_dataset_test.f90

@@ -9,7 +9,7 @@ program QMCChem_dataset_test
     real(c_double), dimension(:,:), allocatable :: Updates
     real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_t
 
-    call Read_dataset("datasets/update_cycle_13.dat", &
+    call Read_dataset("datasets/update_cycle_8169.dat", &
                        cycle_id, &
                        dim, &
                        n_updates, &