From 288bc83e1952f773cd17c2d7ea4dde4721af181e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <mail@fmgjcoppens.nl>
Date: Fri, 4 Jun 2021 16:50:49 +0200
Subject: [PATCH 1/6] Created Woodbery kernel infrastructure and started
 wrinting Woodbery 2x2 kernel.

---
 Makefile                                |  4 +-
 include/{SM_Helpers.hpp => Helpers.hpp} |  2 +-
 include/SMWB.hpp                        |  4 ++
 include/Woodbury.hpp                    |  7 +++
 src/{SM_Helpers.cpp => Helpers.cpp}     |  2 +-
 src/SMWB.cpp                            | 21 ++++++++
 src/SM_Maponi.cpp                       |  2 +-
 src/SM_Standard.cpp                     |  2 +-
 src/Woodbury.cpp                        | 69 +++++++++++++++++++++++++
 src/finterface_mod.f90                  |  9 +++-
 tests/cMaponiA3_test_3x3_3.cpp          |  2 +-
 tests/test_h5.cpp                       |  5 +-
 tests/vfc_test_h5.cpp                   |  2 +-
 tools/restyle.sh                        |  4 +-
 14 files changed, 124 insertions(+), 11 deletions(-)
 rename include/{SM_Helpers.hpp => Helpers.hpp} (99%)
 create mode 100644 include/SMWB.hpp
 create mode 100644 include/Woodbury.hpp
 rename src/{SM_Helpers.cpp => Helpers.cpp} (98%)
 create mode 100644 src/SMWB.cpp
 create mode 100644 src/Woodbury.cpp

diff --git a/Makefile b/Makefile
index d69bd40..cb7ad9b 100644
--- a/Makefile
+++ b/Makefile
@@ -42,7 +42,9 @@ FFLAGS = $(CXXFLAGS)
 INCLUDE = -I $(INC_DIR)/
 DEPS_CXX = $(OBJ_DIR)/SM_Maponi.o \
 		   $(OBJ_DIR)/SM_Standard.o \
-		   $(OBJ_DIR)/SM_Helpers.o
+		   $(OBJ_DIR)/Woodbury.o \
+		   $(OBJ_DIR)/SMWB.o \
+		   $(OBJ_DIR)/Helpers.o
 DEPS_F = $(DEPS_CXX) \
 		 $(OBJ_DIR)/finterface_mod.o \
 		 $(OBJ_DIR)/helpers_mod.o
diff --git a/include/SM_Helpers.hpp b/include/Helpers.hpp
similarity index 99%
rename from include/SM_Helpers.hpp
rename to include/Helpers.hpp
index b569d18..015c2ee 100644
--- a/include/SM_Helpers.hpp
+++ b/include/Helpers.hpp
@@ -1,4 +1,4 @@
-// SM_Helpers.hpp
+// Helpers.hpp
 // Some usefull helper functions to support the Maponi algorithm.
 #include <cmath>
 #include <cstring>
diff --git a/include/SMWB.hpp b/include/SMWB.hpp
new file mode 100644
index 0000000..0de1a49
--- /dev/null
+++ b/include/SMWB.hpp
@@ -0,0 +1,4 @@
+// Sherman-Morrison-Woodbury kernel 1
+// WB2, WB3, SM2 mixing scheme 1
+void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
+         double *Updates, unsigned int *Updates_index);
diff --git a/include/Woodbury.hpp b/include/Woodbury.hpp
new file mode 100644
index 0000000..001d912
--- /dev/null
+++ b/include/Woodbury.hpp
@@ -0,0 +1,7 @@
+// Woodbury 2x2 kernel
+void WB2(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index);
+
+// Woodbury 3x3 kernel
+void WB3(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index);
diff --git a/src/SM_Helpers.cpp b/src/Helpers.cpp
similarity index 98%
rename from src/SM_Helpers.cpp
rename to src/Helpers.cpp
index 4d084c0..12ce105 100644
--- a/src/SM_Helpers.cpp
+++ b/src/Helpers.cpp
@@ -1,4 +1,4 @@
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 
 // Set common break-down threshold
 double threshold() {
diff --git a/src/SMWB.cpp b/src/SMWB.cpp
new file mode 100644
index 0000000..0590ff5
--- /dev/null
+++ b/src/SMWB.cpp
@@ -0,0 +1,21 @@
+#include "SMWB.hpp"
+#include "SM_Standard.hpp"
+#include "Woodbury.hpp"
+#include "Helpers.hpp"
+
+// Sherman-Morrison-Woodbury kernel 1
+// WB2, WB3, SM2 mixing scheme 1
+void SMWB1(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;
+           WB3(Slater_inv, Dim, Updates, Updates_index);
+           WB2(Slater_inv, Dim, Updates, Updates_index);
+           SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
+         }
+
+extern "C" {
+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);
+}
+}
diff --git a/src/SM_Maponi.cpp b/src/SM_Maponi.cpp
index 7b65fde..65afefa 100644
--- a/src/SM_Maponi.cpp
+++ b/src/SM_Maponi.cpp
@@ -2,7 +2,7 @@
 // Algorithm 3 from P. Maponi,
 // p. 283, doi:10.1016/j.laa.2006.07.007
 #include "SM_Maponi.hpp"
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 
 // #define DEBUG
 
diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp
index 91b69c9..c68dd81 100644
--- a/src/SM_Standard.cpp
+++ b/src/SM_Standard.cpp
@@ -1,7 +1,7 @@
 // SM-Standard.cpp
 // Standard Sherman Morrison with multiple updates
 #include "SM_Standard.hpp"
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 
 // Naïve Sherman Morrison
 void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
new file mode 100644
index 0000000..21b9d9d
--- /dev/null
+++ b/src/Woodbury.cpp
@@ -0,0 +1,69 @@
+// Woodbury.cpp
+// Woodbury 2x2 and 3x3 kernels
+//
+// Woodbury matrix identity:
+// (S + U * V)^{-1} = S^{-1} - S^{-1} * U * B^{-1} * V * S^{-1}
+// B := 1 + V * C,  2 x 2
+// C := S^{-1} * U, dim x 2
+#include "Woodbury.hpp"
+#include "Helpers.hpp"
+
+// Woodbury 2x2 kernel: 
+void WB2(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index) {
+  std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
+
+  // Construct V from Updates_index
+  unsigned int *V = new unsigned int[2 * Dim]{0};
+  for (unsigned int i = 0; i < Dim; i++) {
+    if (i == Updates_index[0] - 1) V[i] = 1;
+    if (Dim + i == Updates_index[1] - 1) V[Dim + i] = 1;
+  }
+
+  // Compute B from U, V and Slater_inv
+  double B[4]{0};
+  double Binv[4]{0};
+  double C[4]{0};
+  double D[2*Dim]{0};
+    // Compute C
+    for (unsigned int i = 0; i < Dim; i++) {
+      for (unsigned int j = 0; j < 2; j++) {
+        for (unsigned int k = 0; k < Dim; k++) {
+          C[i * Dim + j] += Slater_inv[i * Dim + k] * Updates[k * Dim + j];
+        }
+      }
+    }
+    // Compute B
+    for (unsigned int i = 0; i < 2; i++) {
+      for (unsigned int j = 0; j < 2; j++) {
+        for (unsigned int k = 0; k < Dim; k++) {
+          B[i * Dim + j] += (i == j) + V[i * Dim + k] * C[k * Dim + j];
+        }
+      }
+    }
+
+  // Invert B with explicit formula for 2x2 inversion
+  double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]);
+  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 (S + U * V)^{-1} with Woobury identity
+  
+}
+
+// Woodbury 3x3 kernel
+void WB3(double *Slater_inv, unsigned int Dim,
+         double *Updates, unsigned int *Updates_index) {
+           std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
+
+           // Construct V from Updates_index
+
+           // Compute B from U, V and Slater_inv
+
+           // Invert B with explicit formula for 3x3 inversion
+
+           // Compute (S + U * V)^{-1} with Woobury identity
+           
+         }
diff --git a/src/finterface_mod.f90 b/src/finterface_mod.f90
index bcf072c..09de1e8 100644
--- a/src/finterface_mod.f90
+++ b/src/finterface_mod.f90
@@ -43,5 +43,12 @@ module Sherman_Morrison
             real(c_double), dimension(:,:), allocatable, intent(in) :: Updates
             real(c_double), dimension(:,:), allocatable, intent(in out) :: Slater_inv
         end subroutine SM4
-    end interface
+        subroutine SMWB1(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SMWB1_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 SMWB1
+        end interface
 end module Sherman_Morrison
diff --git a/tests/cMaponiA3_test_3x3_3.cpp b/tests/cMaponiA3_test_3x3_3.cpp
index bcd363a..7f6b2af 100644
--- a/tests/cMaponiA3_test_3x3_3.cpp
+++ b/tests/cMaponiA3_test_3x3_3.cpp
@@ -1,5 +1,5 @@
 // main.cpp
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 #include "SM_Maponi.hpp"
 
 int main() {
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index fb0ceda..f750f10 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -1,9 +1,10 @@
 #include "hdf5/serial/H5Cpp.h"
 #include "hdf5/serial/hdf5.h"
 
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 #include "SM_Maponi.hpp"
 #include "SM_Standard.hpp"
+#include "SMWB.hpp"
 
 using namespace H5;
 // #define DEBUG
@@ -81,6 +82,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
     SM3(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "sm4") {
     SM4(slater_inverse, dim, nupdates, u, col_update_index);
+  } else if (version == "smwb1") {
+    SMWB1(slater_inverse, dim, nupdates, u, col_update_index);
 #ifdef MKL
   } else if (version == "lapack") {
     memcpy(slater_inverse, slater_matrix, dim * dim * sizeof(double));
diff --git a/tests/vfc_test_h5.cpp b/tests/vfc_test_h5.cpp
index 920e069..d7f3c8b 100644
--- a/tests/vfc_test_h5.cpp
+++ b/tests/vfc_test_h5.cpp
@@ -10,7 +10,7 @@
 #include <sstream>
 #include <vector>
 
-#include "SM_Helpers.hpp"
+#include "Helpers.hpp"
 #include "SM_Maponi.hpp"
 #include "SM_Standard.hpp"
 #include "vfc_probe.h"
diff --git a/tools/restyle.sh b/tools/restyle.sh
index ba07f81..b16069b 100755
--- a/tools/restyle.sh
+++ b/tools/restyle.sh
@@ -9,7 +9,7 @@ then
 	exit 1
 fi
 
-for ext in c cpp h hpp
+for ext in c cc cpp h hpp
 do
     find $SMROOT -type f -iname "*.${ext}" -exec echo "$FORMATER $STYLE" {} \;
-done
\ No newline at end of file
+done

From efe96cbeeaccce41a3c37c70e709f71f166f9d91 Mon Sep 17 00:00:00 2001
From: Francois Coppens <francois.coppens@uvsq.fr>
Date: Thu, 10 Jun 2021 08:46:40 +0200
Subject: [PATCH 2/6] First implementation of Woodbury 2x2 and 3x3 kernels.

---
 include/Helpers.hpp  |  51 ++++++++++++++++
 include/Woodbury.hpp |   4 +-
 smvars.sh            |   4 ++
 src/SMWB.cpp         |  13 ++--
 src/Woodbury.cpp     | 139 +++++++++++++++++++++++++++++++------------
 todo.txt             |  10 ----
 6 files changed, 165 insertions(+), 56 deletions(-)
 delete mode 100644 todo.txt

diff --git a/include/Helpers.hpp b/include/Helpers.hpp
index 015c2ee..2384db1 100644
--- a/include/Helpers.hpp
+++ b/include/Helpers.hpp
@@ -56,6 +56,45 @@ void showMatrix(T *matrix, unsigned int M, std::string name) {
   std::cout << std::endl;
 }
 
+template <typename T>
+void showMatrix2(T *matrix, unsigned int M, unsigned int N, std::string name) {
+  std::cout.precision(17);
+  std::cout << name << " = [" << std::endl;
+  for (unsigned int i = 0; i < M; i++) {
+    std::cout << "[";
+    for (unsigned int j = 0; j < N; j++) {
+      if (matrix[i * N + j] >= 0) {
+        std::cout << "  " << matrix[i * N + j] << ",";
+      } else {
+        std::cout << " " << matrix[i * N + j] << ",";
+      }
+    }
+    std::cout << " ]," << std::endl;
+  }
+  std::cout << "]" << std::endl;
+  std::cout << std::endl;
+}
+
+
+template <typename T>
+void showMatrixNS(T *matrix, unsigned int M, unsigned int N, std::string name) {
+  std::cout.precision(17);
+  std::cout << name << " = [" << std::endl;
+  for (unsigned int i = 0; i < M; i++) {
+    std::cout << "[";
+    for (unsigned int j = 0; j < N; j++) {
+      if (matrix[i * N + j] >= 0) {
+        std::cout << "  " << matrix[i * N + j] << ",";
+      } else {
+        std::cout << " " << matrix[i * N + j] << ",";
+      }
+    }
+    std::cout << " ]," << std::endl;
+  }
+  std::cout << "]" << std::endl;
+  std::cout << std::endl;
+}
+
 template <typename T> T *transpose(T *A, unsigned int M) {
   T *B = new T[M * M];
   for (unsigned int i = 0; i < M; i++) {
@@ -77,6 +116,18 @@ template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
   }
 }
 
+template <typename T1, typename T2, typename T3>
+static inline 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++) {
+        C[i * P + j] += A[i * N + k] * B[k * P + j];
+      }
+    }
+  }
+}
+
 template <typename T1, typename T2>
 T1 *outProd(T1 *vec1, T2 *vec2, unsigned int M) {
   T1 *C = new T1[M * M];
diff --git a/include/Woodbury.hpp b/include/Woodbury.hpp
index 001d912..c03e509 100644
--- a/include/Woodbury.hpp
+++ b/include/Woodbury.hpp
@@ -1,7 +1,7 @@
 // Woodbury 2x2 kernel
-void WB2(double *Slater_inv, unsigned int Dim,
+bool WB2(double *Slater_inv, unsigned int Dim,
          double *Updates, unsigned int *Updates_index);
 
 // Woodbury 3x3 kernel
-void WB3(double *Slater_inv, unsigned int Dim,
+bool WB3(double *Slater_inv, unsigned int Dim,
          double *Updates, unsigned int *Updates_index);
diff --git a/smvars.sh b/smvars.sh
index 6e846ee..2e9e868 100644
--- a/smvars.sh
+++ b/smvars.sh
@@ -3,6 +3,10 @@ unset THRESHOLD
 unset MKL
 ENV=$1
 
+echo
+echo "Sherman-Morrison-Woodbury parameters"
+echo "------------------------------------"
+
 ## Set Sherman-Morrison root dir
 PWD=$(pwd)
 SRCDIR=$(dirname $BASH_SOURCE)
diff --git a/src/SMWB.cpp b/src/SMWB.cpp
index 0590ff5..670920b 100644
--- a/src/SMWB.cpp
+++ b/src/SMWB.cpp
@@ -5,12 +5,15 @@
 
 // Sherman-Morrison-Woodbury kernel 1
 // WB2, WB3, SM2 mixing scheme 1
-void SMWB1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-         double *Updates, unsigned int *Updates_index) {
+void SMWB1(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;
-           WB3(Slater_inv, Dim, Updates, Updates_index);
-           WB2(Slater_inv, Dim, Updates, Updates_index);
-           SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
+
+           bool ok;
+           ok = WB2(Slater_inv, Dim, Updates, Updates_index);
+           if (!ok) {
+             std::cerr << "Woodbury kernel failed!" << std::endl;
+             SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
+           }
          }
 
 extern "C" {
diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index 21b9d9d..980464b 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -5,65 +5,126 @@
 // (S + U * V)^{-1} = S^{-1} - S^{-1} * U * B^{-1} * V * S^{-1}
 // B := 1 + V * C,  2 x 2
 // C := S^{-1} * U, dim x 2
+// D := V * S^{-1}, 2 x dim
+//
+// All matrices are stored in row-major order
+
 #include "Woodbury.hpp"
 #include "Helpers.hpp"
 
-// Woodbury 2x2 kernel: 
-void WB2(double *Slater_inv, unsigned int Dim,
-         double *Updates, unsigned int *Updates_index) {
+// Woodbury 2x2 kernel:
+bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
+         unsigned int *Updates_index) {
   std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
 
   // Construct V from Updates_index
-  unsigned int *V = new unsigned int[2 * Dim]{0};
-  for (unsigned int i = 0; i < Dim; i++) {
-    if (i == Updates_index[0] - 1) V[i] = 1;
-    if (Dim + i == Updates_index[1] - 1) V[Dim + i] = 1;
-  }
+  unsigned int V[2 * Dim]{0}; // 2 x Dim matrix stored in row-major order
+  V[Updates_index[0] - 1] = 1;
+  V[Dim + Updates_index[1] - 1] = 1;
 
-  // Compute B from U, V and Slater_inv
+  // Compute C
+  double C[2 * Dim]{0};
+  matMul2(Slater_inv, Updates, C, Dim, Dim, 2);
+  // Compute B
   double B[4]{0};
-  double Binv[4]{0};
-  double C[4]{0};
-  double D[2*Dim]{0};
-    // Compute C
-    for (unsigned int i = 0; i < Dim; i++) {
-      for (unsigned int j = 0; j < 2; j++) {
-        for (unsigned int k = 0; k < Dim; k++) {
-          C[i * Dim + j] += Slater_inv[i * Dim + k] * Updates[k * Dim + j];
-        }
-      }
-    }
-    // Compute B
-    for (unsigned int i = 0; i < 2; i++) {
-      for (unsigned int j = 0; j < 2; j++) {
-        for (unsigned int k = 0; k < Dim; k++) {
-          B[i * Dim + j] += (i == j) + V[i * Dim + k] * C[k * Dim + j];
-        }
-      }
-    }
+  matMul2(V, C, B, 2, Dim, 2);
+  // Compute 1 + B
+  B[0] += 1;
+  B[3] += 1;
 
-  // Invert B with explicit formula for 2x2 inversion
+  // Invert 1 + B with explicit formula for 2x2 inversion
   double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]);
+  double Binv[4]{0};
   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 = B[0] * B[3] - B[1] * B[2];
+  if (std::fabs(det) < threshold()) {
+    std::cerr << "Determinant approached 0!" << std::endl;
+    return false;
+  }
+
   // Compute (S + U * V)^{-1} with Woobury identity
-  
+  double D[2 * Dim]{0};
+  matMul2(V, Slater_inv, D, 2, Dim, Dim);
+  double tmp[2 * Dim]{0};
+  matMul2(Binv, D, tmp, 2, 2, Dim);
+  double tmp2[Dim * Dim]{0};
+  matMul2(C, tmp, tmp2, Dim, 2, Dim);
+  for (unsigned int i = 0; i < Dim * Dim; i++) {
+    Slater_inv[i] -= tmp2[i];
+  }
+  return true;
 }
 
 // Woodbury 3x3 kernel
-void WB3(double *Slater_inv, unsigned int Dim,
-         double *Updates, unsigned int *Updates_index) {
-           std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
+bool WB3(double *Slater_inv, unsigned int Dim, double *Updates,
+         unsigned int *Updates_index) {
+  std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
+  
+  // Construct V from Updates_index
+  unsigned int V[3 * Dim]{0}; // 2 x Dim matrix stored in row-major order
+  V[Updates_index[0] - 1] = 1;
+  V[Dim + Updates_index[1] - 1] = 1;
+  V[2 * Dim + Updates_index[2] - 1] = 1;
 
-           // Construct V from Updates_index
+  // Compute C
+  double C[3 * Dim]{0};
+  matMul2(Slater_inv, Updates, C, Dim, Dim, 3);
+  // Compute B
+  double B[9]{0};
+  matMul2(V, C, B, 3, Dim, 3);
+  // Compute 1 + B
+  B[0] += 1;
+  B[4] += 1;
+  B[8] += 1;
 
-           // Compute B from U, V and Slater_inv
+  double Binv[9];
+  Binv[0] = B[4] * B[8] - B[5] * B[7];
+  Binv[3] = B[5] * B[6] - B[3] * B[8];
+  Binv[6] = B[3] * B[7] - B[4] * B[6];
 
-           // Invert B with explicit formula for 3x3 inversion
+  Binv[1] = B[2] * B[7] - B[1] * B[8];
+  Binv[4] = B[0] * B[8] - B[2] * B[6];
+  Binv[7] = B[1] * B[6] - B[0] * B[7];
 
-           // Compute (S + U * V)^{-1} with Woobury identity
-           
-         }
+  Binv[2] = B[1] * B[5] - B[2] * B[4];
+  Binv[5] = B[2] * B[3] - B[0] * B[5];
+  Binv[8] = B[0] * B[4] - B[1] * B[3];
+
+  // Check if determinant of inverted matrix is not zero
+  // If so, exigt and return false.
+  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]);
+  if (std::fabs(det) < threshold()) {
+    std::cerr << "Determinant approached 0!" << std::endl;
+    return false;
+  }
+
+  // Compute (S + U * V)^{-1} with Woobury identity
+  double D[3 * Dim]{0};
+  matMul2(V, Slater_inv, D, 3, Dim, Dim);
+  double tmp[3 * Dim]{0};
+  matMul2(Binv, D, tmp, 3, 3, Dim);
+  double tmp2[Dim * Dim]{0};
+  matMul2(C, tmp, tmp2, Dim, 3, Dim);
+  for (unsigned int i = 0; i < Dim * Dim; i++) {
+    Slater_inv[i] -= tmp2[i];
+  }
+  return true;
+}
+
+extern "C" {
+bool WB2_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
+           unsigned int **Updates_index) {
+  WB2(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
+}
+bool WB3_f(double **linSlater_inv, unsigned int *Dim, double **linUpdates,
+           unsigned int **Updates_index) {
+  WB3(*linSlater_inv, *Dim, *linUpdates, *Updates_index);
+}
+}
diff --git a/todo.txt b/todo.txt
deleted file mode 100644
index 5a579c0..0000000
--- a/todo.txt
+++ /dev/null
@@ -1,10 +0,0 @@
--=={TODO}==-
-
-* Looking at the eigenvalues of the intermediate matrix after a problematic update is
-* plotting the ratio of the moduli of the maximum and the minimum eigenvalues of
-
-* Use valgrind to find the problem with MaponiA3S
-* Contact Claudia if she would be interested in sharing datasets from Champ to test with SM
-
-* Make a representative selection of update cycles for presentatoin puroses
-

From 61844da5d347e9e7c07e7c50826a74d50ab9d524 Mon Sep 17 00:00:00 2001
From: Francois Coppens <francois.coppens@uvsq.fr>
Date: Thu, 10 Jun 2021 11:25:31 +0200
Subject: [PATCH 3/6] Corrected array initialisations.

---
 src/Woodbury.cpp | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index 980464b..1c92a68 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -18,15 +18,15 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
 
   // Construct V from Updates_index
-  unsigned int V[2 * Dim]{0}; // 2 x Dim matrix stored in row-major order
+  unsigned int V[2 * Dim]; // 2 x Dim matrix stored in row-major order
   V[Updates_index[0] - 1] = 1;
   V[Dim + Updates_index[1] - 1] = 1;
 
   // Compute C
-  double C[2 * Dim]{0};
+  double C[2 * Dim];
   matMul2(Slater_inv, Updates, C, Dim, Dim, 2);
   // Compute B
-  double B[4]{0};
+  double B[4];
   matMul2(V, C, B, 2, Dim, 2);
   // Compute 1 + B
   B[0] += 1;
@@ -34,7 +34,7 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
 
   // Invert 1 + B with explicit formula for 2x2 inversion
   double idet = 1.0 / (B[0] * B[3] - B[1] * B[2]);
-  double Binv[4]{0};
+  double Binv[4];
   Binv[0] = idet * B[3];
   Binv[1] = -1.0 * idet * B[1];
   Binv[2] = -1.0 * idet * B[2];
@@ -48,11 +48,11 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   }
 
   // Compute (S + U * V)^{-1} with Woobury identity
-  double D[2 * Dim]{0};
+  double D[2 * Dim];
   matMul2(V, Slater_inv, D, 2, Dim, Dim);
-  double tmp[2 * Dim]{0};
+  double tmp[2 * Dim];
   matMul2(Binv, D, tmp, 2, 2, Dim);
-  double tmp2[Dim * Dim]{0};
+  double tmp2[Dim * Dim];
   matMul2(C, tmp, tmp2, Dim, 2, Dim);
   for (unsigned int i = 0; i < Dim * Dim; i++) {
     Slater_inv[i] -= tmp2[i];
@@ -66,16 +66,16 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates,
   std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
   
   // Construct V from Updates_index
-  unsigned int V[3 * Dim]{0}; // 2 x Dim matrix stored in row-major order
+  unsigned int V[3 * Dim]; // 2 x Dim matrix stored in row-major order
   V[Updates_index[0] - 1] = 1;
   V[Dim + Updates_index[1] - 1] = 1;
   V[2 * Dim + Updates_index[2] - 1] = 1;
 
   // Compute C
-  double C[3 * Dim]{0};
+  double C[3 * Dim];
   matMul2(Slater_inv, Updates, C, Dim, Dim, 3);
   // Compute B
-  double B[9]{0};
+  double B[9];
   matMul2(V, C, B, 3, Dim, 3);
   // Compute 1 + B
   B[0] += 1;
@@ -106,11 +106,11 @@ bool WB3(double *Slater_inv, unsigned int Dim, double *Updates,
   }
 
   // Compute (S + U * V)^{-1} with Woobury identity
-  double D[3 * Dim]{0};
+  double D[3 * Dim];
   matMul2(V, Slater_inv, D, 3, Dim, Dim);
-  double tmp[3 * Dim]{0};
+  double tmp[3 * Dim];
   matMul2(Binv, D, tmp, 3, 3, Dim);
-  double tmp2[Dim * Dim]{0};
+  double tmp2[Dim * Dim];
   matMul2(C, tmp, tmp2, Dim, 3, Dim);
   for (unsigned int i = 0; i < Dim * Dim; i++) {
     Slater_inv[i] -= tmp2[i];

From 573947fe2df601d84e983b9a38b362864edc93d7 Mon Sep 17 00:00:00 2001
From: Francois Coppens <francois.coppens@uvsq.fr>
Date: Fri, 11 Jun 2021 08:46:39 +0200
Subject: [PATCH 4/6] Woodbury debugging...

---
 include/Helpers.hpp |  2 +-
 src/Woodbury.cpp    | 14 ++++++++++++++
 tests/test_h5.cpp   |  2 ++
 3 files changed, 17 insertions(+), 1 deletion(-)

diff --git a/include/Helpers.hpp b/include/Helpers.hpp
index 2384db1..cd6fdbb 100644
--- a/include/Helpers.hpp
+++ b/include/Helpers.hpp
@@ -117,7 +117,7 @@ template <typename T> void matMul(T *A, T *B, T *C, unsigned int M) {
 }
 
 template <typename T1, typename T2, typename T3>
-static inline void matMul2(T1 *A, T2 *B, T3 *C, unsigned int M, unsigned int N, unsigned int P) {
+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;
diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index 1c92a68..3b8e26a 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -19,12 +19,26 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
 
   // Construct V from Updates_index
   unsigned int V[2 * Dim]; // 2 x Dim matrix stored in row-major order
+  std::memset(V, 0, 2 * Dim * sizeof(unsigned int));
   V[Updates_index[0] - 1] = 1;
   V[Dim + Updates_index[1] - 1] = 1;
 
+  showMatrix2(Slater_inv, Dim, Dim, "Slater_inv");
+  showMatrix2(Updates, 2, Dim, "Updates");
+  showMatrix2(Updates_index, 1, 2, "Updates_index");
+  showMatrix2(V, 2, Dim, "V");
+
   // Compute C
   double C[2 * Dim];
   matMul2(Slater_inv, Updates, C, Dim, Dim, 2);
+
+  int A[6] = {1,2,3,4,5,6};
+  int Y[12] = {7,8,9,10,11,12,13,14,15,16,17,18};
+  int Z[8] = {0,0,0,0,0,0,0,0};
+  matMul2(A,Y,Z,2,3,4);
+  showMatrix2(Z, 2,4,"Z");
+
+  showMatrix2(C, 2, Dim, "C = S_inv * Updates");
   // Compute B
   double B[4];
   matMul2(V, C, B, 2, Dim, 2);
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index f750f10..9b49a1c 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -52,6 +52,8 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
   showMatrix(slater_inverse, dim, "OLD Inverse");
 #endif
 
+  showMatrix2(slater_matrix, dim, dim, "Slater");
+  
   // Transform replacement updates in 'updates[]' into additive updates in 'u[]'
   for (j = 0; j < nupdates; j++) {
     for (i = 0; i < dim; i++) {

From b6efc972338b1bd8e75a17b3a931b52737a5328d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <mail@fmgjcoppens.nl>
Date: Mon, 14 Jun 2021 15:05:39 +0200
Subject: [PATCH 5/6] Woodbury 2x2 kernel fixed and working.

---
 src/Woodbury.cpp | 62 ++++++++++++++++++++++++++----------------------
 1 file changed, 33 insertions(+), 29 deletions(-)

diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index 3b8e26a..43eb301 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -1,13 +1,15 @@
-// Woodbury.cpp
-// Woodbury 2x2 and 3x3 kernels
-//
-// Woodbury matrix identity:
-// (S + U * V)^{-1} = S^{-1} - S^{-1} * U * B^{-1} * V * S^{-1}
-// B := 1 + V * C,  2 x 2
-// C := S^{-1} * U, dim x 2
-// D := V * S^{-1}, 2 x dim
-//
-// All matrices are stored in row-major order
+/*  Woodbury.cpp
+    Woodbury 2x2 and 3x3 kernels
+
+    Woodbury matrix identity:
+
+    (S + U V)^{-1}  =  S^{-1} - C B^{-1} D
+
+    C := S^{-1} U,    dim x 2
+    B := 1 + V C,     2 x 2
+    D := V S^{-1},    2 x dim
+
+    All matrices are stored in row-major order */
 
 #include "Woodbury.hpp"
 #include "Helpers.hpp"
@@ -23,30 +25,29 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   V[Updates_index[0] - 1] = 1;
   V[Dim + Updates_index[1] - 1] = 1;
 
-  showMatrix2(Slater_inv, Dim, Dim, "Slater_inv");
-  showMatrix2(Updates, 2, Dim, "Updates");
-  showMatrix2(Updates_index, 1, 2, "Updates_index");
-  showMatrix2(V, 2, Dim, "V");
-
-  // Compute C
+  // Compute C = S_inv U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE 
+  // OF LAYOUT OF 'Updates' !!
   double C[2 * Dim];
-  matMul2(Slater_inv, Updates, C, Dim, Dim, 2);
+  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++) {
+        C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[4 * j + k];
+      }
+    }
+  }
 
-  int A[6] = {1,2,3,4,5,6};
-  int Y[12] = {7,8,9,10,11,12,13,14,15,16,17,18};
-  int Z[8] = {0,0,0,0,0,0,0,0};
-  matMul2(A,Y,Z,2,3,4);
-  showMatrix2(Z, 2,4,"Z");
+  // Compute D = V x S^{-1}
+  double D[2 * Dim];
+  matMul2(V, Slater_inv, D, 2, Dim, Dim);
 
-  showMatrix2(C, 2, Dim, "C = S_inv * Updates");
-  // Compute B
+  // Compute B = 1 + V C
   double B[4];
   matMul2(V, C, B, 2, Dim, 2);
-  // Compute 1 + B
   B[0] += 1;
   B[3] += 1;
 
-  // Invert 1 + B with explicit formula for 2x2 inversion
+  // 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];
@@ -61,16 +62,19 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
     return false;
   }
 
-  // Compute (S + U * V)^{-1} with Woobury identity
-  double D[2 * Dim];
-  matMul2(V, Slater_inv, D, 2, Dim, Dim);
+  // Compute B^{-1} x D
   double tmp[2 * Dim];
   matMul2(Binv, D, tmp, 2, 2, Dim);
+
+  // Compute C x B^{-1} x D
   double tmp2[Dim * Dim];
   matMul2(C, tmp, tmp2, Dim, 2, Dim);
+
+  // 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];
   }
+
   return true;
 }
 

From 22590b7684a256fc5da9599e1d38a678f3fb7269 Mon Sep 17 00:00:00 2001
From: Francois Coppens <francois.coppens@uvsq.fr>
Date: Tue, 15 Jun 2021 11:53:04 +0200
Subject: [PATCH 6/6] * Woodburry 3x3 kernel fixed * Written Unified
 Sherman-Morrison-Woodbury kernel that partitions   the updates in blocks of 3
 and tries them with Woodbury 3x3.   The remainder of 2 or one are attempted
 with Woodbury 2x2 and SM2.

For now the unified kernel gives fails where pure SM2 does not.
I suspect there is something going wrong in how the updates are partitioned.
---
 include/Helpers.hpp            |  20 -----
 src/SMWB.cpp                   |  61 ++++++++++++--
 src/Woodbury.cpp               | 145 ++++++++++++++++++++++++---------
 tests/QMCChem_dataset_test.f90 |  28 +++++--
 tests/test_h5.cpp              |   2 -
 5 files changed, 181 insertions(+), 75 deletions(-)

diff --git a/include/Helpers.hpp b/include/Helpers.hpp
index cd6fdbb..f2a4c5d 100644
--- a/include/Helpers.hpp
+++ b/include/Helpers.hpp
@@ -75,26 +75,6 @@ void showMatrix2(T *matrix, unsigned int M, unsigned int N, std::string name) {
   std::cout << std::endl;
 }
 
-
-template <typename T>
-void showMatrixNS(T *matrix, unsigned int M, unsigned int N, std::string name) {
-  std::cout.precision(17);
-  std::cout << name << " = [" << std::endl;
-  for (unsigned int i = 0; i < M; i++) {
-    std::cout << "[";
-    for (unsigned int j = 0; j < N; j++) {
-      if (matrix[i * N + j] >= 0) {
-        std::cout << "  " << matrix[i * N + j] << ",";
-      } else {
-        std::cout << " " << matrix[i * N + j] << ",";
-      }
-    }
-    std::cout << " ]," << std::endl;
-  }
-  std::cout << "]" << std::endl;
-  std::cout << std::endl;
-}
-
 template <typename T> T *transpose(T *A, unsigned int M) {
   T *B = new T[M * M];
   for (unsigned int i = 0; i < M; i++) {
diff --git a/src/SMWB.cpp b/src/SMWB.cpp
index 670920b..a112dd7 100644
--- a/src/SMWB.cpp
+++ b/src/SMWB.cpp
@@ -8,13 +8,60 @@
 void SMWB1(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;
 
-           bool ok;
-           ok = WB2(Slater_inv, Dim, Updates, Updates_index);
-           if (!ok) {
-             std::cerr << "Woodbury kernel failed!" << std::endl;
-             SM2(Slater_inv, Dim, N_updates, Updates, Updates_index);
-           }
-         }
+  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;
+
+  // std::cerr << "Number of blocks: " << n_of_3blocks << ". Remainder: " << remainder << "." << std::endl;
+
+  // 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[length_3block];
+      unsigned int Updates_index_3block[3];
+      Updates_index_3block[0] = Updates_index[3 * i + 0];
+      Updates_index_3block[1] = Updates_index[3 * i + 1];
+      Updates_index_3block[2] = Updates_index[3 * i + 2];
+      for (unsigned int j = 0; j < length_3block; j++) {
+        Updates_3block[j] = Updates[i * length_3block + j];
+      }
+      bool ok;
+      ok = WB3(Slater_inv, Dim, Updates_3block, Updates_index_3block);
+      if (!ok) { // Send the entire block to SM2
+        std::cerr << "Woodbury 3x3 kernel failed! Sending block to SM2" << std::endl;
+        SM2(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[length_2block];
+    unsigned int Updates_index_2block[2];
+    Updates_index_2block[0] = Updates_index[3 * n_of_3blocks + 0];
+    Updates_index_2block[1] = Updates_index[3 * n_of_3blocks + 1];
+    for (unsigned int i = 0; i < length_2block; i++) {
+      Updates_2block[i] = Updates[n_of_3blocks * length_3block + i];
+    }
+    bool ok;
+    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;
+      SM2(Slater_inv, Dim, 2, Updates_2block, Updates_index_2block);
+    }
+  } else if (remainder == 1) { // Apply last remaining update with SM2
+    double Updates_1block[length_1block];
+    unsigned int Updates_index_1block[1];
+    Updates_index_1block[0] = Updates_index[3 * n_of_3blocks + 0];
+    for (unsigned int i = 0; i < length_1block; i++) {
+      Updates_1block[i] = Updates[n_of_3blocks * length_3block + i];
+    }
+    SM2(Slater_inv, Dim, 1, Updates_1block, Updates_index_1block);
+  } else { // remainder == 0
+    // Nothing left to do.
+  }
+}
 
 extern "C" {
 void SMWB1_f(double **linSlater_inv, unsigned int *Dim, unsigned int *N_updates,
diff --git a/src/Woodbury.cpp b/src/Woodbury.cpp
index 43eb301..dacf676 100644
--- a/src/Woodbury.cpp
+++ b/src/Woodbury.cpp
@@ -2,21 +2,22 @@
     Woodbury 2x2 and 3x3 kernels
 
     Woodbury matrix identity:
-
     (S + U V)^{-1}  =  S^{-1} - C B^{-1} D
 
-    C := S^{-1} U,    dim x 2
-    B := 1 + V C,     2 x 2
-    D := V S^{-1},    2 x dim
-
-    All matrices are stored in row-major order */
+    All matrices are stored in row-major order
+*/
 
 #include "Woodbury.hpp"
 #include "Helpers.hpp"
 
-// Woodbury 2x2 kernel:
+// 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
+  */
   std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
 
   // Construct V from Updates_index
@@ -25,23 +26,23 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   V[Updates_index[0] - 1] = 1;
   V[Dim + Updates_index[1] - 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++) {
       C[i * 2 + j] = 0;
       for(unsigned int k = 0; k < Dim; k++) {
-        C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[4 * j + k];
+        C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k];
       }
     }
   }
 
-  // Compute D = V x S^{-1}
+  // Compute D = V * S^{-1}
   double D[2 * Dim];
   matMul2(V, Slater_inv, D, 2, Dim, Dim);
 
-  // Compute B = 1 + V C
+  // Compute B = 1 + V * C
   double B[4];
   matMul2(V, C, B, 2, Dim, 2);
   B[0] += 1;
@@ -56,9 +57,9 @@ bool WB2(double *Slater_inv, unsigned int Dim, double *Updates,
   Binv[3] = idet * B[0];
 
   // Check if determinant of inverted matrix is not zero
-  double det = B[0] * B[3] - B[1] * B[2];
+  double det = Binv[0] * Binv[3] - Binv[1] * Binv[2];
   if (std::fabs(det) < threshold()) {
-    std::cerr << "Determinant approached 0!" << std::endl;
+    std::cerr << "Determinant approaching 0!" << std::endl;
     return false;
   }
 
@@ -79,60 +80,130 @@ 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
+  */
   std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
-  
+
+#ifdef DEBUG
+  showMatrix2(Slater_inv, Dim, Dim, "Slater_inv BEFORE update");
+  showMatrix2(Updates, 3, Dim, "Updates");
+  showMatrix2(Updates_index, 1, 3, "Updates_index");
+#endif
+
   // Construct V from Updates_index
-  unsigned int V[3 * Dim]; // 2 x Dim matrix stored in row-major order
+  unsigned int V[3 * Dim]; // 3 x Dim matrix stored in row-major order
+  std::memset(V, 0, 3 * Dim * sizeof(unsigned int));
   V[Updates_index[0] - 1] = 1;
   V[Dim + Updates_index[1] - 1] = 1;
   V[2 * Dim + Updates_index[2] - 1] = 1;
 
-  // Compute C
+#ifdef DEBUG
+  showMatrix2(V, 3, Dim, "V");
+#endif
+
+  // Compute C = S_inv * U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE 
+  // OF LAYOUT OF 'Updates' !!
   double C[3 * Dim];
-  matMul2(Slater_inv, Updates, C, Dim, Dim, 3);
-  // Compute B
+  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++) {
+        C[i * 3 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k];
+      }
+    }
+  }
+
+#ifdef DEBUG
+  showMatrix2(C, Dim, 3, "C = S_inv * U");
+#endif
+  // Compute D = V * S^{-1}
+  double D[3 * Dim];
+  matMul2(V, Slater_inv, D, 3, Dim, Dim);
+
+#ifdef DEBUG
+  showMatrix2(D, 3, Dim, "D = V * S_inv");
+#endif
+
+  // Compute B = 1 + V * C
   double B[9];
   matMul2(V, C, B, 3, Dim, 3);
-  // Compute 1 + B
   B[0] += 1;
   B[4] += 1;
   B[8] += 1;
 
-  double Binv[9];
-  Binv[0] = B[4] * B[8] - B[5] * B[7];
-  Binv[3] = B[5] * B[6] - B[3] * B[8];
-  Binv[6] = B[3] * B[7] - B[4] * B[6];
+#ifdef DEBUG
+  showMatrix2(B, 3, 3, "B = 1 + V * C");
+#endif
 
-  Binv[1] = B[2] * B[7] - B[1] * B[8];
-  Binv[4] = B[0] * B[8] - B[2] * B[6];
-  Binv[7] = B[1] * B[6] - B[0] * B[7];
+  // Compute B^{-1} with explicit formula for 3x3 inversion
+  double Binv[9], det, idet;
+  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;
 
-  Binv[2] = B[1] * B[5] - B[2] * B[4];
-  Binv[5] = B[2] * B[3] - B[0] * B[5];
-  Binv[8] = B[0] * B[4] - B[1] * B[3];
+#ifdef DEBUG 
+  std::cerr << "Determinant of B = " << det << std::endl;
+#endif
+
+  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 DEBUG
+  showMatrix2(Binv, 3, 3, "Binv");
+#endif
 
   // Check if determinant of inverted matrix is not zero
-  // If so, exigt and return false.
-  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]);
+  // 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 DEBUG
+  std::cerr << "Determinant of Binv = " << det << std::endl;
+#endif
+
   if (std::fabs(det) < threshold()) {
     std::cerr << "Determinant approached 0!" << std::endl;
     return false;
   }
 
-  // Compute (S + U * V)^{-1} with Woobury identity
-  double D[3 * Dim];
-  matMul2(V, Slater_inv, D, 3, Dim, Dim);
+  // Compute B^{-1} x D
   double tmp[3 * Dim];
   matMul2(Binv, D, tmp, 3, 3, Dim);
+
+#ifdef DEBUG
+  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 DEBUG
+  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];
   }
+
+#ifdef DEBUG
+  showMatrix2(Slater_inv, Dim, Dim, "Slater_inv AFTER update");
+#endif
   return true;
 }
 
diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90
index b8fac8f..3aae867 100644
--- a/tests/QMCChem_dataset_test.f90
+++ b/tests/QMCChem_dataset_test.f90
@@ -33,15 +33,15 @@ program QMCChem_dataset_test
     close(2000)
     close(3000)
 
-    !! Write Updates to file to check
-    open(unit = 2000, file = "Updates.dat")
-    do i=1,dim
-        do j=1,n_updates
-            write(2000,"(E23.15, 1X)", advance="no") Updates(i,j)
-        end do
-        write(2000,*)
-    end do
-    close(2000)
+    ! !! Write Updates to file to check
+    ! open(unit = 2000, file = "Updates.dat")
+    ! do i=1,dim
+    !     do j=1,n_updates
+    !         write(2000,"(E23.15, 1X)", advance="no") Updates(i,j)
+    !     end do
+    !     write(2000,*)
+    ! end do
+    ! close(2000)
 
     !! Update S & transform replacement updates 'Updates'
     !! into additive updates 'U' to compute the inverse
@@ -53,6 +53,16 @@ program QMCChem_dataset_test
         end do
     end do
 
+    !! Write Updates to file to check
+    open(unit = 2000, file = "Updates.dat")
+    do i=1,dim
+        do j=1,n_updates
+            write(2000,"(E23.15, 1X)", advance="no") U(i,j)
+        end do
+        write(2000,*)
+    end do
+    close(2000)
+
     !! Update S_inv
     !! S_inv needs to be transposed first before it
     !! goes to MaponiA3
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index 9b49a1c..f750f10 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -52,8 +52,6 @@ int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
   showMatrix(slater_inverse, dim, "OLD Inverse");
 #endif
 
-  showMatrix2(slater_matrix, dim, dim, "Slater");
-  
   // Transform replacement updates in 'updates[]' into additive updates in 'u[]'
   for (j = 0; j < nupdates; j++) {
     for (i = 0; i < dim; i++) {