From 71de52c39fdcad574e3289fb8bb743e52f04feb7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 13 Dec 2024 03:00:01 +0100 Subject: [PATCH] Introduced many tests in Jastrow --- configure.ac | 7 ++ org/qmckl_jastrow_champ.org | 218 ++++++++++++++++++++++++++---------- 2 files changed, 168 insertions(+), 57 deletions(-) diff --git a/configure.ac b/configure.ac index 9ad7633..ae01256 100644 --- a/configure.ac +++ b/configure.ac @@ -359,6 +359,13 @@ AS_IF([test "x$ok" = "xyes"], [ ]) +AC_ARG_ENABLE(sanitizer, [AS_HELP_STRING([--enable-sanitizer],[enable sanitizer debug flags])], ok=$enableval, ok=no) +if test "$ok" = "yes"; then + CFLAGS="${CFLAGS} -fsanitize=address -fsanitize=undefined -fsanitize=leak -fsanitize=pointer-compare -fsanitize=pointer-subtract -fsanitize=bounds -fsanitize=bounds-strict" + FCFLAGS="${FCFLAGS} -fsanitize=address -fsanitize=undefined -fsanitize=leak -fsanitize=pointer-compare -fsanitize=pointer-subtract -fsanitize=bounds -fsanitize=bounds-strict" +fi + + ## diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org index ec6bb97..43a497a 100644 --- a/org/qmckl_jastrow_champ.org +++ b/org/qmckl_jastrow_champ.org @@ -93,6 +93,7 @@ #include #include #include "n2.h" +#include "qmckl_context_private_type.h" #include "qmckl_jastrow_champ_private_func.h" int main() { @@ -2554,7 +2555,7 @@ assert(qmckl_electron_provided(context)); double ee_distance_rescaled_hpc[walk_num * elec_num * elec_num * (cord_num+1)]; memset(ee_distance_rescaled_hpc, 0, sizeof(ee_distance_rescaled_hpc)); - rc = qmckl_compute_een_rescaled_e_doc (context, walk_num, + rc = qmckl_compute_een_rescaled_e_hpc (context, walk_num, elec_num, cord_num, rescale_factor_ee, &(ee_distance[0]), @@ -7183,7 +7184,7 @@ print(" een_rescaled_e[1, 5, 2] = ",een_rescaled_e[1, 5, 2]) #+begin_src c :tangle (eval c_test) assert(qmckl_electron_provided(context)); - +{ double een_rescaled_e[walk_num][(cord_num + 1)][elec_num][elec_num]; rc = qmckl_get_jastrow_champ_een_distance_rescaled_e(context, &(een_rescaled_e[0][0][0][0]),elec_num*elec_num*(cord_num+1)*walk_num); @@ -7195,6 +7196,36 @@ assert(fabs(een_rescaled_e[0][1][0][4]- 0.0884012350763747 ) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][3]- 0.1016685507354656 ) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][4]- 0.0113118073246869 ) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][5]- 0.5257156022077619 ) < 1.e-12); +} + +{ + printf("een_rescaled_e_hpc\n"); + double ee_distance[walk_num * elec_num * elec_num]; + rc = qmckl_get_electron_ee_distance(context, &(ee_distance[0]), walk_num*elec_num*elec_num); + assert(rc == QMCKL_SUCCESS); + + double een_rescaled_e_doc[walk_num][cord_num+1][elec_num][elec_num]; + memset(&(een_rescaled_e_doc[0][0][0][0]), 0, sizeof(een_rescaled_e_doc)); + rc = qmckl_compute_een_rescaled_e(context, walk_num, elec_num, cord_num, + 0.6, &(ee_distance[0]), &(een_rescaled_e_doc[0][0][0][0])); + assert(rc == QMCKL_SUCCESS); + + double een_rescaled_e_hpc[walk_num][cord_num+1][elec_num][elec_num]; + memset(&(een_rescaled_e_hpc[0][0][0][0]), 0, sizeof(een_rescaled_e_hpc)); + rc = qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, + 0.6, &(ee_distance[0]), &(een_rescaled_e_hpc[0][0][0][0])); + assert(rc == QMCKL_SUCCESS); + + for (int64_t i = 0; i < walk_num; i++) { + for (int64_t j = 0; j < cord_num+1; j++) { + for (int64_t k = 0; k < elec_num; k++) { + for (int64_t l = 0; l < elec_num; l++) { + assert(fabs(een_rescaled_e_doc[i][j][k][l] - een_rescaled_e_hpc[i][j][k][l]) < 1.e-8); + } + } + } + } +} #+end_src *** Electron-electron rescaled distances derivatives in $J_\text{eeN}$ @@ -7367,10 +7398,12 @@ integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( & double precision , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) double precision , intent(out) :: een_rescaled_e_gl(elec_num,4,elec_num,0:cord_num,walk_num) - double precision,dimension(:,:,:),allocatable :: elec_dist_gl - double precision :: x, rij_inv, kappa_l + double precision, allocatable :: elec_dist_gl(:,:,:) + double precision :: x, kappa_l integer*8 :: i, j, k, l, nw, ii + double precision :: rij_inv(elec_num) + allocate(elec_dist_gl(elec_num, 4, elec_num)) info = QMCKL_SUCCESS @@ -7398,22 +7431,15 @@ integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( & ! Prepare table of exponentiated distances raised to appropriate power do nw = 1, walk_num do j = 1, elec_num - do i = 1, j-1 - rij_inv = 1.0d0 / ee_distance(i, j, nw) + do i = 1, elec_num + rij_inv(i) = 1.0d0 / (ee_distance(i, j, nw) + 1.d-30) + enddo + rij_inv(j) = 0.0d0 + do i = 1, elec_num do ii = 1, 3 - elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv + elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv(i) end do - elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv - end do - - elec_dist_gl(j, :, j) = 0.0d0 - - do i = j+1, elec_num - rij_inv = 1.0d0 / ee_distance(i, j, nw) - do ii = 1, 3 - elec_dist_gl(i, ii, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv - end do - elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv + elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv(i) end do end do @@ -7428,17 +7454,21 @@ integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( & een_rescaled_e_gl(i, 2, j, l, nw) = kappa_l * elec_dist_gl(i, 2, j) een_rescaled_e_gl(i, 3, j, l, nw) = kappa_l * elec_dist_gl(i, 3, j) een_rescaled_e_gl(i, 4, j, l, nw) = kappa_l * elec_dist_gl(i, 4, j) + end do - een_rescaled_e_gl(i, 4, j, l, nw) = een_rescaled_e_gl(i, 4, j, l, nw) & - + een_rescaled_e_gl(i, 1, j, l, nw) * een_rescaled_e_gl(i, 1, j, l, nw) & - + een_rescaled_e_gl(i, 2, j, l, nw) * een_rescaled_e_gl(i, 2, j, l, nw) & - + een_rescaled_e_gl(i, 3, j, l, nw) * een_rescaled_e_gl(i, 3, j, l, nw) - - een_rescaled_e_gl(i,1,j,l,nw) = een_rescaled_e_gl(i,1,j,l,nw) * een_rescaled_e(i,j,l,nw) - een_rescaled_e_gl(i,2,j,l,nw) = een_rescaled_e_gl(i,2,j,l,nw) * een_rescaled_e(i,j,l,nw) - een_rescaled_e_gl(i,3,j,l,nw) = een_rescaled_e_gl(i,3,j,l,nw) * een_rescaled_e(i,j,l,nw) - een_rescaled_e_gl(i,4,j,l,nw) = een_rescaled_e_gl(i,4,j,l,nw) * een_rescaled_e(i,j,l,nw) + do i = 1, elec_num + een_rescaled_e_gl(i, 4, j, l, nw) = een_rescaled_e_gl(i, 4, j, l, nw) & + + een_rescaled_e_gl(i, 1, j, l, nw) * een_rescaled_e_gl(i, 1, j, l, nw) & + + een_rescaled_e_gl(i, 2, j, l, nw) * een_rescaled_e_gl(i, 2, j, l, nw) & + + een_rescaled_e_gl(i, 3, j, l, nw) * een_rescaled_e_gl(i, 3, j, l, nw) end do + + do i = 1, elec_num + een_rescaled_e_gl(i,1,j,l,nw) = een_rescaled_e_gl(i,1,j,l,nw) * een_rescaled_e(i,j,l,nw) + een_rescaled_e_gl(i,2,j,l,nw) = een_rescaled_e_gl(i,2,j,l,nw) * een_rescaled_e(i,j,l,nw) + een_rescaled_e_gl(i,3,j,l,nw) = een_rescaled_e_gl(i,3,j,l,nw) * een_rescaled_e(i,j,l,nw) + een_rescaled_e_gl(i,4,j,l,nw) = een_rescaled_e_gl(i,4,j,l,nw) * een_rescaled_e(i,j,l,nw) + end do end do end do end do @@ -7606,7 +7636,6 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc ( double kappa_l = - (double)l * rescale_factor_ee; for (int64_t j = 0; j < elec_num; ++j) { double* restrict eegl = &een_rescaled_e_gl[ elec_num * 4 * (j + elec_num * (l + (cord_num + 1) * nw))]; - const double* restrict ee = &een_rescaled_e [ elec_num * (j + elec_num * (l + (cord_num + 1) * nw))]; #ifdef HAVE_OPENMP #pragma omp simd #endif @@ -7625,6 +7654,7 @@ qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc ( eegl[i + elec_num*1] * eegl[i + elec_num*1] + eegl[i + elec_num*2] * eegl[i + elec_num*2]; } + const double* restrict ee = &een_rescaled_e[ elec_num * (j + elec_num * (l + (cord_num + 1) * nw))]; #ifdef HAVE_OPENMP #pragma omp simd #endif @@ -7743,27 +7773,62 @@ print(" een_rescaled_e_gl[2, 1, 6, 2] = ",een_rescaled_e_gl[1, 0, 5, 2]) #+end_src #+begin_src c :tangle (eval c_test) -double een_rescaled_e_gl[walk_num][(cord_num + 1)][elec_num][4][elec_num]; -size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num; -rc = qmckl_get_jastrow_champ_een_distance_rescaled_e_gl(context, - &(een_rescaled_e_gl[0][0][0][0][0]),size_max); +{ + double een_rescaled_e_gl[walk_num][(cord_num + 1)][elec_num][4][elec_num]; + size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num; + rc = qmckl_get_jastrow_champ_een_distance_rescaled_e_gl(context, + &(een_rescaled_e_gl[0][0][0][0][0]),size_max); -assert(fabs(een_rescaled_e_gl[0][1][0][0][2] + 0.09831391870751387 ) < 1.e-12); -assert(fabs(een_rescaled_e_gl[0][1][0][0][3] + 0.017204157459682526 ) < 1.e-12); -assert(fabs(een_rescaled_e_gl[0][1][0][0][4] + 0.013345768421098641 ) < 1.e-12); -assert(fabs(een_rescaled_e_gl[0][2][1][0][3] + 0.03733086358273962 ) < 1.e-12); -assert(fabs(een_rescaled_e_gl[0][2][1][0][4] + 0.004922634822943517 ) < 1.e-12); -assert(fabs(een_rescaled_e_gl[0][2][1][0][5] + 0.5416751547830984 ) < 1.e-12); - #+end_src + assert(fabs(een_rescaled_e_gl[0][1][0][0][2] + 0.09831391870751387 ) < 1.e-12); + assert(fabs(een_rescaled_e_gl[0][1][0][0][3] + 0.017204157459682526 ) < 1.e-12); + assert(fabs(een_rescaled_e_gl[0][1][0][0][4] + 0.013345768421098641 ) < 1.e-12); + assert(fabs(een_rescaled_e_gl[0][2][1][0][3] + 0.03733086358273962 ) < 1.e-12); + assert(fabs(een_rescaled_e_gl[0][2][1][0][4] + 0.004922634822943517 ) < 1.e-12); + assert(fabs(een_rescaled_e_gl[0][2][1][0][5] + 0.5416751547830984 ) < 1.e-12); + #+end_src - #+begin_src c :tangle (eval c_test) -assert(qmckl_electron_provided(context)); + #+begin_src c :tangle (eval c_test) + assert(qmckl_electron_provided(context)); + + qmckl_context_struct* ctx = (qmckl_context_struct*) context; + double een_rescaled_e_gl_doc[walk_num*(cord_num + 1)*elec_num*4*elec_num]; + memset(een_rescaled_e_gl_doc, 0, sizeof(een_rescaled_e_gl_doc)); + rc = qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_doc(context, + ctx->electron.walker.num, + ctx->electron.num, + ctx->jastrow_champ.cord_num, + ctx->jastrow_champ.rescale_factor_ee, + ctx->electron.walker.point.coord.data, + ctx->electron.ee_distance, + ctx->jastrow_champ.een_rescaled_e, + een_rescaled_e_gl_doc); + assert(rc == QMCKL_SUCCESS); + + double een_rescaled_e_gl_hpc[walk_num*(cord_num + 1)*elec_num*4*elec_num]; + memset(een_rescaled_e_gl_hpc, 0, sizeof(een_rescaled_e_gl_hpc)); + rc = qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_hpc(context, + ctx->electron.walker.num, + ctx->electron.num, + ctx->jastrow_champ.cord_num, + ctx->jastrow_champ.rescale_factor_ee, + ctx->electron.walker.point.coord.data, + ctx->electron.ee_distance, + ctx->jastrow_champ.een_rescaled_e, + een_rescaled_e_gl_hpc); + assert(rc == QMCKL_SUCCESS); + + for (int64_t i = 0; i < walk_num*(cord_num + 1)*elec_num*4*elec_num; i++) { + printf("i = %ld, doc = %e, hpc = %e\n", i, een_rescaled_e_gl_doc[i], een_rescaled_e_gl_hpc[i]); + assert(fabs(een_rescaled_e_gl_doc[i] - een_rescaled_e_gl_hpc[i]) < 1.e-12); + } + +} { printf("een_distance_rescaled_e_gl\n"); double fd[walk_num][cord_num+1][elec_num][4][elec_num]; - double delta_x = 0.0001; + double delta_x = 0.001; // Finite difference coefficients for gradients double coef[9] = { 1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0 }; @@ -7795,8 +7860,6 @@ assert(qmckl_electron_provided(context)); memset(&(fd[0][0][0][0]), 0, sizeof(fd)); - printf("%lu %lu\n", sizeof(function_values)/sizeof(double), walk_num*(cord_num+1)*elec_num*elec_num); - for (int64_t k = 0; k < 3; k++) { for (int64_t i = 0; i < elec_num; i++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement @@ -7810,6 +7873,8 @@ assert(qmckl_electron_provided(context)); &(temp_coord[0][0][0]), walk_num*3*elec_num); assert(rc == QMCKL_SUCCESS); + rc = qmckl_context_touch(context); + assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_jastrow_champ_een_distance_rescaled_e(context, @@ -7827,9 +7892,9 @@ assert(qmckl_electron_provided(context)); } } - for (int64_t nw=0 ; nwelectron.walker.num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow_champ.cord_num, + ctx->jastrow_champ.dim_c_vector, + ctx->jastrow_champ.c_vector_full, + ctx->jastrow_champ.lkpm_combined_index, + ctx->jastrow_champ.tmp_c, + ctx->jastrow_champ.dtmp_c, + ctx->jastrow_champ.een_rescaled_n, + ctx->jastrow_champ.een_rescaled_n_gl, + factor_een_gl_doc); + assert(rc == QMCKL_SUCCESS); + + double factor_een_gl_hpc[walk_num*4*elec_num]; + memset(&(factor_een_gl_hpc[0]), 0, walk_num*4*elec_num*sizeof(double)); + rc = qmckl_compute_jastrow_champ_factor_een_gl(context, + ctx->electron.walker.num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow_champ.cord_num, + ctx->jastrow_champ.dim_c_vector, + ctx->jastrow_champ.c_vector_full, + ctx->jastrow_champ.lkpm_combined_index, + ctx->jastrow_champ.tmp_c, + ctx->jastrow_champ.dtmp_c, + ctx->jastrow_champ.een_rescaled_n, + ctx->jastrow_champ.een_rescaled_n_gl, + factor_een_gl_hpc); + + for (int64_t i = 0; i < walk_num*4*elec_num; ++i) { + assert(fabs(factor_een_gl_doc[i] - factor_een_gl_hpc[i]) < 1e-12); + } +} #+end_src ** Total Jastrow