From 7cef89e2b59cb370d63a2e53cb980abcddeff32e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <mail@fmgjcoppens.nl>
Date: Fri, 16 Apr 2021 17:19:55 +0200
Subject: [PATCH 1/5] Slagel-splitting implemented but not tested yet.

---
 src/SM_MaponiA3S.cpp | 46 +++++++++++++++++++++++---------------------
 1 file changed, 24 insertions(+), 22 deletions(-)

diff --git a/src/SM_MaponiA3S.cpp b/src/SM_MaponiA3S.cpp
index 66efd37..bb02872 100644
--- a/src/SM_MaponiA3S.cpp
+++ b/src/SM_MaponiA3S.cpp
@@ -19,6 +19,10 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   double *next = new double[Dim * Dim]{0};
   double *last = Slater_inv, *tmp;
 
+  double later_updates[Dim * N_updates];
+  unsigned int later_index[N_updates];
+  unsigned int later = 0;
+
   // Populate update-order vector
   for (i = 0; i < N_updates; i++) {
     p[i + 1] = i + 1;
@@ -79,18 +83,15 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
       std::cerr << "Breakdown condition triggered at " << component
                 << std::endl;
       
-    //   for (unsigned int i = 0; i < Dim; i++) {
-    //     later_updates[later * Dim + i] = Updates[l * Dim + i] / 2.0;
-    //     ylk[l][p[l + 1]][i] /= 2.0;
-    //   }
-    //   later_index[later] = Updates_index[l];
-    //   later++;
-
-    //   den = 1 + C[Updates_index[l] - 1];
-    // }
-    // double iden = 1 / den;
-
+      for (unsigned int i = 0; i < Dim; i++) {
+        later_updates[later * Dim + i] = Updates[l * Dim + i] * 0.5;
+        ylk[l][p[l + 1]][i] *= 0.5;
+      }
+      later_index[later] = Updates_index[p[l + 1]];
+      later++;
+      beta = 1 + ylk[l][p[l + 1]][component];
     }
+    double ibeta = 1.0 / beta;
 
 // Compute intermediate update to Slater_inv
 #ifdef DEBUG
@@ -105,20 +106,17 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     for (i = 0; i < Dim; i++) {
       for (j = 0; j < Dim; j++) {
         Al[i * Dim + j] =
-            (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
+            (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta;
       }
     }
-    matMul(Al, last, next, Dim);
-    tmp = next;
-    next = last;
-    last = tmp;
+    matMul(Al, last, next, Dim); tmp = next; next = last; last = tmp;
 #ifdef DEBUG
     showMatrix(last, Dim, "last");
 #endif
 
     // For given l != 0 compute the next {y_{l,k}}
     for (k = l + 2; k < N_updates + 1; k++) {
-      alpha = ylk[l][p[k]][component] / beta;
+      alpha = ylk[l][p[k]][component] * ibeta;
 #ifdef DEBUG
       std::cout << "Inside k-loop: k = " << k << std::endl;
       std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl;
@@ -131,6 +129,10 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   }
   memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
 
+  if (later > 0) {
+    MaponiA3S(Slater_inv, Dim, later, later_updates, later_index);
+  }
+
   /*
   CLEANUP MEMORY
   */
@@ -145,9 +147,9 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 }
 
 extern "C" {
-void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
-                 unsigned int *N_updates, double **linUpdates,
-                 unsigned int **Updates_index) {
-  MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
-}
+  void MaponiA3S_f(double **linSlater_inv, unsigned int *Dim,
+                  unsigned int *N_updates, double **linUpdates,
+                  unsigned int **Updates_index) {
+    MaponiA3S(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
 }

From af64b91ca6d78d50b4c23994e95dd6954e2ff8b8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <mail@fmgjcoppens.nl>
Date: Fri, 16 Apr 2021 17:33:15 +0200
Subject: [PATCH 2/5] Started debugging Slagel-splitting in MaponiA3.

---
 src/SM_MaponiA3S.cpp | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/SM_MaponiA3S.cpp b/src/SM_MaponiA3S.cpp
index bb02872..634ba3a 100644
--- a/src/SM_MaponiA3S.cpp
+++ b/src/SM_MaponiA3S.cpp
@@ -129,10 +129,6 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   }
   memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
 
-  if (later > 0) {
-    MaponiA3S(Slater_inv, Dim, later, later_updates, later_index);
-  }
-
   /*
   CLEANUP MEMORY
   */
@@ -144,6 +140,11 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     delete[] ylk[l];
   }
   delete[] Al, next, p;
+
+  if (later > 0) {
+    std::cout << "Entering recursive loop with " << l << " updates" << std::endl;
+    MaponiA3S(Slater_inv, Dim, later, later_updates, later_index);
+  }
 }
 
 extern "C" {

From fb3ed9ce0c011b548a0d6a659d39460f1513c3ad Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <mail@fmgjcoppens.nl>
Date: Mon, 19 Apr 2021 08:28:19 +0200
Subject: [PATCH 3/5] In Maponi A3, converted divisions into multiplications.

---
 src/SM_MaponiA3.cpp | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp
index 8d62e42..71e1bf2 100644
--- a/src/SM_MaponiA3.cpp
+++ b/src/SM_MaponiA3.cpp
@@ -79,6 +79,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     if (fabs(beta) < threshold()) {
       std::cerr << "Break-down occured." << std::endl;
     }
+    double ibeta = 1.0 / beta; 
 
 // Compute intermediate update to Slater_inv
 #ifdef DEBUG
@@ -93,7 +94,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     for (i = 0; i < Dim; i++) {
       for (j = 0; j < Dim; j++) {
         Al[i * Dim + j] =
-            (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] / beta;
+            (i == j) - (j == component - 1) * ylk[l][p[l + 1]][i + 1] * ibeta;
       }
     }
     matMul(Al, last, next, Dim);
@@ -106,7 +107,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 
     // For given l != 0 compute the next {y_{l,k}}
     for (k = l + 2; k < N_updates + 1; k++) {
-      alpha = ylk[l][p[k]][component] / beta;
+      alpha = ylk[l][p[k]][component] * ibeta;
 #ifdef DEBUG
       std::cout << "Inside k-loop: k = " << k << std::endl;
       std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl;

From e3815267fbeab1fd6073890bc54de1cc108b7295 Mon Sep 17 00:00:00 2001
From: Francois Coppens <francois.coppens@uvsq.fr>
Date: Tue, 20 Apr 2021 16:08:09 +0200
Subject: [PATCH 4/5] * Added SM{1,2,3} algorithms to Fortran module * Deleted
 final ylk array in MaponiA3(S) * Output of Valgrind seem to suggest there are
 no major memory problems.

---
 Makefile                       |  2 +-
 Makefile.verificarlo           |  2 +-
 src/SM_MaponiA3.cpp            | 12 +++++-----
 src/SM_MaponiA3S.cpp           |  4 ++--
 src/SM_MaponiA3_mod.f90        | 12 ----------
 src/SM_Standard.cpp            | 20 +++++++++++++++++
 src/SM_mod.f90                 | 40 ++++++++++++++++++++++++++++++++++
 tests/QMCChem_dataset_test.f90 |  4 ++--
 tests/test_h5.cpp              |  2 +-
 todo.txt                       | 10 +++++++++
 10 files changed, 83 insertions(+), 25 deletions(-)
 delete mode 100644 src/SM_MaponiA3_mod.f90
 create mode 100644 src/SM_mod.f90
 create mode 100644 todo.txt

diff --git a/Makefile b/Makefile
index fc39771..439b5e1 100644
--- a/Makefile
+++ b/Makefile
@@ -34,7 +34,7 @@ DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
 		   $(OBJ_DIR)/SM_Standard.o \
 		   $(OBJ_DIR)/SM_Helpers.o
 DEPS_F = $(DEPS_CXX) \
-		 $(OBJ_DIR)/SM_MaponiA3_mod.o \
+		 $(OBJ_DIR)/SM_mod.o \
 		 $(OBJ_DIR)/Helpers_mod.o
 
 ## Directory structure
diff --git a/Makefile.verificarlo b/Makefile.verificarlo
index bd10702..7e5c347 100644
--- a/Makefile.verificarlo
+++ b/Makefile.verificarlo
@@ -21,7 +21,7 @@ DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o \
 		   $(OBJ_DIR)/SM_Standard.o \
 		   $(OBJ_DIR)/SM_Helpers.o
 DEPS_F = $(DEPS_CXX) \
-		 $(OBJ_DIR)/SM_MaponiA3_mod.o \
+		 $(OBJ_DIR)/SM_mod.o \
 		 $(OBJ_DIR)/Helpers_mod.o
 FLIBS = -lstdc++
 
diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp
index 71e1bf2..f16f625 100644
--- a/src/SM_MaponiA3.cpp
+++ b/src/SM_MaponiA3.cpp
@@ -130,13 +130,13 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     }
     delete[] ylk[l];
   }
-  delete[] Al, next, p;
+  delete[] Al, next, p, ylk;
 }
 
 extern "C" {
-void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
-                unsigned int *N_updates, double **linUpdates,
-                unsigned int **Updates_index) {
-  MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
-}
+  void MaponiA3_f(double **linSlater_inv, unsigned int *Dim,
+                  unsigned int *N_updates, double **linUpdates,
+                  unsigned int **Updates_index) {
+    MaponiA3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
 }
diff --git a/src/SM_MaponiA3S.cpp b/src/SM_MaponiA3S.cpp
index 634ba3a..8c58ef4 100644
--- a/src/SM_MaponiA3S.cpp
+++ b/src/SM_MaponiA3S.cpp
@@ -4,7 +4,7 @@
 #include "SM_MaponiA3S.hpp"
 #include "SM_Helpers.hpp"
 
-// #define DEBUG
+#define DEBUG
 
 void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
                double *Updates, unsigned int *Updates_index) {
@@ -139,7 +139,7 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     }
     delete[] ylk[l];
   }
-  delete[] Al, next, p;
+  delete[] Al, next, p, ylk;
 
   if (later > 0) {
     std::cout << "Entering recursive loop with " << l << " updates" << std::endl;
diff --git a/src/SM_MaponiA3_mod.f90 b/src/SM_MaponiA3_mod.f90
deleted file mode 100644
index a7105b2..0000000
--- a/src/SM_MaponiA3_mod.f90
+++ /dev/null
@@ -1,12 +0,0 @@
-module Sherman_Morrison
-    use, intrinsic :: iso_c_binding, only : c_int, c_double
-    interface
-        subroutine MaponiA3(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_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 MaponiA3
-    end interface
-end module Sherman_Morrison
diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp
index 2591b07..d4ba046 100644
--- a/src/SM_Standard.cpp
+++ b/src/SM_Standard.cpp
@@ -170,3 +170,23 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     SM3(Slater_inv, Dim, later, later_updates, later_index);
   }
 }
+
+extern "C" {
+  void SM1_f(double **linSlater_inv, unsigned int *Dim,
+                  unsigned int *N_updates, double **linUpdates,
+                  unsigned int **Updates_index) {
+    SM1(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
+
+  void SM2_f(double **linSlater_inv, unsigned int *Dim,
+                  unsigned int *N_updates, double **linUpdates,
+                  unsigned int **Updates_index) {
+    SM2(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
+  
+  void SM3_f(double **linSlater_inv, unsigned int *Dim,
+                  unsigned int *N_updates, double **linUpdates,
+                  unsigned int **Updates_index) {
+    SM3(*linSlater_inv, *Dim, *N_updates, *linUpdates, *Updates_index);
+  }
+}
diff --git a/src/SM_mod.f90 b/src/SM_mod.f90
new file mode 100644
index 0000000..f947458
--- /dev/null
+++ b/src/SM_mod.f90
@@ -0,0 +1,40 @@
+module Sherman_Morrison
+    use, intrinsic :: iso_c_binding, only : c_int, c_double
+    interface
+        subroutine MaponiA3(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3_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 MaponiA3
+        subroutine MaponiA3S(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="MaponiA3S_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 MaponiA3S
+        subroutine SM1(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM1_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 SM1
+        subroutine SM2(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM2_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 SM2
+        subroutine SM3(Slater_inv, dim, n_updates, Updates, Updates_index) bind(C, name="SM3_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 SM3
+    end interface
+end module Sherman_Morrison
diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90
index 303b578..b8fac8f 100644
--- a/tests/QMCChem_dataset_test.f90
+++ b/tests/QMCChem_dataset_test.f90
@@ -57,12 +57,12 @@ program QMCChem_dataset_test
     !! S_inv needs to be transposed first before it
     !! goes to MaponiA3
     call Transpose(S_inv, S_inv_t, dim)
-    call MaponiA3(S_inv_t, dim, n_updates, U, Updates_index)
+    ! call MaponiA3S(S_inv_t, dim, n_updates, U, Updates_index)
+    call SM2(S_inv_t, dim, n_updates, U, Updates_index)
     !! S_inv_t needs to be transposed back before it
     !! can be multiplied with S to test unity
     call Transpose(S_inv_t, S_inv, dim)
 
-
     !! Write new S and S_inv to file for check in Octave
     open(unit = 4000, file = "Slater.dat")
     open(unit = 5000, file = "Slater_inv.dat")
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index f965b91..9e3f49a 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -61,7 +61,6 @@ int test_cycle(H5File file, int cycle, std::string version) {
 
   double * slater_inverse = new double[dim*dim];
   read_double(file, group + "/slater_inverse", slater_inverse);
-  //slater_inverse = transpose(slater_inverse, dim);
 
   unsigned int * col_update_index = new unsigned int[nupdates];
   read_int(file, group + "/col_update_index", col_update_index);
@@ -80,6 +79,7 @@ int test_cycle(H5File file, int cycle, std::string version) {
   showMatrix(slater_inverse, dim, "OLD Inverse");
 #endif
 
+  // Transform replacement updates in 'updates[]' into additive updates in 'u[]'
   for (j = 0; j < nupdates; j++) {
     for (i = 0; i < dim; i++) {
       col = col_update_index[j];
diff --git a/todo.txt b/todo.txt
new file mode 100644
index 0000000..5a579c0
--- /dev/null
+++ b/todo.txt
@@ -0,0 +1,10 @@
+-=={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 3ee4d3b56daf1e03c5d2aed1e27d013945755406 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= <mail@fmgjcoppens.nl>
Date: Wed, 5 May 2021 11:51:22 +0200
Subject: [PATCH 5/5] - Added a small 3x3 example that will always become
 singular regardles the order of the applied updates. This example is added to
 test Maponi A3 augmented with Slagel-splitting. - Slagel-splitting in Maponi
 A3 is now working. A thourough analysis on the QMC=Chem datasets still has te
 be done. - Acceptance tollerance  in 'test_h5' on the residuals can be set at
 runtime now.

---
 datasets/small/dataset.dat  |  32 +++++++++++++++--
 datasets/small/dataset.hdf5 | Bin 30583 -> 41231 bytes
 src/SM_MaponiA3S.cpp        |  66 +++++++++++++++++++++---------------
 src/SM_Standard.cpp         |   3 +-
 tests/run_test.sh           |   7 ++--
 tests/test_h5.cpp           |  27 +++++++++------
 6 files changed, 92 insertions(+), 43 deletions(-)

diff --git a/datasets/small/dataset.dat b/datasets/small/dataset.dat
index a393c0b..1d0d8f0 100644
--- a/datasets/small/dataset.dat
+++ b/datasets/small/dataset.dat
@@ -1,5 +1,5 @@
 #START_PACKET
-#CYCLE_ID:   13
+#CYCLE_ID:   1
 #SLATER_MATRIX_DIM:  21
 #NUPDATES:  3
 #SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
@@ -513,7 +513,7 @@
 #END_PACKET
 
 #START_PACKET
-#CYCLE_ID:   8169
+#CYCLE_ID:   2
 #SLATER_MATRIX_DIM:  4
 #NUPDATES:  2
 #SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
@@ -544,3 +544,31 @@
 #COL_UPDATE_COMP_(03):   0.000000000000000E+00
 #COL_UPDATE_COMP_(04):  -1.000000000000000E+00
 #END_PACKET
+
+#START_PACKET
+#CYCLE_ID:   3
+#SLATER_MATRIX_DIM:  3
+#NUPDATES:  3
+#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet
+(01,01)    1.000000000000000E+00    1.000000000000000E+00
+(01,02)    0.000000000000000E+00    0.000000000000000E+00
+(01,03)    0.000000000000000E+00    0.000000000000000E+00
+(02,01)    0.000000000000000E+00    0.000000000000000E+00
+(02,02)    1.000000000000000E+00    1.000000000000000E+00
+(02,03)    0.000000000000000E+00    0.000000000000000E+00
+(03,01)    0.000000000000000E+00    0.000000000000000E+00
+(03,02)    0.000000000000000E+00    0.000000000000000E+00
+(03,03)   -1.000000000000000E+00   -1.000000000000000E+00
+#COL_UPDATE_INDEX: 1
+#COL_UPDATE_COMP_(01):   1.000000000000000E+00
+#COL_UPDATE_COMP_(02):   1.000000000000000E+00
+#COL_UPDATE_COMP_(03):  -1.000000000000000E+00
+#COL_UPDATE_INDEX: 2
+#COL_UPDATE_COMP_(01):   1.000000000000000E+00
+#COL_UPDATE_COMP_(02):   1.000000000000000E+00
+#COL_UPDATE_COMP_(03):   0.000000000000000E+00
+#COL_UPDATE_INDEX: 3
+#COL_UPDATE_COMP_(01):   1.000000000000000E+00
+#COL_UPDATE_COMP_(02):   0.000000000000000E+00
+#COL_UPDATE_COMP_(03):  -1.000000000000000E+00
+#END_PACKET
diff --git a/datasets/small/dataset.hdf5 b/datasets/small/dataset.hdf5
index 872ecb8a576bc63a3db6f31d9de65452c7c11057..a583fe978f8219eb33a9827505ca7e4163b8b1e3 100644
GIT binary patch
delta 735
zcmezVj<Np`(*zAB{)HR0`j{9MChugDlVV7&OwLJ-H-gf}42%pAU@>u{{$>T{2}~0c
zB;?D>A%fE@p|Vga*w5dEk%4*gb5>18fyshw>Piw&nGLKE2Fy5w2@^M-o7^B6G`Wgh
ziTU=U`pxs%ty$ujq1s?{MO<<MLnHsU1TJ+2js5Zr6-Qb4u!(_;X98l7K7_pEd-uoA
znGy#cyk}<!H{kTz{E*v&6=-P*$im6<*=J6kAan%Gj&n}jyhB(*6zr0B@q&z-pJ=K}
zfO&jQzKmc-oO3TPSn+)Jm7Avrt^tceO;MQqKX?gP5@Zs4=o2ommUr=6xxu1*PAez-
z=O9#s=PGbbULede`LVAIlYI4NMn5qort_T>H^za&gy|OdWXDjE$rIRlU=}iMpE7Zy
z9G3#ne;|-vGx=bk8q@v`ke~rXjG?@CvS5%J6Mq#*Fas)heCp)HpaV>g&`kvSU19Qa
zM<I5ILgvXAoU(Y^>!3U!#XQ+TSOn%LMv;jd<(Xb_Pk!hu#>l++qwsT1tOnF~VKHDc
zC`2J9{6jIJK~QM&OveBwA)d*v9TlJ<wQ=%#M=lgsU7v=<(8<7%W|RWCN{$D~RRCPh
B;2{72

delta 245
zcmeA_#Pt0g;{*++^74&ZeN2oRlXo)7DH$_BL2_ksPHMb`p_wH_fDtNgFma>$W(DR6
zOcS>oOqQ=^V`Q5AoK=%iVzMBc`o!%YCN~HMO|D{BV&1`9v3WkbHOpp3k!7qv(Grlz
z<oWC~Cr=PM0%pfKCvM&$EFlV3`Yv9Oaq|;RbqO$!&&ih&%!qUD<pnFA&%Sc=^uRS>
pQK%^jlm7=V0ZW2RVh?@71=jK|ek(Uvl+S79Wd9t5itt<oE&!0WQGx&f

diff --git a/src/SM_MaponiA3S.cpp b/src/SM_MaponiA3S.cpp
index 8c58ef4..0b8c4f3 100644
--- a/src/SM_MaponiA3S.cpp
+++ b/src/SM_MaponiA3S.cpp
@@ -4,14 +4,14 @@
 #include "SM_MaponiA3S.hpp"
 #include "SM_Helpers.hpp"
 
-#define DEBUG
+// #define DEBUG
 
 void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
                double *Updates, unsigned int *Updates_index) {
   /*
   DECLARE AND INITIALISE ARRAYS
   */
-
+  showVector(Updates_index, Dim, "Updates_index");
   unsigned int k, l, i, j, component;
   unsigned int *p = new unsigned int[N_updates + 1]{0};
   double alpha, beta;
@@ -40,13 +40,14 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   /*
   START ALGORITHM
   */
+ showMatrix(Slater_inv, Dim, "S^{-1}");
 
-  // Calculate the {y_{0,k}}
+  // Calculate the y_{0,k}  = S_0^{-1} u_k
   for (k = 1; k < N_updates + 1; k++) {
 #ifdef DEBUG
-    std::cout << "Compute y0k: " << std::endl;
-    std::cout << "ylk[0][" << k << "][:]";
-    std::cout << std::endl;
+    std::cerr << "Compute y0k: " << std::endl;
+    std::cerr << "ylk[0][" << k << "][:]";
+    std::cerr << std::endl;
 #endif
     for (i = 1; i < Dim + 1; i++) {
       for (j = 1; j < Dim + 1; j++) {
@@ -55,15 +56,15 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
       }
     }
 #ifdef DEBUG
-    showVector(ylk[0][k], Dim, "");
+    showVector(ylk[0][k], Dim+1, "");
 #endif
   }
 
   // Calculate the {y_{l,k}} from the {y_{0,k}}
   for (l = 0; l < N_updates; l++) {
 #ifdef DEBUG
-    std::cout << "In outer compute-ylk-loop: l = " << l << std::endl;
-    std::cout << std::endl;
+    std::cerr << "In outer compute-ylk-loop: l = " << l << std::endl;
+    std::cerr << std::endl;
 #endif
 
     // For given l select intermediate update with largest break-down val
@@ -73,35 +74,44 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     component = Updates_index[p[l + 1] - 1];
     beta = 1 + ylk[l][p[l + 1]][component];
 #ifdef DEBUG
-    std::cout << "p[l+1] = " << p[l + 1] << std::endl;
-    std::cout << "component = " << component << std::endl;
-    std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
+    std::cerr << "p[l+1] = " << p[l + 1] << std::endl;
+    std::cerr << "component = " << component << std::endl;
+    std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
               << "] = " << beta << std::endl;
-    std::cout << std::endl;
+    std::cerr << std::endl;
 #endif
     if (fabs(beta) < threshold()) {
-      std::cerr << "Breakdown condition triggered at " << component
+#ifdef DEBUG
+      std::cerr << "Breakdown condition triggered at l = " << l << " and component "<< component
                 << std::endl;
-      
-      for (unsigned int i = 0; i < Dim; i++) {
-        later_updates[later * Dim + i] = Updates[l * Dim + i] * 0.5;
+#endif
+      for (unsigned int i = 1; i < Dim + 1; i++) {
+        std::cout << "Let's split..." << std::endl;
         ylk[l][p[l + 1]][i] *= 0.5;
+        // later_updates[later * Dim + i - 1] =  ylk[l][p[l + 1]][i]; // Updates[l * Dim + i] * 0.5;
+        later_updates[later * Dim + i - 1] = Updates[l * Dim + i - 1] * 0.5;
       }
-      later_index[later] = Updates_index[p[l + 1]];
+      later_index[later] = Updates_index[p[l + 1] - 1];
+      std::cout << "Later = " << later << std::endl;
       later++;
+      std::cout << "Later = " << later << std::endl;
       beta = 1 + ylk[l][p[l + 1]][component];
+      std::cout << "New beta is: " << beta << std::endl;
     }
     double ibeta = 1.0 / beta;
+    showVector(ylk[l][p[l + 1]], Dim+1, "");
+    showVector(later_updates+(Dim*0), Dim, "later_updateds");
+    
 
 // Compute intermediate update to Slater_inv
 #ifdef DEBUG
-    std::cout << "Compute intermediate update to Slater_inv" << std::endl;
-    std::cout << "component = " << component << std::endl;
-    std::cout << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
+    std::cerr << "Compute intermediate update to Slater_inv" << std::endl;
+    std::cerr << "component = " << component << std::endl;
+    std::cerr << "beta = 1 + ylk[" << l << "][" << p[l + 1] << "][" << component
               << "]" << std::endl;
-    std::cout << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
+    std::cerr << "ylk[l][p[k]][:] = ylk[" << l << "][" << p[l + 1] << "][:]"
               << std::endl;
-    std::cout << std::endl;
+    std::cerr << std::endl;
 #endif
     for (i = 0; i < Dim; i++) {
       for (j = 0; j < Dim; j++) {
@@ -118,13 +128,14 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     for (k = l + 2; k < N_updates + 1; k++) {
       alpha = ylk[l][p[k]][component] * ibeta;
 #ifdef DEBUG
-      std::cout << "Inside k-loop: k = " << k << std::endl;
-      std::cout << "ylk[" << l + 1 << "][" << p[k] << "][:]" << std::endl;
-      std::cout << std::endl;
+      std::cerr << "Inside k-loop: k = " << k << std::endl;
+      std::cerr << "ylk[" << l + 1 << "][" << p[k] << "][:]";
+      std::cerr << std::endl;
 #endif
       for (i = 1; i < Dim + 1; i++) {
         ylk[l + 1][p[k]][i] = ylk[l][p[k]][i] - alpha * ylk[l][p[l + 1]][i];
       }
+      showVector(ylk[l + 1][p[k]], Dim+1, "");
     }
   }
   memcpy(Slater_inv, last, Dim * Dim * sizeof(double));
@@ -142,7 +153,8 @@ void MaponiA3S(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
   delete[] Al, next, p, ylk;
 
   if (later > 0) {
-    std::cout << "Entering recursive loop with " << l << " updates" << std::endl;
+    std::cerr << "Entering recursive loop with " << later << " updates" << std::endl;
+    std::cout << later_index[0] << std::endl;
     MaponiA3S(Slater_inv, Dim, later, later_updates, later_index);
   }
 }
diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp
index d4ba046..b188a2a 100644
--- a/src/SM_Standard.cpp
+++ b/src/SM_Standard.cpp
@@ -162,7 +162,8 @@ 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) {
-    std::cerr << "SM3 cannot invert this matrix" << std::endl;
+    std::cerr << "SM3 cannot invert matrix." << std::endl;
+    showMatrix(Slater_inv, Dim, "Slater_inverse");
     return;
   }
   // If some have failed, make a recursive call
diff --git a/tests/run_test.sh b/tests/run_test.sh
index 01c425b..c23676d 100755
--- a/tests/run_test.sh
+++ b/tests/run_test.sh
@@ -12,6 +12,7 @@ ALGO=$1
 START=$2
 STOP=$3
 BLOCKSIZE=329
+TOLERANCE=1e-3
 
 BLOCKS=$((STOP / BLOCKSIZE))
 LEFT=$((STOP % BLOCKSIZE))
@@ -25,17 +26,17 @@ then
     do
         BSTART=$(((i-1)*BLOCKSIZE+1))
         BSTOP=$((i*BLOCKSIZE-1))
-        $TEST $ALGO $BSTART $BSTOP
+        $TEST $ALGO $BSTART $BSTOP $TOLERANCE
         LAST=$i
     done
     LSTART=$((LAST*BLOCKSIZE+1))
     LSTOP=$((LSTART+LEFT-1))
-    $TEST $ALGO $LSTART $LSTOP
+    $TEST $ALGO $LSTART $LSTOP $TOLERANCE
 else
     for ((i=1; i<=$BLOCKS; i++))
     do
         BSTART=$(((i-1)*BLOCKSIZE+1))
         BSTOP=$((i*BLOCKSIZE-1))
-        $TEST $ALGO $BSTART $BSTOP
+        $TEST $ALGO $BSTART $BSTOP $TOLERANCE
     done
 fi
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index 9e3f49a..8f5c5d2 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -7,7 +7,7 @@
 #include "SM_Helpers.hpp"
 
 using namespace H5;
-// #define DEBUG
+#define DEBUG
 
 const H5std_string FILE_NAME( "dataset.hdf5" );
 
@@ -46,7 +46,7 @@ void read_double(H5File file, std::string key, double * data) {
   ds.close();
 }
 
-int test_cycle(H5File file, int cycle, std::string version) {
+int test_cycle(H5File file, int cycle, std::string version, double tolerance) {
 
   /* Read the data */
 
@@ -71,10 +71,6 @@ int test_cycle(H5File file, int cycle, std::string version) {
   double * u = new double[nupdates*dim];
 
   /* Test */
-#ifdef DEBUG
-  showMatrix(slater_matrix, dim, "OLD Slater");
-#endif
-
 #ifdef DEBUG
   showMatrix(slater_inverse, dim, "OLD Inverse");
 #endif
@@ -88,6 +84,16 @@ int test_cycle(H5File file, int cycle, std::string version) {
     }
   }
 
+#ifdef DEBUG
+  showMatrix(slater_matrix, dim, "OLD Slater");
+#endif
+
+#ifdef DEBUG
+  showMatrix(u, dim, "Updates");
+#endif
+
+
+
   if (version == "maponia3") {
     MaponiA3(slater_inverse, dim, nupdates, u, col_update_index);
   } else if (version == "maponia3s") {
@@ -113,7 +119,7 @@ int test_cycle(H5File file, int cycle, std::string version) {
 
   double * res = new double[dim*dim] {0};
   matMul(slater_matrix, slater_inverse, res, dim);
-  bool ok = is_identity(res, dim, 1e-3);
+  bool ok = is_identity(res, dim, tolerance);
 
   double res_max = residual_max(res, dim);
   double res2 = residual2(res, dim);
@@ -130,19 +136,20 @@ int test_cycle(H5File file, int cycle, std::string version) {
 }
 
 int main(int argc, char **argv) {
-  if (argc != 4) {
+  if (argc != 5) {
     std::cerr << "Execute from within 'datasets/'" << std::endl;
-    std::cerr << "usage: test_h5 <version> <start cycle> <stop cycle>" << std::endl;
+    std::cerr << "usage: test_h5 <version> <start cycle> <stop cycle> <tolerance>" << std::endl;
     return 1;
   }
   std::string version(argv[1]);
   int start_cycle = std::stoi(argv[2]);
   int stop_cycle = std::stoi(argv[3]);
+  double tolerance = std::stod(argv[4]);
   H5File file(FILE_NAME, H5F_ACC_RDONLY);
 
   bool ok;
   for (int cycle = start_cycle; cycle < stop_cycle+1; cycle++) {
-    ok = test_cycle(file, cycle, version);
+    ok = test_cycle(file, cycle, version, tolerance);
     if (ok) {
       std::cerr << "ok -- cycle " << std::to_string(cycle)
       << std::endl;