From 983f87d504ab7da3770511903864e1dc17442aba Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <francois.coppens@.uvsq.fr>
Date: Mon, 1 Mar 2021 16:08:14 +0100
Subject: [PATCH] Rewrote matMul function so it doesn't declare memory inside.

---
 Makefile                       |  8 ++++----
 include/Helpers.hpp            | 13 +++++++++++++
 src/SM_MaponiA3.cpp            | 32 +++++++++++++-------------------
 tests/QMCChem_dataset_test.f90 |  2 +-
 4 files changed, 31 insertions(+), 24 deletions(-)

diff --git a/Makefile b/Makefile
index 12deb07..c4fc58b 100644
--- a/Makefile
+++ b/Makefile
@@ -1,8 +1,8 @@
 ## Compilers
 ARCH =
-H5CXX = h5c++
-CXX = clang++
-FC = flang
+H5CXX = h5c++ -std=gnu++11
+CXX = g++
+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)
diff --git a/include/Helpers.hpp b/include/Helpers.hpp
index 466478d..9f3426b 100644
--- a/include/Helpers.hpp
+++ b/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];
diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp
index c09f40f..dfdd093 100644
--- a/src/SM_MaponiA3.cpp
+++ b/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
-  // 'Slater_inv' gets reassigned by 'matMul(...)' in the next line,
-  // by creating a new pointer 'copy' that points to whereever
-  // 'Slater_inv' points to now.
-  double *copy = Slater_inv;
-
+  double *next = new double[Dim*Dim] {0};
+  double *last;
+  last =  Slater_inv;
   // Construct A-inverse from A0-inverse and the ylk
-  for (l = 0; l < N_updates; l++) {       // l    = 0, 1
-    k = l + 1;                            // k    = 1, 2
-    component = Updates_index[p[k] - 1];  // comp = 2, 4
+  for (l = 0; l < N_updates; l++) {
+    k = l + 1;
+    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
-  // in 'copy[][]'
-  for (i = 0; i < Dim; i++) {
-    for (j = 0; j < Dim; j++) {
-      copy[i * Dim + j] = Slater_inv[i * Dim + j];
-    }
+    matMul2(Al, last, next, Dim);
+    double *tmp = next;
+    next = last;
+    last = tmp;
   }
-
+  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;
 }
 
diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90
index e4c4986..508c0f3 100644
--- a/tests/QMCChem_dataset_test.f90
+++ b/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, &