Add a standard Sherman-Morisson implementation

- It can serve as a baseline reference

- It can serve as a starting point for including the pivot
  and splitting techniques from Maponi and Slaggel without the full
  complexity of the MaponiA3 algorithm

Makefile

@@ -1,16 +1,16 @@
 ## Compilers
-ARCH = -xCORE-AVX2
+ARCH =
 H5CXX = h5c++
-CXX = icpc
-FC = ifort
+CXX = clang++
+FC = flang
 
 ## Compiler flags
 H5CXXFLAGS = -O0 -g
-CXXFLAGS = -O0 -g -traceback
-FFLAGS = -O0 -g -traceback
+CXXFLAGS = -O0 -g
+FFLAGS = -O0 -g
 
 INCLUDE = -I $(INC_DIR)/
-DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o
+DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o
 DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o
 FLIBS = -lstdc++
 
@@ -66,6 +66,8 @@ $(OBJ_DIR)/%.o: $(TST_DIR)/%.f90 | $(OBJ_DIR)
 $(OBJ_DIR)/SM_MaponiA3.o: $(SRC_DIR)/SM_MaponiA3.cpp $(INC_DIR)/* | $(OBJ_DIR)
 	$(CXX) $(CXXFLAGS) -fPIC $(ARCH) $(INCLUDE) -c -o $@ $<
 
+$(OBJ_DIR)/SM_Standard.o: $(SRC_DIR)/SM_Standard.cpp $(INC_DIR)/* | $(OBJ_DIR)
+	$(CXX) $(CXXFLAGS) -fPIC $(ARCH) $(INCLUDE) -c -o $@ $<
 
 
 #### LINKING

include/Helpers.hpp

@@ -35,20 +35,21 @@ void showVector(T *vector, unsigned int size, string name) {
 
 template<typename T>
 void showMatrix(T *matrix, unsigned int M, string name) {
-    cout << name << " = " << endl;
+    cout.precision(17);
+    cout << name << " = [" << endl;
     for (unsigned int i = 0; i < M; i++) {
         cout << "[";
         for (unsigned int j = 0; j < M; j++) {
             if (matrix[i*M + j] >= 0) {
-                cout << "  " << matrix[i*M + j];
+                cout << "  " << matrix[i*M + j] << ",";
             }
             else {
-                cout << " " << matrix[i*M + j];
+                cout << " " << matrix[i*M + j] << "," ;
             }
         }
-        cout << " ]" << endl;
+        cout << " ]," << endl;
     }
-    cout << endl;
+    cout << "]" << endl;
 }
 
 template<typename T>

include/SM_Standard.hpp

@@ -0,0 +1,3 @@
+void SM(double *Slater_inv, unsigned int Dim,
+              unsigned int N_updates, double *Updates,
+              unsigned int *Updates_index);

src/SM_MaponiA3.cpp

@@ -5,7 +5,7 @@
 #include "Helpers.hpp"
 
 void MaponiA3(double *Slater_inv, unsigned int Dim,
-              unsigned int N_updates, double *Updates, 
+              unsigned int N_updates, double *Updates,
               unsigned int *Updates_index) {
 
   unsigned int k, l, lbar, i, j, tmp, component;
@@ -27,7 +27,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
       ylk[l][k] = new double[Dim + 1] {0};
     }
   }
-  
+
   // Calculate the y0k
   for (k = 1; k < N_updates + 1; k++) {
     for (i = 1; i < Dim + 1; i++) {
@@ -55,9 +55,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     p[lbar] = tmp;
     component = Updates_index[p[l] - 1];
     beta = 1 + ylk[l - 1][p[l]][component];
-    if (beta == 0) {
+    if (fabs(beta) < 1e-6) {
       cout << "Break-down occured. Exiting..." << endl;
-      exit;
+      exit(1);
     }
     for (k = l + 1; k < N_updates + 1; k++) {
       alpha = ylk[l - 1][p[k]][component] / beta;
@@ -110,7 +110,7 @@ extern "C" {
   void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
                   unsigned int *N_updates, double **linUpdates,
                   unsigned int **Updates_index) {
-    MaponiA3(*linSlater_inv, *Dim, 
+    MaponiA3(*linSlater_inv, *Dim,
              *N_updates, *linUpdates,
              *Updates_index);
   }

src/SM_Standard.cpp

@@ -0,0 +1,44 @@
+// SM-Standard.cpp
+// Standard Sherman Morrison with multiple updates
+#include "SM_Standard.hpp"
+#include "Helpers.hpp"
+#include <iostream>
+
+void SM(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+        double *Updates, unsigned int *Updates_index) {
+
+  double C[Dim];
+  double D[Dim];
+
+  // For each update
+  for (unsigned int l = 0; l < N_updates; l++) {
+
+    // C = A^{-1} x U_l
+    for (unsigned int i = 0; i < Dim; i++) {
+      C[i] = 0;
+      for (unsigned int j = 0; j < Dim; j++) {
+        C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
+      }
+    }
+
+    // Denominator
+    double den = 1 + C[Updates_index[l] - 1];
+    if (den < 1e-6) {
+      std::cerr << "Breakdown condition triggered" << std::endl;
+    }
+    double iden = 1 / den;
+
+    // D = v^T x A^{-1}
+    for (unsigned int j = 0; j < Dim; j++) {
+      D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
+    }
+
+    // A^{-1} = A^{-1} - C x D / den
+    for (unsigned int i = 0; i < Dim; i++) {
+      for (unsigned int j = 0; j < Dim; j++) {
+        double update = C[i] * D[j] * iden;
+        Slater_inv[i * Dim + j] -= update;
+      }
+    }
+  }
+}

tests/test_external_h5.cpp

@@ -4,10 +4,11 @@
 #include "hdf5/serial/H5Cpp.h"
 
 #include "SM_MaponiA3.hpp"
+#include "SM_Standard.hpp"
 #include "Helpers.hpp"
 
 using namespace H5;
-//#define DEBUG 0
+//#define DEBUG
 
 const H5std_string FILE_NAME( "datasets.hdf5" );
 
@@ -62,7 +63,8 @@ int test_cycle(H5File file, int cycle) {
     }
   }
 
-  MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
+  //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
+  SM(slater_inverse, dim, nupdates, updates, col_update_index);
 
 #ifdef DEBUG
   showMatrix(slater_matrix, dim, "NEW Slater");

tests/test_internal_h5.cpp

@@ -4,6 +4,7 @@
 #include "hdf5/serial/H5Cpp.h"
 
 #include "SM_MaponiA3.hpp"
+#include "SM_Standard.hpp"
 #include "Helpers.hpp"
 
 using namespace H5;
@@ -63,6 +64,7 @@ int test_cycle(H5File file, int cycle) {
   }
 
   MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index);
+  //SM(slater_inverse, dim, nupdates, updates, col_update_index);
 
 #ifdef DEBUG
   showMatrix(slater_matrix, dim, "NEW Slater");