Problems with MaponiA3 when calling from C++ with HDF5 datasets is fixed. Bug was caused by uninitialised values of NEWed variables at several places. There are still a few remaining memory leaks. Explore with:

valgrind --track-origins=yes --leak-check=full --show-leak-kinds=all

Helpers.hpp

@@ -7,9 +7,9 @@ using namespace std;
 
 template<typename T>
 unsigned int getMaxIndex(T *vector, unsigned int size) {
-    unsigned int i;
-    unsigned int max = vector[0];
-    unsigned int maxi = 0;
+    unsigned int i = 0;
+    unsigned int maxi = i;
+    unsigned int max = vector[maxi];
     for (i = 1; i < size; i++) {
         if (vector[i] > max) {
             max = vector[i];
@@ -64,7 +64,7 @@ T *transpose(T *A, unsigned int M) {
 
 template<typename T>
 T *matMul(T *A, T *B, unsigned int M) {
-    T *C = new T[M*M];
+    T *C = new T[M*M] {0};
     for (unsigned int i = 0; i < M; i++) {
         for (unsigned int j = 0; j < M; j++) {
             for (unsigned int k = 0; k < M; k++) {
@@ -75,20 +75,6 @@ T *matMul(T *A, T *B, unsigned int M) {
     return C;
 }
 
-template<typename T>
-T *matMul2(T *A, T *B, unsigned int M) {
-    T *C = new T[M*M];
-    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];
-            }
-        }
-    }
-    showMatrix(C,M,"C");
-    return C;
-}
-
 template<typename T1, typename T2>
 T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
     T1 *C = new T1[M*M];

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_8169.dat", &
+    call Read_dataset("datasets/update_cycle_13.dat", &
                        cycle_id, &
                        dim, &
                        n_updates, &

SM_MaponiA3.cpp

@@ -4,15 +4,17 @@
 #include "SM_MaponiA3.hpp"
 #include "Helpers.hpp"
 
-void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-              double *Updates, unsigned int *Updates_index) {
+void MaponiA3(double *Slater_inv, unsigned int Dim,
+              unsigned int N_updates, double *Updates, 
+              unsigned int *Updates_index) {
 
   unsigned int k, l, lbar, i, j, tmp, component;
-  unsigned int *p = new unsigned int[N_updates + 1];
+  unsigned int *p = new unsigned int[N_updates + 1] {0};
   double alpha, beta;
-  double *breakdown = new double[N_updates + 1];
+  double *breakdown = new double[N_updates + 1] {0};
   double *Al = new double[Dim * Dim];
-  p[0] = 0;
+
+  // Populate update-order vector
   for (i = 0; i < N_updates; i++) {
     p[i + 1] = i + 1;
   }
@@ -43,8 +45,8 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
       breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]);
     }
     lbar = getMaxIndex(breakdown, N_updates + 1);
-    // Reset breakdown back to 0 for next round to avoid case where its
-    // first element is always the largest
+    // Reset breakdown back to 0 for next round to avoid case where
+    // its first element is always the largest
     for (i = 0; i < N_updates + 1; i++) {
       breakdown[i] = 0;
     }
@@ -60,14 +62,16 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     for (k = l + 1; k < N_updates + 1; k++) {
       alpha = ylk[l - 1][p[k]][component] / beta;
       for (i = 1; i < Dim + 1; i++) {
-        ylk[l][p[k]][i] = ylk[l - 1][p[k]][i] - alpha * ylk[l - 1][p[l]][i];
+        ylk[l][p[k]][i] = ylk[l - 1][p[k]][i]
+                        - alpha * ylk[l - 1][p[l]][i];
       }
     }
   }
 
   // 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.
+  // '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;
 
   // Construct A-inverse from A0-inverse and the ylk
@@ -77,13 +81,15 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     beta = 1 + ylk[l][p[k]][component];
     for (i = 0; i < Dim; i++) {
       for (j = 0; j < Dim; j++) {
-        Al[i * Dim + j] = (i == j) - (j == component - 1) * ylk[l][p[k]][i + 1] / beta;
+        Al[i*Dim + j] = (i == j) - (j == component-1)
+                      * 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[][]'
+  // 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];

tests/test.cpp

@@ -48,11 +48,11 @@ int test_cycle(H5File file, int cycle) {
 
   /* Test */
 #ifdef DEBUG
-  showMatrix(slater_matrix, dim, "Slater");
+  showMatrix(slater_matrix, dim, "OLD Slater");
 #endif
 
 #ifdef DEBUG
-  showMatrix(slater_inverse, dim, "Inverse");
+  showMatrix(slater_inverse, dim, "OLD Inverse");
 #endif
 
   for (j = 0; j < nupdates; j++) {
@@ -65,21 +65,22 @@ int test_cycle(H5File file, int cycle) {
   MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
 
 #ifdef DEBUG
-  showMatrix(slater_matrix, dim, "Slater");
+  showMatrix(slater_matrix, dim, "NEW Slater");
 #endif
 
 #ifdef DEBUG
-  showMatrix(slater_inverse, dim, "Inverse");
+  showMatrix(slater_inverse, dim, "NEW Inverse");
 #endif
 
-  double * res = matMul2(slater_matrix, slater_inverse, dim);
-  bool ok = is_identity(res, dim, 1.0e-8);
+  double * res = matMul(slater_matrix, slater_inverse, dim);
+  bool ok = is_identity(res, dim, 0.5e-4);
 
 #ifdef DEBUG
   showMatrix(res, dim, "Result");
 #endif
 
-  delete [] res, updates, col_update_index, slater_matrix, slater_inverse;
+  delete [] res, updates, col_update_index,
+            slater_matrix, slater_inverse;
 
   return ok;
 }
@@ -95,10 +96,12 @@ int main(int argc, char **argv) {
   bool ok = test_cycle(file, cycle);
 
   if (ok) {
-    std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl;
+    std::cerr << "ok -- cycle " << std::to_string(cycle)
+              << std::endl;
   }
   else {
-    std::cerr << "failed -- cycle " << std::to_string(cycle) << std::endl;
+    std::cerr << "failed -- cycle " << std::to_string(cycle)
+              << std::endl;
   }
   return ok;
 }