Cleanup and compiler flags.

Makefile

@@ -16,7 +16,7 @@ else ifeq ($(ENV),GNU)
 	FC = gfortran
 	# ARCH = -mavx
 	ARCH =
-	OPT = -O1
+	OPT = -O0
 	DEBUG = -g
 else
     $(error No valid compiler environment set in $$ENV. \

include/Helpers.hpp

@@ -10,14 +10,6 @@
 
 #include <cstdio>
 
-
-// #define DEBUG
-#ifndef THRESHOLD
-#define THRESHOLD 1e-3
-#endif
-
-double threshold();
-
 void Switch(unsigned int *p, unsigned int l, unsigned int lbar);
 
 void selectLargestDenominator(unsigned int l, unsigned int N_updates,

include/SM_Maponi.hpp

@@ -1,8 +1,8 @@
 // P. Maponi Algorithm 3
 void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-              double *Updates, unsigned int *Updates_index);
+              double *Updates, unsigned int *Updates_index, double breakdown);
 
 // P. Maponi Algorithm 3 with J. Slagel splitting
 // http://hdl.handle.net/10919/52966
 void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-               double *Updates, unsigned int *Updates_index);
+               double *Updates, unsigned int *Updates_index, double breakdown);

qmckl

@@ -1 +1 @@
-Subproject commit a5e58c80d70978c12d24ff47b2362c9145af325c
+Subproject commit 978e20ac428e2232384c5024f123bc93605da485

src/Helpers.cpp

@@ -1,14 +1,5 @@
 #include "Helpers.hpp"
 
-// Set common break-down threshold
-double threshold() {
-  const double threshold = THRESHOLD;
-#ifdef DEBUG2
-  std::cerr << "Break-down threshold set to: " << threshold << std::endl;
-#endif
-  return threshold;
-}
-
 void Switch(unsigned int *p, unsigned int l, unsigned int lbar) {
   unsigned int tmp = p[l + 1];
   p[l + 1] = p[lbar];

src/SM_Maponi.cpp

@@ -5,7 +5,7 @@
 #include "Helpers.hpp"
 
 void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-              double *Updates, unsigned int *Updates_index) {
+              double *Updates, unsigned int *Updates_index, double breakdown) {
   std::cerr << "Called MaponiA3 with " << N_updates << " updates" << std::endl;
 
   /*
@@ -74,7 +74,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
               << "] = " << beta << std::endl;
     std::cerr << std::endl;
 #endif
-    if (std::fabs(beta) < threshold()) {
+    if (std::fabs(beta) < breakdown) {
 #ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
@@ -135,7 +135,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 }
 
 void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-               double *Updates, unsigned int *Updates_index) {
+               double *Updates, unsigned int *Updates_index, double breakdown) {
   std::cerr << "Called MaponiA3S with " << N_updates << " updates" << std::endl;
 
   /*
@@ -208,7 +208,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
               << "] = " << beta << std::endl;
     std::cerr << std::endl;
 #endif
-    if (std::fabs(beta) < threshold()) {
+    if (std::fabs(beta) < breakdown) {
 #ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
@@ -275,22 +275,22 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   delete[] Al, next, p, ylk;
 
   if (later > 0) {
-    MaponiA3S(Slater_inv, Dim, later, later_updates, later_index);
+    MaponiA3S(Slater_inv, Dim, later, later_updates, later_index, breakdown);
   }
 }
 
 extern "C" {
 void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
                 unsigned int *N_updates, double **linUpdates,
-                unsigned int **Updates_index) {
+                unsigned int **Updates_index, double breakdown) {
-  MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
 }
 }
 
 extern "C" {
 void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
                  unsigned int *N_updates, double **linUpdates,
-                 unsigned int **Updates_index) {
+                 unsigned int **Updates_index, double breakdown) {
-  MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index, breakdown);
 }
 }

tests/qmckl_test_h5.cpp

@@ -53,6 +53,7 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
   for (i = 0; i < nupdates; i++) {
 	  col_update_index[i] = temp[i];
   }
+  delete[] temp;
 
   double *updates = new double[nupdates * dim];
   read_double(file, group + "/updates", updates);
@@ -170,69 +171,10 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
   }
 #endif // PERF
   delete[] u, col_update_index;
+  rc = qmckl_context_destroy(context);
 
-  showMatrix(slater_matrix, dim, "Slater Matrix");
-  showMatrix(slater_inverse, dim, "Slater Inverse");
   double *res = new double[dim * dim]{0};
-  {
+  matMul2(slater_matrix, slater_inverse, res, dim_32, dim_32, dim_32);
-  for (unsigned int i = 0; i < dim; i++) {
-    for (unsigned int j = 0; j < dim; j++) {
-      for (unsigned int k = 0; k < dim; k++) {
-        res[i * dim + j] += slater_matrix[i * dim + k] * slater_inverse[k * dim + j];
-      }
-    }
-  }
-  }
-
-  //matMul2(slater_matrix, slater_inverse, res, dim_32, dim_32, dim_32);
-  //
-  //
-  for (unsigned int i = 0; i < dim; i++) {
-    printf("[");
-    for (unsigned int j = 0; j < dim; j++) {
-      if (slater_matrix[i * dim + j] >= 0) {
-        printf("  %17.10e,", slater_matrix[i * dim + j]);
-      } else {
-        printf(" %17.10e,", slater_matrix[i * dim + j]);
-      }
-    }
-    printf(" ],\n");
-  }
-  printf("\n\n");
-  //
-  //
-  //
-  //
-  for (unsigned int i = 0; i < dim; i++) {
-    printf("[");
-    for (unsigned int j = 0; j < dim; j++) {
-      if (slater_inverse[i * dim + j] >= 0) {
-        printf("  %17.10e,", slater_inverse[i * dim + j]);
-      } else {
-        printf(" %17.10e,", slater_inverse[i * dim + j]);
-      }
-    }
-    printf(" ],\n");
-  }
-  printf("\n\n");
-  //
-  //
-  //
-  //
-  for (unsigned int i = 0; i < dim; i++) {
-    printf("[");
-    for (unsigned int j = 0; j < dim; j++) {
-      if (res[i * dim + j] >= 0) {
-        printf("  %17.10e,", res[i * dim + j]);
-      } else {
-        printf(" %17.10e,", res[i * dim + j]);
-      }
-    }
-    printf(" ],\n");
-  }
-  printf("\n\n");
-  //
-  //
   bool ok = is_identity(res, dim, tolerance);
   double res_max = residual_max(res, dim);
   double res2 = residual_frobenius2(res, dim);

tests/test_h5.cpp

@@ -35,7 +35,6 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
 
   unsigned int dim, nupdates, col, i, j;
 
-
   read_int(file, group + "/slater_matrix_dim", &dim);
   read_int(file, group + "/nupdates", &nupdates);
 
@@ -175,68 +174,8 @@ int test_cycle(H5::H5File file, int cycle, std::string version, double breakdown
 #endif // PERF
   delete[] u, col_update_index;
 
-  showMatrix(slater_matrix, dim, "Slater Matrix");
-  showMatrix(slater_inverse, dim, "Slater Inverse");
   double *res = new double[dim * dim]{0};
-  {
+  matMul2(slater_matrix, slater_inverse, res, dim, dim, dim);
-  for (unsigned int i = 0; i < dim; i++) {
-    for (unsigned int j = 0; j < dim; j++) {
-      for (unsigned int k = 0; k < dim; k++) {
-        res[i * dim + j] += slater_matrix[i * dim + k] * slater_inverse[k * dim + j];
-      }
-    }
-  }
-  }
-
-  //matMul2(slater_matrix, slater_inverse, res, dim, dim, dim);
-  //
-  //
-  for (unsigned int i = 0; i < dim; i++) {
-    printf("[");
-    for (unsigned int j = 0; j < dim; j++) {
-      if (slater_matrix[i * dim + j] >= 0) {
-        printf("  %17.10e,", slater_matrix[i * dim + j]);
-      } else {
-        printf(" %17.10e,", slater_matrix[i * dim + j]);
-      }
-    }
-    printf(" ],\n");
-  }
-  printf("\n\n");
-  //
-  //
-  //
-  //
-  for (unsigned int i = 0; i < dim; i++) {
-    printf("[");
-    for (unsigned int j = 0; j < dim; j++) {
-      if (slater_inverse[i * dim + j] >= 0) {
-        printf("  %17.10e,", slater_inverse[i * dim + j]);
-      } else {
-        printf(" %17.10e,", slater_inverse[i * dim + j]);
-      }
-    }
-    printf(" ],\n");
-  }
-  printf("\n\n");
-  //
-  //
-  //
-  //
-  for (unsigned int i = 0; i < dim; i++) {
-    printf("[");
-    for (unsigned int j = 0; j < dim; j++) {
-      if (res[i * dim + j] >= 0) {
-        printf("  %17.10e,", res[i * dim + j]);
-      } else {
-        printf(" %17.10e,", res[i * dim + j]);
-      }
-    }
-    printf(" ],\n");
-  }
-  printf("\n\n");
-  //
-  //
   bool ok = is_identity(res, dim, tolerance);
   double res_max = residual_max(res, dim);
   double res2 = residual_frobenius2(res, dim);