* 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

@@ -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

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++
 

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,
+  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);
-}
+  }
 }

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;

src/SM_MaponiA3_mod.f90

@@ -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

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);
+  }
+}

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

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")

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];

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
+