diff --git a/include/Helpers.hpp b/include/Helpers.hpp
index c8eb3f8..41b320a 100644
--- a/include/Helpers.hpp
+++ b/include/Helpers.hpp
@@ -98,11 +98,12 @@ template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
 }
 
 template <typename T1, typename T2, typename T3>
-void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N, unsigned int P) {
-  for(unsigned int i = 0; i < M; i++) {
-    for(unsigned int j = 0; j < P; j++) {
+void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N,
+             unsigned int P) {
+  for (unsigned int i = 0; i < M; i++) {
+    for (unsigned int j = 0; j < P; j++) {
       C[i * P + j] = 0;
-      for(unsigned int k = 0; k < N; k++) {
+      for (unsigned int k = 0; k < N; k++) {
         C[i * P + j] += A[i * N + k] * B[k * P + j];
       }
     }
@@ -265,13 +266,12 @@ template <typename T> T residual_frobenius2(T *A, unsigned int Dim) {
   return res;
 }
 
-
-template <typename T> T residual2(T * A, unsigned int Dim) {
+template <typename T> T residual2(T *A, unsigned int Dim) {
   double res = 0.0;
   for (unsigned int i = 0; i < Dim; i++) {
     for (unsigned int j = 0; j < Dim; j++) {
       T delta = (A[i * Dim + j] - (i == j));
-      res += delta*delta;
+      res += delta * delta;
     }
   }
   return res;
diff --git a/include/SMWB.hpp b/include/SMWB.hpp
index 7ba4fdd..6d1d0f9 100644
--- a/include/SMWB.hpp
+++ b/include/SMWB.hpp
@@ -1,7 +1,8 @@
 // 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);
+void SMWB1(double *Slater_inv, const unsigned int Dim,
+           const unsigned int N_updates, double *Updates,
+           unsigned int *Updates_index);
 
 // // Sherman-Morrison-Woodbury kernel 2
 // // WB2, WB3, SM3 mixing scheme
@@ -15,5 +16,6 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 
 // Sherman-Morrison-Woodbury kernel 4
 // WB2, SM2 mixing scheme
-void SMWB4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-         double *Updates, unsigned int *Updates_index);
+void SMWB4(double *Slater_inv, const unsigned int Dim,
+           const unsigned int N_updates, double *Updates,
+           unsigned int *Updates_index);
diff --git a/include/SM_Standard.hpp b/include/SM_Standard.hpp
index f68b122..c004980 100644
--- a/include/SM_Standard.hpp
+++ b/include/SM_Standard.hpp
@@ -8,7 +8,9 @@ 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);
+             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,
diff --git a/include/Woodbury.hpp b/include/Woodbury.hpp
index c03e509..3457790 100644
--- a/include/Woodbury.hpp
+++ b/include/Woodbury.hpp
@@ -1,7 +1,7 @@
 // Woodbury 2x2 kernel
-bool WB2(double *Slater_inv, unsigned int Dim,
-         double *Updates, unsigned int *Updates_index);
+bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
+         unsigned int *Updates_index);
 
 // Woodbury 3x3 kernel
-bool WB3(double *Slater_inv, unsigned int Dim,
-         double *Updates, unsigned int *Updates_index);
+bool WB3(double *Slater_inv, unsigned int Dim, double *Updates,
+         unsigned int *Updates_index);
diff --git a/src/SMWB.cpp b/src/SMWB.cpp
index fe29d34..18c7d54 100644
--- a/src/SMWB.cpp
+++ b/src/SMWB.cpp
@@ -1,16 +1,19 @@
 #include "SMWB.hpp"
+#include "Helpers.hpp"
 #include "SM_Standard.hpp"
 #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) {
+void SMWB1(double *Slater_inv, const unsigned int Dim,
+           const 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;
+  std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates
+            << " updates" << std::endl;
   showMatrix2(Updates_index, 1, N_updates, "Updates_index");
   showMatrix2(Updates, N_updates, Dim, "Updates");
 #endif
@@ -21,7 +24,8 @@ void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates, double
   unsigned int length_2block = 2 * Dim;
   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
+  // 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;
@@ -32,50 +36,58 @@ 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      
-        std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM2" << std::endl;
+#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
+#endif
         unsigned int l = 0;
-        SM2star(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block, later_updates + (Dim * later), later_index + later, &l);
+        SM2star(Slater_inv, Dim, 3, Updates_3block, Updates_index_3block,
+                later_updates + (Dim * later), later_index + later, &l);
         later = later + l;
       }
     }
   }
 
-  if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
+  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 SM2
-      #ifdef DEBUG2    
-      std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2" << std::endl;
-      #endif      
+#ifdef DEBUG2
+      std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2"
+                << std::endl;
+#endif
       unsigned int l = 0;
-      SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block, later_updates+(Dim * later), later_index + later, &l);
+      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];
     unsigned int l = 0;
-    SM2star(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block, later_updates+(Dim * later), later_index + later, &l);
+    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 4
 // WB2, SM2 mixing scheme
-void SMWB4(double *Slater_inv, const unsigned int Dim, const unsigned int N_updates, double *Updates, unsigned int *Updates_index) {
+void SMWB4(double *Slater_inv, const unsigned int Dim,
+           const 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;
+  std::cerr << "Called Sherman-Morrison-Woodbury kernel 1 with " << N_updates
+            << " updates" << std::endl;
   showMatrix2(Updates_index, 1, N_updates, "Updates_index");
   showMatrix2(Updates, N_updates, Dim, "Updates");
 #endif
@@ -85,7 +97,8 @@ void SMWB4(double *Slater_inv, const unsigned int Dim, const unsigned int N_upda
   unsigned int length_2block = 2 * Dim;
   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
+  // 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;
@@ -96,15 +109,17 @@ void SMWB4(double *Slater_inv, const unsigned int Dim, const unsigned int N_upda
       bool ok;
       ok = WB2(Slater_inv, Dim, Updates_2block, Updates_index_2block);
       if (!ok) { // Send the entire block to SM2
-        #ifdef DEBUG1
-        std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2star" << std::endl;
-        #endif
-        #ifdef DEBUG2
+#ifdef DEBUG1
+        std::cerr << "Woodbury 2x2 kernel failed! Sending block to SM2star"
+                  << std::endl;
+#endif
+#ifdef DEBUG2
         showMatrix2(Updates_2block, 2, Dim, "Updates_2block");
         showMatrix2(Updates_index_2block, 1, 2, "Updates_index_2block");
-        #endif
+#endif
         unsigned int l = 0;
-        SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block, later_updates + (Dim * later), later_index + later, &l);
+        SM2star(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block,
+                later_updates + (Dim * later), later_index + later, &l);
         later = later + l;
       }
     }
@@ -113,24 +128,24 @@ void SMWB4(double *Slater_inv, const unsigned int Dim, const unsigned int N_upda
   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];
-        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;
+    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" {
 void SMWB1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
-           double **linUpdates, unsigned int **Updates_index) {
+             double **linUpdates, unsigned int **Updates_index) {
   SMWB1(*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) {
+             double **linUpdates, unsigned int **Updates_index) {
   SMWB4(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
 }
 }
diff --git a/src/SM_Maponi.cpp b/src/SM_Maponi.cpp
index 2348abc..a95a441 100644
--- a/src/SM_Maponi.cpp
+++ b/src/SM_Maponi.cpp
@@ -75,10 +75,10 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     std::cerr << std::endl;
 #endif
     if (std::fabs(beta) < threshold()) {
-      #ifdef DEBUG1
+#ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
-      #endif
+#endif
     }
     double ibeta = 1.0 / beta;
 
@@ -209,10 +209,10 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     std::cerr << std::endl;
 #endif
     if (std::fabs(beta) < threshold()) {
-      #ifdef DEBUG1
+#ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
-      #endif
+#endif
       for (unsigned int i = 1; i < Dim + 1; i++) {
         ylk[l][p[l + 1]][i] *= 0.5;
         later_updates[later * Dim + i - 1] = Updates[l * Dim + i - 1] * 0.5;
diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp
index dc9cc54..84b5c82 100644
--- a/src/SM_Standard.cpp
+++ b/src/SM_Standard.cpp
@@ -8,9 +8,9 @@
 // Naïve Sherman Morrison
 void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
-  #ifdef DEBUG1
+#ifdef DEBUG1
   std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
-  #endif
+#endif
 
   double C[Dim];
   double D[Dim];
@@ -29,10 +29,10 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
     if (std::fabs(den) < threshold()) {
-      #ifdef DEBUG1
+#ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
-      #endif
+#endif
     }
     double iden = 1 / den;
 
@@ -57,9 +57,9 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 // http://hdl.handle.net/10919/52966
 void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
-  #ifdef DEBUG1
+#ifdef DEBUG1
   std::cerr << "Called SM2 with " << N_updates << " updates" << std::endl;
-  #endif
+#endif
 
   double later_updates[Dim * N_updates];
   unsigned int later_index[N_updates];
@@ -82,11 +82,11 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
     if (std::fabs(den) < threshold()) {
-      #ifdef DEBUG1
+#ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
       std::cerr << "Denominator = " << den << std::endl;
-      #endif
+#endif
 
       // U_l = U_l / 2 (do the split)
       for (unsigned int i = 0; i < Dim; i++) {
@@ -123,10 +123,12 @@ 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
+             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
+#endif
 
   double C[Dim];
   double D[Dim];
@@ -145,11 +147,11 @@ void SM2star(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
     if (std::fabs(den) < threshold()) {
-      #ifdef DEBUG1
+#ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
       std::cerr << "Denominator = " << den << std::endl;
-      #endif
+#endif
 
       // U_l = U_l / 2 (do the split)
       for (unsigned int i = 0; i < Dim; i++) {
@@ -179,13 +181,12 @@ void SM2star(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) {
-  #ifdef DEBUG1
+#ifdef DEBUG1
   std::cerr << "Called SM3 with " << N_updates << " updates" << std::endl;
-  #endif
+#endif
 
   double C[Dim];
   double D[Dim];
@@ -208,10 +209,10 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
     if (std::fabs(den) < threshold()) {
-      #ifdef DEBUG1
+#ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
-      #endif
+#endif
       for (unsigned int j = 0; j < Dim; j++) {
         later_updates[later * Dim + j] = Updates[l * Dim + j];
       }
@@ -239,10 +240,10 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 
   // If all the updates have failed, exit early with an error
   if (later == N_updates) {
-    #ifdef DEBUG1
+#ifdef DEBUG1
     std::cerr << "SM3 cannot invert matrix." << std::endl;
     showMatrix(Slater_inv, Dim, "Slater_inverse");
-    #endif
+#endif
     return;
   }
   // If some have failed, make a recursive call
@@ -256,9 +257,9 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 // (SM2)
 void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
          double *Updates, unsigned int *Updates_index) {
-  #ifdef DEBUG1           
+#ifdef DEBUG1
   std::cerr << "Called SM4 with " << N_updates << " updates" << std::endl;
-  #endif
+#endif
 
   double C[Dim];
   double D[Dim];
@@ -281,10 +282,10 @@ void SM4(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
     if (std::fabs(den) < threshold()) {
-      #ifdef DEBUG1
+#ifdef DEBUG1
       std::cerr << "Breakdown condition triggered at " << Updates_index[l]
                 << std::endl;
-      #endif
+#endif
       for (unsigned int j = 0; j < Dim; j++) {
         later_updates[later * Dim + j] = Updates[l * Dim + j];
       }
diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index 37ebbf0..bd04633 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -16,11 +16,11 @@
 // Woodbury 2x2 kernel
 bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
          unsigned int *Updates_index) {
-  /*
-      C := S^{-1} * U,    dim x 2
-      B := 1 + V * C,     2 x 2
-      D := V * S^{-1},    2 x dim
-  */
+/*
+    C := S^{-1} * U,    dim x 2
+    B := 1 + V * C,     2 x 2
+    D := V * S^{-1},    2 x dim
+*/
 #ifdef DEBUG1
   std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
 #endif
@@ -28,13 +28,13 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   const unsigned int row1 = (Updates_index[0] - 1);
   const unsigned int row2 = (Updates_index[1] - 1);
 
-  // Compute C = S_inv * U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE 
+  // Compute C = S_inv * U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE
   // OF LAYOUT OF 'Updates' !!
   double C[2 * Dim];
-  for(unsigned int i = 0; i < Dim; i++) {
-    for(unsigned int j = 0; j < 2; j++) {
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < 2; j++) {
       C[i * 2 + j] = 0;
-      for(unsigned int k = 0; k < Dim; k++) {
+      for (unsigned int k = 0; k < Dim; k++) {
         C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k];
       }
     }
@@ -51,10 +51,11 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   // Check if determinant of inverted matrix is not zero
   double det = B[0] * B[3] - B[1] * B[2];
   if (std::fabs(det) < threshold()) {
-    #ifdef DEBUG1
-    std::cerr << "Determinant too close to zero! No inverse found." << std::endl;
+#ifdef DEBUG1
+    std::cerr << "Determinant too close to zero! No inverse found."
+              << std::endl;
     std::cerr << "Determinant = " << det << std::endl;
-    #endif
+#endif
     return false;
   }
 
@@ -67,17 +68,17 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
 
   // Compute tmp = B^{-1} x (V.S^{-1})
   double tmp[2 * Dim];
-  for(unsigned int i = 0; i < 2; i++) {
-    for(unsigned int j = 0; j < Dim; j++) {
-      tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1*Dim + j];
-      tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2*Dim + j];
+  for (unsigned int i = 0; i < 2; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1 * Dim + j];
+      tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2 * Dim + j];
     }
   }
 
   // Compute (S + U V)^{-1} = S^{-1} - C x tmp
-  for(unsigned int i = 0; i < Dim; i++) {
-    for(unsigned int j = 0; j < Dim; j++) {
-      Slater_inv[i * Dim + j] -= C[i * 2 ] * tmp[j];
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      Slater_inv[i * Dim + j] -= C[i * 2] * tmp[j];
       Slater_inv[i * Dim + j] -= C[i * 2 + 1] * tmp[Dim + j];
     }
   }
@@ -86,41 +87,28 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
 }
 
 // Woodbury 3x3 kernel
-bool  WB3(double *Slater_inv, unsigned int Dim, double *Updates,
+bool WB3(double *Slater_inv, unsigned int Dim, double *Updates,
          unsigned int *Updates_index) {
-  /*
-      C := S^{-1} * U,    dim x 3
-      B := 1 + V * C,     3 x 3
-      D := V * S^{-1},    3 x dim
-  */
+/*
+    C := S^{-1} * U,    dim x 3
+    B := 1 + V * C,     3 x 3
+    D := V * S^{-1},    3 x dim
+*/
 #ifdef DEBUG1
   std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
 #endif
-#ifdef DEBUG2
-  showMatrix2(Slater_inv, Dim, Dim, "Slater_inv BEFORE update");
-  showMatrix2(Updates, 3, Dim, "Updates");
-  showMatrix2(Updates_index, 1, 3, "Updates_index");
-#endif
 
-  // Compute D = V * S^{-1}
-  double D[3 * Dim];
-  unsigned int row1, row2, row3;
-  row1 = Updates_index[0] - 1;
-  row2 = Updates_index[1] - 1;
-  row3 = Updates_index[2] - 1;
-  for (unsigned int i = 0; i < Dim; i++) {
-    D[i] = Slater_inv[row1 * Dim + i];
-    D[Dim + i] = Slater_inv[row2 * Dim + i];
-    D[2 * Dim + i] = Slater_inv[row3 * Dim + i];
-  }
+  const unsigned int row1 = (Updates_index[0] - 1);
+  const unsigned int row2 = (Updates_index[1] - 1);
+  const unsigned int row3 = (Updates_index[2] - 1);
 
-  // Compute C = S_inv * U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE 
+  // Compute C = S_inv * U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE
   // OF LAYOUT OF 'Updates' !!
   double C[3 * Dim];
-  for(unsigned int i = 0; i < Dim; i++) {
-    for(unsigned int j = 0; j < 3; j++) {
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < 3; j++) {
       C[i * 3 + j] = 0;
-      for(unsigned int k = 0; k < Dim; k++) {
+      for (unsigned int k = 0; k < Dim; k++) {
         C[i * 3 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k];
       }
     }
@@ -128,9 +116,6 @@ bool  WB3(double *Slater_inv, unsigned int Dim, double *Updates,
 
 #ifdef DEBUG2
   showMatrix2(C, Dim, 3, "C = S_inv * U");
-#endif
-
-#ifdef DEBUG2
   showMatrix2(D, 3, Dim, "D = V * S_inv");
 #endif
 
@@ -154,55 +139,58 @@ bool  WB3(double *Slater_inv, unsigned int Dim, double *Updates,
   // 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]);
-#ifdef DEBUG2 
+        B[1] * (B[3] * B[8] - B[5] * B[6]) + B[2] * (B[3] * B[7] - B[4] * B[6]);
+#ifdef DEBUG2
   std::cerr << "Determinant of B = " << det << std::endl;
 #endif
   if (std::fabs(det) < threshold()) {
-    #ifdef DEBUG1
-    std::cerr << "Determinant too close to zero! No inverse found." << std::endl;
+#ifdef DEBUG1
+    std::cerr << "Determinant too close to zero! No inverse found."
+              << std::endl;
     std::cerr << "Determinant = " << det << std::endl;
-    #endif
+#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;
-  Binv[3] = - ( B[3] * B[8] - B[6] * B[5] ) * idet;
-  Binv[4] =   ( B[0] * B[8] - B[6] * B[2] ) * idet;
-  Binv[5] = - ( B[0] * B[5] - B[3] * B[2] ) * idet;
-  Binv[6] =   ( B[3] * B[7] - B[6] * B[4] ) * idet;
-  Binv[7] = - ( B[0] * B[7] - B[6] * B[1] ) * idet;
-  Binv[8] =   ( B[0] * B[4] - B[3] * B[1] ) * idet;
+  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;
+  Binv[3] = -(B[3] * B[8] - B[6] * B[5]) * idet;
+  Binv[4] = (B[0] * B[8] - B[6] * B[2]) * idet;
+  Binv[5] = -(B[0] * B[5] - B[3] * B[2]) * idet;
+  Binv[6] = (B[3] * B[7] - B[6] * B[4]) * idet;
+  Binv[7] = -(B[0] * B[7] - B[6] * B[1]) * idet;
+  Binv[8] = (B[0] * B[4] - B[3] * B[1]) * idet;
 
 #ifdef DEBUG2
-  std::cerr << "Conditioning number of B = " << condition1(B, Binv, 3) << std::endl;
+  std::cerr << "Conditioning number of B = " << condition1(B, Binv, 3)
+            << std::endl;
   showMatrix2(Binv, 3, 3, "Binv");
 #endif
 
-  // Compute B^{-1} x D
+  // Compute tmp = B^{-1} x (V.S^{-1})
   double tmp[3 * Dim];
-  matMul2(Binv, D, tmp, 3, 3, Dim);
+  for (unsigned int i = 0; i < 3; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      tmp[i * Dim + j] = Binv[i * 3] * Slater_inv[row1 * Dim + j];
+      tmp[i * Dim + j] += Binv[i * 3 + 1] * Slater_inv[row2 * Dim + j];
+      tmp[i * Dim + j] += Binv[i * 3 + 2] * Slater_inv[row3 * Dim + j];
+    }
+  }
 
 #ifdef DEBUG2
   showMatrix2(tmp, 3, Dim, "tmp = Binv * D");
 #endif
 
-  // Compute C x B^{-1} x D
-  double tmp2[Dim * Dim];
-  matMul2(C, tmp, tmp2, Dim, 3, Dim);
-
-#ifdef DEBUG2
-  showMatrix2(tmp2, Dim, Dim, "tmp2 = C * tmp");
-#endif
-
-  // Compute (S + U V)^{-1} = S^{-1} - C B^{-1} D
-  for (unsigned int i = 0; i < Dim * Dim; i++) {
-    Slater_inv[i] -= tmp2[i];
+  // Compute (S + U V)^{-1} = S^{-1} - C x tmp
+  for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int j = 0; j < Dim; j++) {
+      Slater_inv[i * Dim + j] -= C[i * 3] * tmp[j];
+      Slater_inv[i * Dim + j] -= C[i * 3 + 1] * tmp[Dim + j];
+      Slater_inv[i * Dim + j] -= C[i * 3 + 2] * tmp[2 * Dim + j];
+    }
   }
 
 #ifdef DEBUG2
diff --git a/tests/fnu_test_h5.cpp b/tests/fnu_test_h5.cpp
index f0a171a..3f9763a 100644
--- a/tests/fnu_test_h5.cpp
+++ b/tests/fnu_test_h5.cpp
@@ -2,10 +2,10 @@
 #include "hdf5/serial/hdf5.h"
 
 #include "Helpers.hpp"
+#include "SMWB.hpp"
 #include "SM_Maponi.hpp"
 #include "SM_Standard.hpp"
 #include "Woodbury.hpp"
-#include "SMWB.hpp"
 #include <fstream>
 #include <vector>
 
@@ -59,10 +59,10 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
 
   double *u = new double[nupdates * dim];
 
-  /* Test */
-  #ifdef DEBUG2
+/* Test */
+#ifdef DEBUG2
   showMatrix(slater_inverse, dim, "OLD Inverse");
-  #endif
+#endif
 
   // Transform replacement updates in 'updates[]' into additive updates in 'u[]'
   for (j = 0; j < nupdates; j++) {
@@ -74,70 +74,80 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
     }
   }
 
-  #ifdef DEBUG2
+#ifdef DEBUG2
   showMatrix(slater_matrix, dim, "OLD Slater");
   showMatrix(u, dim, "Updates");
-  #endif
+#endif
 
-  #ifdef PERF
-  #ifdef DEBUG1
+#ifdef PERF
+#ifdef DEBUG1
   std::cerr << "# of reps. = " << repetition_number << std::endl;
-  #endif // DEBUG1
+#endif // DEBUG1
   double *slater_inverse_nonpersistent = new double[dim * dim];
   if (version == "sm1") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       SM1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     }
   } else if (version == "sm2") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       SM2(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     }
   } else if (version == "sm3") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       SM3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     }
   } else if (version == "sm4") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       SM4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     }
   } else if (version == "wb2") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       WB2(slater_inverse_nonpersistent, dim, u, col_update_index);
     }
   } else if (version == "wb3") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       WB3(slater_inverse_nonpersistent, dim, u, col_update_index);
     }
   } else if (version == "smwb1") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       SMWB1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     }
   } else if (version == "smwb4") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                  dim * dim * sizeof(double));
       SMWB4(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     }
-  #ifdef MKL
+#ifdef MKL
   } else if (version == "lapack") {
     for (unsigned int i = 0; i < repetition_number; i++) {
-      std::memcpy(slater_inverse_nonpersistent, slater_matrix, dim * dim * sizeof(double));
+      std::memcpy(slater_inverse_nonpersistent, slater_matrix,
+                  dim * dim * sizeof(double));
       inverse(slater_inverse_nonpersistent, dim);
     }
-  #endif // MKL
+#endif // MKL
   } else {
     std::cerr << "Unknown version " << version << std::endl;
     exit(1);
   }
-  std::memcpy(slater_inverse, slater_inverse_nonpersistent, dim * dim * sizeof(double));
+  std::memcpy(slater_inverse, slater_inverse_nonpersistent,
+              dim * dim * sizeof(double));
   delete[] slater_inverse_nonpersistent;
-  #else //  No performance measurements repetition
+#else //  No performance measurements repetition
   if (version == "maponia3") {
     MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "maponia3s") {
@@ -156,27 +166,27 @@ 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
+#ifdef MKL
   } else if (version == "lapack") {
     memcpy(slater_inverse, slater_matrix, dim * dim * sizeof(double));
     inverse(slater_inverse, dim);
-  #endif // MKL
+#endif // MKL
   } else {
     std::cerr << "Unknown version " << version << std::endl;
     exit(1);
   }
-  #endif // PERF
+#endif // PERF
 
-  #ifdef DEBUG2
+#ifdef DEBUG2
   showMatrix(slater_matrix, dim, "NEW Slater");
   showMatrix(slater_inverse, dim, "NEW Inverse");
-  #endif
+#endif
 
   double *res = new double[dim * dim]{0};
   matMul(slater_matrix, slater_inverse, res, dim);
@@ -184,14 +194,14 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
   double res_max = residual_max(res, dim);
   double res2 = residual_frobenius2(res, dim);
 
-  #ifdef RESIDUAL
+#ifdef RESIDUAL
   std::cout << "Residual = " << version << " " << cycle << " " << res_max << " "
             << res2 << std::endl;
-  #endif
+#endif
 
-  #ifdef DEBUG2
+#ifdef DEBUG2
   showMatrix(res, dim, "Result");
-  #endif
+#endif
 
   delete[] res, updates, u, col_update_index, slater_matrix, slater_inverse;
 
@@ -199,7 +209,7 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
 }
 
 int main(int argc, char **argv) {
-  #ifdef PERF
+#ifdef PERF
   if (argc != 5) {
     std::cerr << "Execute from within 'datasets/'" << std::endl;
     std::cerr
@@ -207,38 +217,38 @@ int main(int argc, char **argv) {
         << std::endl;
     return 1;
   }
-  #else
+#else
   if (argc != 4) {
     std::cerr << "Execute from within 'datasets/'" << std::endl;
-    std::cerr
-        << "usage: test_h5 <version> <cycle file> <tolerance>"
-        << std::endl;
+    std::cerr << "usage: test_h5 <version> <cycle file> <tolerance>"
+              << std::endl;
     return 1;
   }
-  #endif    
+#endif
   std::string version(argv[1]);
   std::string cyclefile_name(argv[2]);
   std::ifstream cyclefile(cyclefile_name);
   std::vector<int> cycles;
   unsigned int cycle;
-  while (cyclefile >> cycle) cycles.push_back(cycle);
+  while (cyclefile >> cycle)
+    cycles.push_back(cycle);
   double tolerance = std::stod(argv[3]);
   H5File file(FILE_NAME, H5F_ACC_RDONLY);
 
-  #ifdef PERF
+#ifdef PERF
   repetition_number = std::stoi(argv[4]);
-  #endif
+#endif
 
   bool ok;
-  for (auto & cycle : cycles) {
+  for (auto &cycle : cycles) {
     ok = test_cycle(file, cycle, version, tolerance);
-    #ifdef STATUS
+#ifdef STATUS
     if (ok) {
       std::cerr << "ok -- cycle " << std::to_string(cycle) << std::endl;
     } else {
       std::cerr << "failed -- cycle " << std::to_string(cycle) << std::endl;
     }
-    #endif
+#endif
   }
 
   return ok;
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index bb83b7f..5763125 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -2,10 +2,10 @@
 #include "hdf5/serial/hdf5.h"
 
 #include "Helpers.hpp"
+#include "SMWB.hpp"
 #include "SM_Maponi.hpp"
 #include "SM_Standard.hpp"
 #include "Woodbury.hpp"
-#include "SMWB.hpp"
 
 #define PERF
 
@@ -82,11 +82,14 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
   std::cout << "# of reps. = " << repetition_number << std::endl;
   double *slater_inverse_nonpersistent = new double[dim * dim];
   for (unsigned int i = 0; i < repetition_number; i++) {
-    std::memcpy(slater_inverse_nonpersistent, slater_inverse, dim * dim * sizeof(double));
+    std::memcpy(slater_inverse_nonpersistent, slater_inverse,
+                dim * dim * sizeof(double));
     if (version == "maponia3") {
-      MaponiA3(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
+      MaponiA3(slater_inverse_nonpersistent, dim, nupdates, u,
+               col_update_index);
     } else if (version == "maponia3s") {
-      MaponiA3S(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
+      MaponiA3S(slater_inverse_nonpersistent, dim, nupdates, u,
+                col_update_index);
     } else if (version == "sm1") {
       SM1(slater_inverse_nonpersistent, dim, nupdates, u, col_update_index);
     } else if (version == "sm2") {
@@ -101,15 +104,18 @@ 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
     } else if (version == "lapack") {
-      memcpy(slater_inverse_nonpersistent, slater_matrix, dim * dim * sizeof(double));
+      memcpy(slater_inverse_nonpersistent, slater_matrix,
+             dim * dim * sizeof(double));
       inverse(slater_inverse_nonpersistent, dim);
 #endif // MKL
     } else {
@@ -117,7 +123,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
       exit(1);
     }
   }
-  std::memcpy(slater_inverse, slater_inverse_nonpersistent, dim * dim * sizeof(double));
+  std::memcpy(slater_inverse, slater_inverse_nonpersistent,
+              dim * dim * sizeof(double));
   delete[] slater_inverse_nonpersistent;
 #else
   if (version == "maponia3") {
@@ -138,10 +145,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
diff --git a/tests/vfc_test_h5.cpp b/tests/vfc_test_h5.cpp
index f30f32d..e02355c 100644
--- a/tests/vfc_test_h5.cpp
+++ b/tests/vfc_test_h5.cpp
@@ -166,7 +166,7 @@ int main(int argc, char **argv) {
   bool ok;
   for (int i = 0; i < cycles_list.size(); i++) {
 
-    for(const auto& version: versions) {
+    for (const auto &version : versions) {
 
       ok = test_cycle(file, cycles_list[i], version, &probes);
       if (ok) {
@@ -174,9 +174,7 @@ int main(int argc, char **argv) {
       } else {
         std::cerr << "failed -- cycle " << std::to_string(i) << std::endl;
       }
-
     }
-
   }
 
   vfc_dump_probes(&probes);
diff --git a/tools/restyle.sh b/tools/restyle.sh
index b16069b..0743ab0 100755
--- a/tools/restyle.sh
+++ b/tools/restyle.sh
@@ -11,5 +11,5 @@ fi
 
 for ext in c cc cpp h hpp
 do
-    find $SMROOT -type f -iname "*.${ext}" -exec echo "$FORMATER $STYLE" {} \;
+    find $SMROOT -type f -iname "*.${ext}" -exec $FORMATER $STYLE {} \;
 done