From ef97b40196e86007ae9562994371174dbbdef355 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Fri, 12 Mar 2021 12:38:50 +0100 Subject: [PATCH 1/3] Changed code to treat the dataset as replacement updates, instead of additive updates. Before the updates are used to compute the inverse they are transformed accordingly. However, this breaks the 4x4 example in update cycle 8169. --- Makefile | 16 ++++++++-------- octave-tests/QMCChem_dataset_test.m | 20 +++++++++++++++++++- smvarsrc | 3 +++ src/det.irp.f | 5 +---- tests/QMCChem_dataset_test.f90 | 17 ++++++++++------- 5 files changed, 41 insertions(+), 20 deletions(-) diff --git a/Makefile b/Makefile index c4fc58b..1249960 100644 --- a/Makefile +++ b/Makefile @@ -1,18 +1,18 @@ ## Compilers -ARCH = -H5CXX = h5c++ -std=gnu++11 -CXX = g++ -FC = gfortran +ARCH = +CXX = icpc +FC = ifort +H5CXX = h5c++ ## Compiler flags -H5CXXFLAGS = -O0 -g -CXXFLAGS = -O0 -g -FFLAGS = -O0 -g +CXXFLAGS = -O0 -debug full -traceback +FFLAGS = $(CXXFLAGS) +H5CXXFLAGS = $(CXXFLAGS) -fPIC +FLIBS = -lstdc++ INCLUDE = -I $(INC_DIR)/ DEPS_CXX = $(OBJ_DIR)/SM_MaponiA3.o $(OBJ_DIR)/SM_Standard.o DEPS_F = $(DEPS_CXX) $(OBJ_DIR)/SM_MaponiA3_mod.o $(OBJ_DIR)/Helpers_mod.o -FLIBS = -lstdc++ SRC_DIR := src TST_DIR := tests diff --git a/octave-tests/QMCChem_dataset_test.m b/octave-tests/QMCChem_dataset_test.m index 2861551..b43ab74 100644 --- a/octave-tests/QMCChem_dataset_test.m +++ b/octave-tests/QMCChem_dataset_test.m @@ -17,4 +17,22 @@ printf("--------------------------------------\n") printf("Determinant of S x S_inv : %f\n", det(S*S_inv)) printf("Trace of S x S_inv : %f\n", trace(S*S_inv)) printf("Norm of S x S_inv : %f\n", norm(S*S_inv)) -printf("======================================\n") + +cutoff = 1e-6; +printf("\n") +printf("Cutoff set to %e: S x S_inv = \n", cutoff) +printf(" -----------------------------------------\n") +dim = columns(S); +res=S*S_inv; +for i = 1:dim + for j = 1:dim + if (res(i,j) < cutoff) + res(i,j) = 0; + elseif (S(i,j)-1 < cutoff) + res(i,j) = 1; + endif + endfor +endfor +format free; +disp(res); +printf(" =========================================\n") diff --git a/smvarsrc b/smvarsrc index 0265f11..a4434e3 100644 --- a/smvarsrc +++ b/smvarsrc @@ -1,3 +1,6 @@ #!/usr/bin/env bash +export HDF5_CXX=icpc +export HDF5_CXXLINKER=icpc + export PATH=$PWD/bin:$PATH diff --git a/src/det.irp.f b/src/det.irp.f index 10a0968..e83d64a 100644 --- a/src/det.irp.f +++ b/src/det.irp.f @@ -970,7 +970,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d) !DIR$ ASSUME (MOD(LDS,$IRP_ALIGN/8) == 0) !DIR$ ASSUME (n>150) - integer :: i,j,n4 + integer :: i,j,n4 double precision :: zl !DIR$ NOPREFETCH @@ -1309,9 +1309,6 @@ END_PROVIDER write(501,10007) write(502,10007) cycle_id = cycle_id + 1 - !close(501) !! Close file - !close(502) !! Close file - !stop else diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90 index 508c0f3..2f537a8 100644 --- a/tests/QMCChem_dataset_test.f90 +++ b/tests/QMCChem_dataset_test.f90 @@ -6,7 +6,7 @@ program QMCChem_dataset_test integer :: i, j, col !! Iterators integer :: cycle_id, dim, n_updates integer(c_int), dimension(:), allocatable :: Updates_index - real(c_double), dimension(:,:), allocatable :: Updates + real(c_double), dimension(:,:), allocatable :: Updates, U real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_t call Read_dataset("datasets/update_cycle_8169.dat", & @@ -17,7 +17,7 @@ program QMCChem_dataset_test S_inv, & Updates_index, & Updates) - allocate(S_inv_t(dim,dim)) + allocate(S_inv_t(dim,dim), U(dim,n_updates)) !! Write current S and S_inv to file for check in Octave open(unit = 2000, file = "Slater_old.dat") @@ -43,11 +43,13 @@ program QMCChem_dataset_test end do close(2000) - !! Update S + !! Update S & transform replacement updates 'Updates' + !! into additive updates 'U' to compute the inverse do j=1,n_updates do i=1,dim col = Updates_index(j) - S(i,col) = S(i,col) + Updates(i,j) + U(i,j) = Updates(i,j) - S(i,col) + S(i,col) = Updates(i,j) end do end do @@ -55,11 +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, Updates, Updates_index) - !! S_inv_t needs to be transposed back before we - !! can multiply it with S to test unity + call MaponiA3(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") From 6012b9decfc324cc984b7d5b70424502c7386095 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Fri, 12 Mar 2021 15:47:47 +0100 Subject: [PATCH 2/3] 4x4 with 2 updates example is working again after adding update cycle 8169 dataset but with replacement updates instead of additive updates. --- datasets/update_cycle_8169_repl.dat | 32 +++++++++++++++++++++++++++++ src/SM_MaponiA3.cpp | 9 +++++--- tests/QMCChem_dataset_test.f90 | 2 +- 3 files changed, 39 insertions(+), 4 deletions(-) create mode 100644 datasets/update_cycle_8169_repl.dat diff --git a/datasets/update_cycle_8169_repl.dat b/datasets/update_cycle_8169_repl.dat new file mode 100644 index 0000000..46c8522 --- /dev/null +++ b/datasets/update_cycle_8169_repl.dat @@ -0,0 +1,32 @@ +#START_PACKET +#CYCLE_ID: 8169 +#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 +(01,01) 0.100000000000000E+01 0.100000000000000E+01 +(01,02) 0.000000000000000E+00 0.100000000000000E+01 +(01,03) 0.100000000000000E+01 -0.100000000000000E+01 +(01,04) -0.100000000000000E+01 -0.100000000000000E+01 +(02,01) 0.000000000000000E+00 -0.100000000000000E+01 +(02,02) 0.100000000000000E+01 0.000000000000000E+00 +(02,03) 0.100000000000000E+01 0.100000000000000E+01 +(02,04) 0.000000000000000E+00 0.000000000000000E+00 +(03,01) -0.100000000000000E+01 0.100000000000000E+01 +(03,02) 0.000000000000000E+00 0.200000000000000E+01 +(03,03) -0.100000000000000E+01 -0.200000000000000E+01 +(03,04) 0.000000000000000E+00 -0.100000000000000E+01 +(04,01) 0.100000000000000E+01 0.100000000000000E+01 +(04,02) 0.100000000000000E+01 0.100000000000000E+01 +(04,03) 0.100000000000000E+01 -0.100000000000000E+01 +(04,04) 0.100000000000000E+01 0.000000000000000E+00 +#COL_UPDATE_INDEX: 2 +#COL_UPDATE_COMP_(01): 0.000000000000000E+00 +#COL_UPDATE_COMP_(02): -1.000000000000000E+00 +#COL_UPDATE_COMP_(03): 0.000000000000000E+00 +#COL_UPDATE_COMP_(04): 1.000000000000000E+00 +#COL_UPDATE_INDEX: 4 +#COL_UPDATE_COMP_(01): -1.000000000000000E+00 +#COL_UPDATE_COMP_(02): -1.000000000000000E+00 +#COL_UPDATE_COMP_(03): 0.000000000000000E+00 +#COL_UPDATE_COMP_(04): -1.000000000000000E+00 +#END_PACKET diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp index ed64470..786da30 100644 --- a/src/SM_MaponiA3.cpp +++ b/src/SM_MaponiA3.cpp @@ -17,7 +17,7 @@ void selectBestUpdate(unsigned int l, unsigned int N_updates, breakdown = abs(1 + ylk[l][index][component]); #ifdef DEBUG cout << "Inside selectBestUpdate()" << endl; - cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "])" << endl; + cout << "breakdown = abs(1 + ylk[" << l << "][" << index << "][" << component << "]) = " << breakdown << endl; cout << endl; #endif if (breakdown > max) { @@ -69,7 +69,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, for (k = 1; k < N_updates + 1; k++) { #ifdef DEBUG cout << "Compute y0k: " << endl; - cout << "ylk[0][" << k << "][:]" << endl; + cout << "ylk[0][" << k << "][:]"; cout << endl; #endif for (i = 1; i < Dim + 1; i++) { @@ -78,6 +78,9 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, * Updates[(k-1)*Dim + (j-1)]; } } + #ifdef DEBUG + showVector(ylk[0][k], Dim, ""); + #endif } // Calculate the {y_{l,k}} from the {y_{0,k}} @@ -96,7 +99,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim, #ifdef DEBUG cout << "p[l+1] = " << p[l+1] << endl; cout << "component = " << component << endl; - cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "]" << endl; + cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << endl; cout << endl; #endif if (fabs(beta) < 1e-6) { diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90 index 2f537a8..e7b507d 100644 --- a/tests/QMCChem_dataset_test.f90 +++ b/tests/QMCChem_dataset_test.f90 @@ -9,7 +9,7 @@ program QMCChem_dataset_test real(c_double), dimension(:,:), allocatable :: Updates, U real(c_double), dimension(:,:), allocatable :: S, S_inv, S_inv_t - call Read_dataset("datasets/update_cycle_8169.dat", & + call Read_dataset("datasets/update_cycle_8169_repl.dat", & cycle_id, & dim, & n_updates, & From 4e7c4de0e85576ba40369c7b34aed447d42878c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Coppens?= Date: Fri, 12 Mar 2021 16:18:03 +0100 Subject: [PATCH 3/3] tests/test_internal_h5 and test_external_h5 updated to use replacement updates instead of additive updates. --- tests/QMCChem_dataset_test.f90 | 2 +- tests/test_external_h5.cpp | 10 +++++++--- tests/test_internal_h5.cpp | 11 +++++++---- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/tests/QMCChem_dataset_test.f90 b/tests/QMCChem_dataset_test.f90 index e7b507d..a3811db 100644 --- a/tests/QMCChem_dataset_test.f90 +++ b/tests/QMCChem_dataset_test.f90 @@ -77,5 +77,5 @@ program QMCChem_dataset_test close(4000) close(5000) - deallocate(S, S_inv, S_inv_t, Updates, Updates_index) + deallocate(S, S_inv, S_inv_t, Updates, U, Updates_index) end program diff --git a/tests/test_external_h5.cpp b/tests/test_external_h5.cpp index f775510..44c4a5d 100644 --- a/tests/test_external_h5.cpp +++ b/tests/test_external_h5.cpp @@ -47,6 +47,9 @@ int test_cycle(H5File file, int cycle) { double * updates = new double[nupdates*dim]; read_double(file, group + "/updates", updates); + double * u = new double[nupdates*dim]; + + /* Test */ #ifdef DEBUG showMatrix(slater_matrix, dim, "OLD Slater"); @@ -59,12 +62,13 @@ int test_cycle(H5File file, int cycle) { for (j = 0; j < nupdates; j++) { for (i = 0; i < dim; i++) { col = col_update_index[j]; - slater_matrix[i*dim + (col - 1)] += updates[i + j*dim]; + u[i + j*dim] = updates[i + j*dim] - slater_matrix[i*dim + (col - 1)]; + slater_matrix[i*dim + (col - 1)] = updates[i + j*dim]; } } //MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index); - SM(slater_inverse, dim, nupdates, updates, col_update_index); + SM(slater_inverse, dim, nupdates, u, col_update_index); #ifdef DEBUG showMatrix(slater_matrix, dim, "NEW Slater"); @@ -82,7 +86,7 @@ int test_cycle(H5File file, int cycle) { showMatrix(res, dim, "Result"); #endif - delete [] res, updates, col_update_index, + delete [] res, updates, u, col_update_index, slater_matrix, slater_inverse; return ok; diff --git a/tests/test_internal_h5.cpp b/tests/test_internal_h5.cpp index 1ed3075..5d55804 100644 --- a/tests/test_internal_h5.cpp +++ b/tests/test_internal_h5.cpp @@ -8,7 +8,7 @@ #include "Helpers.hpp" using namespace H5; -//#define DEBUG +// #define DEBUG const H5std_string FILE_NAME( "datasets.hdf5" ); @@ -47,6 +47,8 @@ int test_cycle(H5File file, int cycle) { double * updates = new double[nupdates*dim]; read_double(file, group + "/updates", updates); + double * u = new double[nupdates*dim]; + /* Test */ #ifdef DEBUG showMatrix(slater_matrix, dim, "OLD Slater"); @@ -59,11 +61,12 @@ int test_cycle(H5File file, int cycle) { for (j = 0; j < nupdates; j++) { for (i = 0; i < dim; i++) { col = col_update_index[j]; - slater_matrix[i*dim + (col - 1)] += updates[i + j*dim]; + u[i + j*dim] = updates[i + j*dim] - slater_matrix[i*dim + (col - 1)]; + slater_matrix[i*dim + (col - 1)] = updates[i + j*dim]; } } - MaponiA3(slater_inverse, dim, nupdates, updates, col_update_index); + MaponiA3(slater_inverse, dim, nupdates, u, col_update_index); //SM(slater_inverse, dim, nupdates, updates, col_update_index); #ifdef DEBUG @@ -82,7 +85,7 @@ int test_cycle(H5File file, int cycle) { showMatrix(res, dim, "Result"); #endif - delete [] res, updates, col_update_index, + delete [] res, updates, u, col_update_index, slater_matrix, slater_inverse; return ok;