First implementation of SM4

SM4: Sherman Morrison, mix between SM3 + SM2
Leave zero denominators for later (SM3), and when none are left then split (SM2)

include/SM_Standard.hpp

@@ -10,3 +10,8 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 // Sherman Morrison, leaving zero denominators for later
 void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index);
+
+// Sherman Morrison (SM3+SM2), leaving zero denominators for later (SM3), and
+// when none are left falling back on Splitting (SM2)
+void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+         double *Updates, unsigned int *Updates_index);

src/SM_Standard.cpp

@@ -172,6 +172,70 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   }
 }
 
+// Sherman Morrison, mix between SM3 + SM2
+// Leave zero denominators for later (SM3), and when none are left then split (SM2)
+void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+         double *Updates, unsigned int *Updates_index) {
+  std::cerr << "Called SM4 with updates " << N_updates << std::endl;
+  double C[Dim];
+  double D[Dim];
+
+  double later_updates[Dim * N_updates];
+  unsigned int later_index[N_updates];
+  unsigned int later = 0;
+
+  unsigned int l = 0;
+  // For each update
+  while (l < N_updates) {
+    // 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 (fabs(den) < threshold()) {
+      std::cerr << "Breakdown condition triggered at " << Updates_index[l]
+                << std::endl;
+
+      for (unsigned int j = 0; j < Dim; j++) {
+        later_updates[later * Dim + j] = Updates[l * Dim + j];
+      }
+      later_index[later] = Updates_index[l];
+      later++;
+      l += 1;
+      continue;
+    }
+    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;
+      }
+    }
+    l += 1;
+  }
+
+  // If all the updates have failed, fall back on splitting (SM2)
+  if (later == N_updates) {
+    SM2(Slater_inv, Dim, later, later_updates, later_index);
+  }
+  // If some have failed, make a recursive call
+  else if (later > 0) {
+    SM4(Slater_inv, Dim, later, later_updates, later_index);
+  }
+}
+
 extern "C" {
   void SM1_f(double **linSlater_inv, unsigned int *Dim,
                   unsigned int *N_updates, double **linUpdates,
@@ -184,10 +248,15 @@ extern "C" {
                   unsigned int **Updates_index) {
     SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
   }
-  
+
   void SM3_f(double **linSlater_inv, unsigned int *Dim,
                   unsigned int *N_updates, double **linUpdates,
                   unsigned int **Updates_index) {
     SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
   }
+  void SM4_f(double **linSlater_inv, unsigned int *Dim,
+                  unsigned int *N_updates, double **linUpdates,
+                  unsigned int **Updates_index) {
+    SM4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
 }

src/SM_mod.f90

@@ -36,5 +36,12 @@ module Sherman_Morrison
             real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
             real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
         end subroutine SM3
+        subroutine SM4(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM4_f")
+            use, intrinsic :: iso_c_binding, only : c_int, c_double
+            integer(c_int), intent(in) :: dim, n_updates
+            integer(c_int), dimension(:), allocatable, intent(in) :: Updates_index
+            real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
+            real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
+        end subroutine SM4
     end interface
 end module Sherman_Morrison

tests/test_h5.cpp

@@ -104,6 +104,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
     SM2(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "sm3") {
     SM3(slater_inverse, dim, nupdates, u, col_update_index);
+  } else if (version == "sm4") {
+    SM4(slater_inverse, dim, nupdates, u, col_update_index);
   } else {
     std::cerr << "Unknown version " << version << std::endl;
     exit(1);

tests/vfc_test_h5.cpp

@@ -137,6 +137,8 @@ int test_cycle(H5File file, int cycle, std::string version, vfc_probes * probes)
     SM2(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "sm3") {
     SM3(slater_inverse, dim, nupdates, u, col_update_index);
+  } else if (version == "sm3") {
+    SM4(slater_inverse, dim, nupdates, u, col_update_index);
   } else {
     std::cerr << "Unknown version " << version << std::endl;
     exit(1);