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 +