- Setting break-down threshold trough smvars.sh

- Fixed SM3 to escape infinite update-loop

Makefile

@@ -21,7 +21,7 @@ else
     $(error No valid compiler environment set in $$ENV. \
 	First run: $$ source smvars.sh {intel | llvm | gnu})
 endif
-CXXFLAGS = $(OPT) $(ARCH) $(DEBUG) -fPIC
+CXXFLAGS = $(OPT) $(ARCH) $(DEBUG) -fPIC $(THRESHOLD)
 FFLAGS = $(CXXFLAGS)
 H5CXX = h5c++
 H5CXXFLAGS = $(CXXFLAGS)

include/SM_Helpers.hpp

@@ -5,6 +5,13 @@
 #include <iostream>
 #include <string>
 
+// #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,

smvars.sh

@@ -1,6 +1,7 @@
 #!/usr/bin/env bash
-
+unset THRESHOLD
 ENV=$1
+THRESHOLD=$2
 
 ## Set Sherman-Morrison root dir
 PWD=$(pwd)
@@ -54,3 +55,9 @@ then
 	export PATH=$SMROOT/bin:$PATH
 	export SMVARS=true
 fi
+
+## If a threshold is provided, export compiler flag
+if [[ $# -gt 1 ]]
+then
+    export THRESHOLD="-DTHRESHOLD=$THRESHOLD"
+fi

src/SM_Helpers.cpp

@@ -1,5 +1,14 @@
 #include "SM_Helpers.hpp"
 
+// Set common break-down threshold
+double threshold() {
+  const double threshold = THRESHOLD;
+#ifdef DEBUG
+  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_MaponiA3.cpp

@@ -76,7 +76,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
               << "] = " << beta << std::endl;
     std::cout << std::endl;
 #endif
-    if (fabs(beta) < 1e-3) {
+    if (fabs(beta) < threshold()) {
       std::cerr << "Break-down occured." << std::endl;
     }
 

src/SM_MaponiA3S.cpp

@@ -4,7 +4,7 @@
 #include "SM_MaponiA3S.hpp"
 #include "SM_Helpers.hpp"
 
-#define DEBUG
+// #define DEBUG
 
 void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
                double *Updates, unsigned int *Updates_index) {
@@ -75,8 +75,21 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
               << "] = " << beta << std::endl;
     std::cout << std::endl;
 #endif
-    if (fabs(beta) < 1e-3) {
+    if (fabs(beta) < threshold()) {
-      std::cerr << "Break-down occured." << std::endl;
+      std::cerr << "Breakdown condition triggered at " << component
+                << std::endl;
+      
+    //   for (unsigned int i = 0; i < Dim; i++) {
+    //     later_updates[later * Dim + i] = Updates[l * Dim + i] / 2.0;
+    //     ylk[l][p[l + 1]][i] /= 2.0;
+    //   }
+    //   later_index[later] = Updates_index[l];
+    //   later++;
+
+    //   den = 1 + C[Updates_index[l] - 1];
+    // }
+    // double iden = 1 / den;
+
     }
 
 // Compute intermediate update to Slater_inv

src/SM_Standard.cpp

@@ -3,12 +3,6 @@
 #include "SM_Standard.hpp"
 #include "SM_Helpers.hpp"
 
-// Set common break-down threshold
-double threshold() {
-  double const threshold = 1e-3;
-  return threshold;
-}
-
 // Naïve Sherman Morrison
 void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
@@ -82,9 +76,9 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
                 << std::endl;
 
       // U_l = U_l / 2 (do the split)
-      for (unsigned int j = 0; j < Dim; j++) {
+      for (unsigned int i = 0; i < Dim; i++) {
-        later_updates[later * Dim + j] = Updates[l * Dim + j] / 2.0;
+        later_updates[later * Dim + i] = Updates[l * Dim + i] / 2.0;
-        C[j] /= 2.0;
+        C[i] /= 2.0;
       }
       later_index[later] = Updates_index[l];
       later++;
@@ -166,7 +160,13 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     l += 1;
   }
 
-  if (later > 0) {
+  // If all the updates have failed, exit early with an error
+  if (later == N_updates) {
+    std::cerr << "SM3 cannot invert this matrix" << std::endl;
+    return;
+  }
+  // If some have failed, make a recursive call
+  else if (later > 0) {
     SM3(Slater_inv, Dim, later, later_updates, later_index);
   }
 }