diff --git a/include/SMWB.hpp b/include/SMWB.hpp
index 28adda4..7ba4fdd 100644
--- a/include/SMWB.hpp
+++ b/include/SMWB.hpp
@@ -3,15 +3,15 @@
 void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index);
 
-// Sherman-Morrison-Woodbury kernel 2
-// WB2, WB3, SM3 mixing scheme
-void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-         double *Updates, unsigned int *Updates_index);
+// // Sherman-Morrison-Woodbury kernel 2
+// // WB2, WB3, SM3 mixing scheme
+// void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+//          double *Updates, unsigned int *Updates_index);
 
-// Sherman-Morrison-Woodbury kernel 3
-// WB2, WB3, SM4 mixing scheme
-void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-         double *Updates, unsigned int *Updates_index);
+// // Sherman-Morrison-Woodbury kernel 3
+// // WB2, WB3, SM4 mixing scheme
+// void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+//          double *Updates, unsigned int *Updates_index);
 
 // Sherman-Morrison-Woodbury kernel 4
 // WB2, SM2 mixing scheme
diff --git a/include/SM_Standard.hpp b/include/SM_Standard.hpp
index 4bfabd0..f68b122 100644
--- a/include/SM_Standard.hpp
+++ b/include/SM_Standard.hpp
@@ -7,6 +7,9 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index);
 
+void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+         double *Updates, unsigned int *Updates_index, double *later_updates, unsigned int* later_index, unsigned int *later);
+
 // 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);
diff --git a/src/SMWB.cpp b/src/SMWB.cpp
index b354a2d..b45b093 100644
--- a/src/SMWB.cpp
+++ b/src/SMWB.cpp
@@ -3,6 +3,9 @@
 #include "Woodbury.hpp"
 #include "Helpers.hpp"
 
+// #define DEBUG1
+// #define DEBUG2
+
 // Sherman-Morrison-Woodbury kernel 1
 // WB2, WB3, SM2 mixing scheme
 void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
@@ -19,6 +22,9 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double
   unsigned int length_1block = 1 * Dim;
 
   // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel
+  double later_updates[Dim * N_updates];
+  unsigned int later_index[N_updates];
+  unsigned int later = 0;
   if (n_of_3blocks > 0) {
     for (unsigned int i = 0; i < n_of_3blocks; i++) {
       double *Updates_3block = &Updates[i * length_3block];
@@ -26,12 +32,14 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double
       bool ok;
       ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
       if (!ok) { // Send the entire block to SM2
-#ifdef DEBUG2      
+        #ifdef DEBUG2      
         std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM2" << std::endl;
         showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
         showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
-#endif
-        SM2(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block);
+        #endif
+        unsigned int l = 0;
+        SM2star(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block, later_updates + (Dim * later), later_index + later, &l);
+        later = later + l;
       }
     }
   }
@@ -42,125 +50,132 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double
     bool ok;
     ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
     if (!ok) { // Send the entire block to SM2
-#ifdef DEBUG2    
+      #ifdef DEBUG2    
       std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl;
-#endif      
-      SM2(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
+      #endif      
+      unsigned int l = 0;
+      SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block, later_updates+(Dim * later), later_index + later, &l);
+      later = later + l;
     }
   } else if (remainder == 1) { // Apply last remaining update with SM2
     double *Updates_1block = &Updates[n_of_3blocks * length_3block];
     unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
-    SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
-  } else { // remainder == 0
-    // Nothing left to do.
+    unsigned int l = 0;
+    SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block, later_updates+(Dim * later), later_index + later, &l);
+    later = later + l;
   }
+
+  if (later > 0) {
+    SM2(Slater_inv, Dim, later, later_updates, later_index);
+  }
+
 }
 
-// Sherman-Morrison-Woodbury kernel 2
-// WB2, WB3, SM3 mixing scheme
-void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
-#ifdef DEBUG2
-  std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl;
-#endif
+// // Sherman-Morrison-Woodbury kernel 2
+// // WB2, WB3, SM3 mixing scheme
+// void SMWB2(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
+// #ifdef DEBUG2
+//   std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl;
+// #endif
 
-  unsigned int n_of_3blocks = N_updates / 3;
-  unsigned int remainder = N_updates % 3;
-  unsigned int length_3block = 3 * Dim;
-  unsigned int length_2block = 2 * Dim;
-  unsigned int length_1block = 1 * Dim;
+//   unsigned int n_of_3blocks = N_updates / 3;
+//   unsigned int remainder = N_updates % 3;
+//   unsigned int length_3block = 3 * Dim;
+//   unsigned int length_2block = 2 * Dim;
+//   unsigned int length_1block = 1 * Dim;
 
-#ifdef DEBUG2
-  showMatrix2(Updates_index, 1, N_updates, "Updates_index");
-  showMatrix2(Updates, N_updates, Dim, "Updates");
-#endif
+// #ifdef DEBUG2
+//   showMatrix2(Updates_index, 1, N_updates, "Updates_index");
+//   showMatrix2(Updates, N_updates, Dim, "Updates");
+// #endif
 
-  // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel
-  if (n_of_3blocks > 0) {
-    for (unsigned int i = 0; i < n_of_3blocks; i++) {
-      double *Updates_3block = &Updates[i * length_3block];
-      unsigned int *Updates_index_3block = &Updates_index[i * 3];
-      bool ok;
-      ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
-      if (!ok) { // Send the entire block to SM3
-#ifdef DEBUG2      
-        std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM3" << std::endl;
-        showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
-        showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
-#endif
-        SM3(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block);
-      }
-    }
-  }
+//   // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel
+//   if (n_of_3blocks > 0) {
+//     for (unsigned int i = 0; i < n_of_3blocks; i++) {
+//       double *Updates_3block = &Updates[i * length_3block];
+//       unsigned int *Updates_index_3block = &Updates_index[i * 3];
+//       bool ok;
+//       ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
+//       if (!ok) { // Send the entire block to SM3
+// #ifdef DEBUG2      
+//         std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM3" << std::endl;
+//         showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
+//         showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
+// #endif
+//         SM3(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block);
+//       }
+//     }
+//   }
 
-  if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
-    double *Updates_2block = &Updates[n_of_3blocks * length_3block];
-    unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
-    bool ok;
-    ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
-    if (!ok) { // Send the entire block to SM3
-      std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM3" << std::endl;
-      SM3(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
-    }
-  } else if (remainder == 1) { // Apply last remaining update with SM3
-    double *Updates_1block = &Updates[n_of_3blocks * length_3block];
-    unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
-    SM3(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
-  } else { // remainder == 0
-    // Nothing left to do.
-  }
-}
+//   if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
+//     double *Updates_2block = &Updates[n_of_3blocks * length_3block];
+//     unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
+//     bool ok;
+//     ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
+//     if (!ok) { // Send the entire block to SM3
+//       std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM3" << std::endl;
+//       SM3(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
+//     }
+//   } else if (remainder == 1) { // Apply last remaining update with SM3
+//     double *Updates_1block = &Updates[n_of_3blocks * length_3block];
+//     unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
+//     SM3(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
+//   } else { // remainder == 0
+//     // Nothing left to do.
+//   }
+// }
 
-// Sherman-Morrison-Woodbury kernel 3
-// WB2, WB3, SM4 mixing scheme
-void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
-           std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl;
+// // Sherman-Morrison-Woodbury kernel 3
+// // WB2, WB3, SM4 mixing scheme
+// void SMWB3(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
+//            std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates << " updates" << std::endl;
 
-  unsigned int n_of_3blocks = N_updates / 3;
-  unsigned int remainder = N_updates % 3;
-  unsigned int length_3block = 3 * Dim;
-  unsigned int length_2block = 2 * Dim;
-  unsigned int length_1block = 1 * Dim;
+//   unsigned int n_of_3blocks = N_updates / 3;
+//   unsigned int remainder = N_updates % 3;
+//   unsigned int length_3block = 3 * Dim;
+//   unsigned int length_2block = 2 * Dim;
+//   unsigned int length_1block = 1 * Dim;
 
-#ifdef DEBUG2
-  showMatrix2(Updates_index, 1, N_updates, "Updates_index");
-  showMatrix2(Updates, N_updates, Dim, "Updates");
-#endif
+// #ifdef DEBUG2
+//   showMatrix2(Updates_index, 1, N_updates, "Updates_index");
+//   showMatrix2(Updates, N_updates, Dim, "Updates");
+// #endif
 
-  // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel
-  if (n_of_3blocks > 0) {
-    for (unsigned int i = 0; i < n_of_3blocks; i++) {
-      double *Updates_3block = &Updates[i * length_3block];
-      unsigned int *Updates_index_3block = &Updates_index[i * 3];
-      bool ok;
-      ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
-      if (!ok) { // Send the entire block to SM4
-        std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM4" << std::endl;
-#ifdef DEBUG2        
-        showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
-        showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
-#endif
-        SM4(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block);
-      }
-    }
-  }
+//   // Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with Woodbury 3x3 kernel
+//   if (n_of_3blocks > 0) {
+//     for (unsigned int i = 0; i < n_of_3blocks; i++) {
+//       double *Updates_3block = &Updates[i * length_3block];
+//       unsigned int *Updates_index_3block = &Updates_index[i * 3];
+//       bool ok;
+//       ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
+//       if (!ok) { // Send the entire block to SM4
+//         std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM4" << std::endl;
+// #ifdef DEBUG2        
+//         showMatrix2(Updates_3block, 3, Dim, "Updates_3block");
+//         showMatrix2(Updates_index_3block, 1, 3, "Updates_index_3block");
+// #endif
+//         SM4(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block);
+//       }
+//     }
+//   }
 
-  if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
-    double *Updates_2block = &Updates[n_of_3blocks * length_3block];
-    unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
-    bool ok;
-    ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
-    if (!ok) { // Send the entire block to SM4
-      std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM4" << std::endl;
-      SM4(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
-    }
-  } else if (remainder == 1) { // Apply last remaining update with SM4
-    double *Updates_1block = &Updates[n_of_3blocks * length_3block];
-    unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
-    SM4(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
-  } else { // remainder == 0
-    // Nothing left to do.
-  }
-}
+//   if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
+//     double *Updates_2block = &Updates[n_of_3blocks * length_3block];
+//     unsigned int *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
+//     bool ok;
+//     ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
+//     if (!ok) { // Send the entire block to SM4
+//       std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM4" << std::endl;
+//       SM4(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
+//     }
+//   } else if (remainder == 1) { // Apply last remaining update with SM4
+//     double *Updates_1block = &Updates[n_of_3blocks * length_3block];
+//     unsigned int *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
+//     SM4(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
+//   } else { // remainder == 0
+//     // Nothing left to do.
+//   }
+// }
 
 // Sherman-Morrison-Woodbury kernel 4
 // WB2, SM2 mixing scheme
@@ -177,6 +192,9 @@ void SMWB4(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double
   unsigned int length_1block = 1 * Dim;
 
   // Apply first 2*n_of_2blocks updates in n_of_2blocks blocks of 2 updates with Woodbury 2x2 kernel
+  double later_updates[Dim * N_updates];
+  unsigned int later_index[N_updates];
+  unsigned int later = 0;
   if (n_of_2blocks > 0) {
     for (unsigned int i = 0; i < n_of_2blocks; i++) {
       double *Updates_2block = &Updates[i * length_2block];
@@ -185,22 +203,29 @@ void SMWB4(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double
       ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
       if (!ok) { // Send the entire block to SM2
         std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl;
-#ifdef DEBUG2
+        #ifdef DEBUG2
         showMatrix2(Updates_2block, 2, Dim, "Updates_2block");
         showMatrix2(Updates_index_2block, 1, 2, "Updates_index_2block");
-#endif
-        SM2(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
+        #endif
+        unsigned int l = 0;
+        SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block, later_updates + (Dim * later), later_index + later, &l);
+        later = later + l;
       }
     }
   }
 
-  if (remainder == 1) { // Apply last remaining update with SM4
+  if (remainder == 1) { // Apply last remaining update with SM2
     double *Updates_1block = &Updates[n_of_2blocks * length_2block];
     unsigned int *Updates_index_1block = &Updates_index[2 * n_of_2blocks];
-    SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
-  } else { // remainder == 0
-    // Nothing left to do.
+        unsigned int l = 0;
+        SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block, later_updates + (Dim * later), later_index + later, &l);
+        later = later + l;
   }
+  
+  if (later > 0) {
+    SM2(Slater_inv, Dim, later, later_updates, later_index);
+  }
+
 }
 
 extern "C" {
@@ -208,14 +233,14 @@ void SMWB1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
            double **linUpdates, unsigned int **Updates_index) {
   SMWB1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
 }
-void SMWB2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
-           double **linUpdates, unsigned int **Updates_index) {
-  SMWB2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
-}
-void SMWB3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
-           double **linUpdates, unsigned int **Updates_index) {
-  SMWB3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
-}
+// void SMWB2_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
+//            double **linUpdates, unsigned int **Updates_index) {
+//   SMWB2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+// }
+// void SMWB3_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
+//            double **linUpdates, unsigned int **Updates_index) {
+//   SMWB3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+// }
 void SMWB4_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
            double **linUpdates, unsigned int **Updates_index) {
   SMWB4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp
index 7423ce8..b2a40dd 100644
--- a/src/SM_Standard.cpp
+++ b/src/SM_Standard.cpp
@@ -3,6 +3,8 @@
 #include "SM_Standard.hpp"
 #include "Helpers.hpp"
 
+#define DEBUG1
+
 // Naïve Sherman Morrison
 void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
@@ -57,13 +59,13 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl;
 #endif
 
-  double C[Dim];
-  double D[Dim];
-
   double later_updates[Dim * N_updates];
   unsigned int later_index[N_updates];
   unsigned int later = 0;
 
+  double C[Dim];
+  double D[Dim];
+
   unsigned int l = 0;
   // For each update
   while (l < N_updates) {
@@ -116,6 +118,66 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   }
 }
 
+// Sherman Morrison, with J. Slagel splitting
+// http://hdl.handle.net/10919/52966
+void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+         double *Updates, unsigned int *Updates_index, double *later_updates, unsigned int* later_index, unsigned int *later) {
+#ifdef DEBUG1
+  std::cerr << "Called SM2* with " << N_updates << " updates" << std::endl;
+#endif
+
+  double C[Dim];
+  double D[Dim];
+
+  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 (std::fabs(den) < threshold()) {
+#ifdef DEBUG1
+      std::cerr << "Breakdown condition triggered at " << Updates_index[l]
+                << std::endl;
+      std::cerr << "Denominator = " << den << std::endl;
+#endif
+
+      // U_l = U_l / 2 (do the split)
+      for (unsigned int i = 0; i < Dim; i++) {
+        later_updates[*later * Dim + i] = Updates[l * Dim + i] / 2.0;
+        C[i] /= 2.0;
+      }
+      later_index[*later] = Updates_index[l];
+      (*later)++;
+
+      den = 1 + C[Updates_index[l] - 1];
+    }
+    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;
+  }
+}
+
+
 // 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) {
diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index f6573f2..b5b299a 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -10,6 +10,9 @@
 #include "Woodbury.hpp"
 #include "Helpers.hpp"
 
+// #define DEBUG1
+// #define DEBUG2
+
 // Woodbury 2x2 kernel
 bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
          unsigned int *Updates_index) {
@@ -18,7 +21,7 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
       B := 1 + V * C,     2 x 2
       D := V * S^{-1},    2 x dim
   */
- #ifdef DEBUG1
+#ifdef DEBUG1
   std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
 #endif
 
@@ -43,7 +46,6 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
       }
     }
   }
-  // matMul2(Updates, Slater_inv, C, 2, Dim, Dim);
 
   // Compute B = 1 + V * C
   double B[4];
@@ -53,24 +55,23 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   B[3] = C[row2 * 2 + 1];
   B[0] += 1, B[3] += 1;
 
-  // Compute B^{-1} with explicit formula for 2x2 inversion
-  double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]);
-  double Binv[4];
-  Binv[0] = idet * B[3];
-  Binv[1] = -1.0 * idet * B[1];
-  Binv[2] = -1.0 * idet * B[2];
-  Binv[3] = idet * B[0];
-
   // Check if determinant of inverted matrix is not zero
-  double det = Binv[0] * Binv[3] - Binv[1] * Binv[2];
+  double det = B[0] * B[3] - B[1] * B[2];
   if (std::fabs(det) < threshold()) {
+    std::cerr << "Determinant too close to zero! No inverse found." << std::endl;
 #ifdef DEBUG1
-    std::cerr << "Determinant too close to zero!" << std::endl;
     std::cerr << "Determinant = " << det << std::endl;
 #endif
     return false;
   }
 
+  // Compute B^{-1} with explicit formula for 2x2 inversion
+  double Binv[4], idet = 1.0 / det;
+  Binv[0] = idet * B[3];
+  Binv[1] = -1.0 * idet * B[1];
+  Binv[2] = -1.0 * idet * B[2];
+  Binv[3] = idet * B[0];
+
   // Compute B^{-1} x D
   double tmp[2 * Dim];
   matMul2(Binv, D, tmp, 2, 2, Dim);
@@ -127,7 +128,6 @@ bool  WB3(double *Slater_inv, unsigned int Dim, double *Updates,
       }
     }
   }
-  // matMul2(Updates, Slater_inv, C, 2, Dim, Dim);
 
 #ifdef DEBUG2
   showMatrix2(C, Dim, 3, "C = S_inv * U");
@@ -154,17 +154,25 @@ bool  WB3(double *Slater_inv, unsigned int Dim, double *Updates,
   showMatrix2(B, 3, 3, "B = 1 + V * C");
 #endif
 
-  // Compute B^{-1} with explicit formula for 3x3 inversion
-  double Binv[9], det, idet;
+  // Check if determinant of B is not too close to zero
+  double det;
   det = B[0] * (B[4] * B[8] - B[5] * B[7]) -
         B[1] * (B[3] * B[8] - B[5] * B[6]) +
         B[2] * (B[3] * B[7] - B[4] * B[6]);
-  idet = 1.0 / det;
-
 #ifdef DEBUG2 
   std::cerr << "Determinant of B = " << det << std::endl;
 #endif
+  if (std::fabs(det) < threshold()) {
+  // if (std::fabs(det) < 1000000) {
+    std::cerr << "Determinant too close to zero! No inverse found." << std::endl;
+#ifdef DEBUG1
+    std::cerr << "Determinant = " << det << std::endl;
+#endif
+    return false;
+  }
 
+  // Compute B^{-1} with explicit formula for 3x3 inversion
+  double Binv[9], idet = 1.0 / det;
   Binv[0] =   ( B[4] * B[8] - B[7] * B[5] ) * idet;
   Binv[1] = - ( B[1] * B[8] - B[7] * B[2] ) * idet;
   Binv[2] =   ( B[1] * B[5] - B[4] * B[2] ) * idet;
@@ -180,24 +188,6 @@ bool  WB3(double *Slater_inv, unsigned int Dim, double *Updates,
   showMatrix2(Binv, 3, 3, "Binv");
 #endif
 
-  // Check if determinant of inverted matrix is not zero
-  // double det;
-  det = Binv[0] * (Binv[4] * Binv[8] - Binv[5] * Binv[7]) -
-        Binv[1] * (Binv[3] * Binv[8] - Binv[5] * Binv[6]) +
-        Binv[2] * (Binv[3] * Binv[7] - Binv[4] * Binv[6]);
-
-#ifdef DEBUG2
-  std::cerr << "Determinant of Binv = " << det << std::endl;
-#endif
-
-  if (std::fabs(det) < threshold()) {
-#ifdef DEBUG1
-    std::cerr << "Determinant too close to zero!" << std::endl;
-    std::cerr << "Determinant = " << det << std::endl;
-#endif
-    return false;
-  }
-
   // Compute B^{-1} x D
   double tmp[3 * Dim];
   matMul2(Binv, D, tmp, 3, 3, Dim);
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index 9cb296b..c174f36 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -7,7 +7,7 @@
 #include "Woodbury.hpp"
 #include "SMWB.hpp"
 
-#define PERF
+// #define PERF
 
 #ifdef PERF
 unsigned int repetition_number;
@@ -101,10 +101,10 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
       WB3(slater_inverse_nonpersistent, dim, u, col_update_index);
     } else if (version == "smwb1") {
       SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
-    } else if (version == "smwb2") {
-      SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
-    } else if (version == "smwb3") {
-      SMWB3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
+    // } else if (version == "smwb2") {
+    //   SMWB2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
+    // } else if (version == "smwb3") {
+    //   SMWB3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     } else if (version == "smwb4") {
       SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
 #ifdef MKL
@@ -138,10 +138,10 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
     WB3(slater_inverse, dim, u, col_update_index);
   } else if (version == "smwb1") {
     SMWB1(slater_inverse, dim, nupdates, u, col_update_index);
-  } else if (version == "smwb2") {
-    SMWB2(slater_inverse, dim, nupdates, u, col_update_index);
-  } else if (version == "smwb3") {
-    SMWB3(slater_inverse, dim, nupdates, u, col_update_index);
+  // } else if (version == "smwb2") {
+  //   SMWB2(slater_inverse, dim, nupdates, u, col_update_index);
+  // } else if (version == "smwb3") {
+  //   SMWB3(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "smwb4") {
     SMWB4(slater_inverse, dim, nupdates, u, col_update_index);
 #ifdef MKL