Merge pull request #13 from fmgjcoppens/loop-merge

Loops merged

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,22 +3,9 @@
 #include <iostream>
 #include <cmath>
 #include <string>
+#include <cstring>
 using namespace std;
 
-template<typename T>
-unsigned int getMaxIndex(T *vector, unsigned int size) {
-    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];
-            maxi = i;
-        }
-    }
-    return maxi;
-}
-
 template<typename T>
 void showScalar(T scalar, string name) {
     cout << name << " = " << scalar << endl << endl;
@@ -50,6 +37,7 @@ void showMatrix(T *matrix, unsigned int M, string name) {
         cout << " ]," << endl;
     }
     cout << "]" << endl;
+    cout << endl;
 }
 
 template<typename T>
@@ -64,8 +52,8 @@ T *transpose(T *A, unsigned int M) {
 }
 
 template<typename T>
-T *matMul(T *A, T *B, unsigned int M) {
+void matMul(T *A, T *B, T *C, unsigned int M) {
-    T *C = new T[M*M] {0};
+    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++) {
@@ -73,7 +61,6 @@ T *matMul(T *A, T *B, unsigned int M) {
             }
         }
     }
-    return C;
 }
 
 template<typename T1, typename T2>

src/SM_MaponiA3.cpp

@@ -4,22 +4,53 @@
 #include "SM_MaponiA3.hpp"
 #include "Helpers.hpp"
 
+void selectBestUpdate(unsigned int l, unsigned int N_updates, 
+                      unsigned int *Updates_index, unsigned int *p, 
+                      double ***ylk) {
+  unsigned int lbar = l+1, max =0;
+  unsigned int index = 0, component = 0;
+  unsigned int tmp = 0;
+  double breakdown = 0;
+  for (unsigned int j = lbar; j < N_updates + 1; j++) {
+    index = p[j];
+    component = Updates_index[index - 1];
+    breakdown = abs(1 + ylk[l][index][component]);
+    #ifdef DEBUG
+    cout << "Inside selectBestUpdate()" << endl;
+    cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "])" << endl;
+    cout << endl;
+    #endif
+    if (breakdown > max) {
+      max = breakdown;
+      lbar = j;
+    }
+  }
+  tmp = p[l+1];
+  p[l+1] = p[lbar];
+  p[lbar] = tmp; 
+}
+
 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;
+  /*
+  DECLARE AND INITIALISE ARRAYS
+  */
+
+  unsigned int k, l, i, j, component;
   unsigned int *p = new unsigned int[N_updates + 1] {0};
   double alpha, beta;
-  double *breakdown = new double[N_updates + 1] {0};
   double *Al = new double[Dim * Dim];
+  double *next = new double[Dim*Dim] {0};
+  double *last = Slater_inv, *tmp;
 
   // Populate update-order vector
   for (i = 0; i < N_updates; i++) {
     p[i + 1] = i + 1;
   }
 
-  // Declare auxiliary solution matrix ylk
+  // Declare auxiliary solution matrix ylk[N_updates][N_updates+1][Dim]
   double ***ylk = new double **[N_updates];
   for (l = 0; l < N_updates; l++) {
     ylk[l] = new double *[N_updates + 1];
@@ -28,8 +59,19 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     }
   }
 
-  // Calculate the y0k
+
+
+  /*
+  START ALGORITHM
+  */
+
+  // Calculate the {y_{0,k}}
   for (k = 1; k < N_updates + 1; k++) {
+    #ifdef DEBUG
+    cout << "Compute y0k: " << endl;
+    cout << "ylk[0][" << k << "][:]" << endl;
+    cout << endl;
+    #endif
     for (i = 1; i < Dim + 1; i++) {
       for (j = 1; j < Dim + 1; j++) {
         ylk[0][k][i] += Slater_inv[(i-1)*Dim + (j-1)]
@@ -38,63 +80,73 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     }
   }
 
-  // Calculate all the ylk from the y0k
+  // Calculate the {y_{l,k}} from the {y_{0,k}}
-  for (l = 1; l < N_updates; l++) {
+  for (l = 0; l < N_updates; l++) {
-    for (j = l; j < N_updates + 1; j++) {
+    #ifdef DEBUG
-      component = Updates_index[p[j] - 1];
+    cout << "In outer compute-ylk-loop: l = " << l << endl;
-      breakdown[j] = abs(1 + ylk[l - 1][p[j]][component]);
+    cout << endl;
-    }
+    #endif
-    lbar = getMaxIndex(breakdown, N_updates + 1);
+
-    // Reset breakdown back to 0 for next round to avoid case where
+    // For given l select intermediate update with largest break-down val
-    // its first element is always the largest
+    selectBestUpdate(l, N_updates, Updates_index, p, ylk);
-    for (i = 0; i < N_updates + 1; i++) {
+
-      breakdown[i] = 0;
+    // Select component and comp. bd-condition.
-    }
+    component = Updates_index[p[l+1] - 1];
-    tmp = p[l];
+    beta = 1 + ylk[l][p[l+1]][component];
-    p[l] = p[lbar];
+    #ifdef DEBUG
-    p[lbar] = tmp;
+    cout << "p[l+1] = " << p[l+1] << endl;
-    component = Updates_index[p[l] - 1];
+    cout << "component = " << component << endl;
-    beta = 1 + ylk[l - 1][p[l]][component];
+    cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << endl;
+    cout << endl;
+    #endif
     if (fabs(beta) < 1e-6) {
       cout << "Break-down occured. Exiting..." << endl;
       exit(1);
     }
-    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];
-      }
-    }
-  }
 
-  // Keep the memory location of the passed array 'Slater_inv' before
+    // Compute intermediate update to Slater_inv
-  // 'Slater_inv' gets reassigned by 'matMul(...)' in the next line,
+    #ifdef DEBUG
-  // by creating a new pointer 'copy' that points to whereever
+    cout << "Compute intermediate update to Slater_inv" << endl;
-  // 'Slater_inv' points to now.
+    cout << "component = " << component << endl;
-  double *copy = Slater_inv;
+    cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << endl;
-
+    cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l+1] << "][:]" << endl;
-  // Construct A-inverse from A0-inverse and the ylk
+    cout << endl;
-  for (l = 0; l < N_updates; l++) {       // l    = 0, 1
+    #endif
-    k = l + 1;                            // k    = 1, 2
-    component = Updates_index[p[k] - 1];  // comp = 2, 4
-    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;
+                      * ylk[l][p[l+1]][i + 1] / beta;
       }
     }
-    Slater_inv = matMul(Al, Slater_inv, Dim);
+    matMul(Al, last, next, Dim);
-  }
+    tmp = next;
+    next = last;
+    last = tmp;
+    #ifdef DEBUG
+    showMatrix(last, Dim, "last");
+    #endif
 
-  // Assign the new values of 'Slater_inv' to the old values
+    // For given l != 0 compute the next {y_{l,k}}
-  // in 'copy[][]'
+    for (k = l+2; k < N_updates+1; k++) {
-  for (i = 0; i < Dim; i++) {
+      alpha = ylk[l][p[k]][component] / beta;
-    for (j = 0; j < Dim; j++) {
+      #ifdef DEBUG
-      copy[i * Dim + j] = Slater_inv[i * Dim + j];
+      cout << "Inside k-loop: k = " << k << endl;
+      cout << "ylk[" << l+1 << "][" << p[k] << "][:]" << endl;
+      cout << endl;
+      #endif
+      for (i = 1; i < Dim + 1; i++) {
+        ylk[l+1][p[k]][i] = ylk[l][p[k]][i]
+                        - alpha * ylk[l][p[l+1]][i];
+      }
     }
   }
+  memcpy(Slater_inv, last, Dim*Dim*sizeof(double));
+
+
+
+  /*
+  CLEANUP MEMORY
+  */
   
   for (l = 0; l < N_updates; l++) {
     for (k = 0; k < N_updates + 1; k++) {
@@ -102,8 +154,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     }
     delete[] ylk[l];
   }
-  delete[] Al;
+  delete[] Al, next, p;
-  delete[] p, breakdown;
 }
 
 extern "C" {

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, &