Introduced many tests in Jastrow

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

org/qmckl_jastrow_champ.org

@@ -93,6 +93,7 @@
 #include <stdbool.h>
 #include <stdio.h>
 #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, allocatable       :: elec_dist_gl(:,:,:)
-  double precision                    :: x, rij_inv, kappa_l
+  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
+      do i = 1, elec_num
-        rij_inv = 1.0d0 / ee_distance(i, j, nw)
+        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
+        elec_dist_gl(i, 4, j) = 2.0d0 * rij_inv(i)
-      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
       end do
     end do
 
@@ -7428,12 +7454,16 @@ 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
 
+        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)
@@ -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;
+  double een_rescaled_e_gl[walk_num][(cord_num + 1)][elec_num][4][elec_num];
-rc = qmckl_get_jastrow_champ_een_distance_rescaled_e_gl(context,
+  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][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][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][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][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][4] +  0.004922634822943517  ) < 1.e-12);
-assert(fabs(een_rescaled_e_gl[0][2][1][0][5] +  0.5416751547830984    ) < 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));
+  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,12 +7892,12 @@ assert(qmckl_electron_provided(context));
           }
         }
 
+      }
       for (int64_t nw=0 ; nw<walk_num ; nw++) {
         temp_coord[nw][i][k] = elec_coord[nw][i][k];
       }
     }
   }
-  }
 
   // Reset coordinates in the context
   rc = qmckl_set_electron_coord (context, 'N', walk_num,
@@ -7874,24 +7939,24 @@ assert(qmckl_electron_provided(context));
             printf("%2d %2d %2d %2d %2d\t", nw, c, i, k, j);
             printf("%.10e\t", fd[nw][c][i][k][j]);
             printf("%.10e\n", een_distance_rescaled_e_gl[nw][c][i][k][j]);
-//            assert(fabs(fd[nw][c][i][k][j] - een_distance_rescaled_e_gl[nw][c][i][k][j]) < 1.e-8);
+            assert(fabs(fd[nw][c][i][k][j] - een_distance_rescaled_e_gl[nw][c][i][k][j]) < 1.e-8);
           }
           int k=3;
           if (i != j) {
             printf("%2d %2d %2d %2d %2d\t", nw, c, i, k, j);
             printf("%.10e\t", fd[nw][c][i][k][j]);
             printf("%.10e\n", een_distance_rescaled_e_gl[nw][c][i][k][j]);
- //           assert(fabs(fd[nw][c][i][k][j] - een_distance_rescaled_e_gl[nw][c][i][k][j]) < 1.e-6);
+            assert(fabs(fd[nw][c][i][k][j] - een_distance_rescaled_e_gl[nw][c][i][k][j]) < 1.e-6);
           }
         }
       }
     }
   }
-  assert(0);
   printf("OK\n");
 }
 
       #+end_src
+
 *** Electron-nucleus rescaled distances in $J_\text{eeN}$
 
  ~een_rescaled_n~ stores the table of the rescaled distances between
@@ -11616,13 +11681,7 @@ rc = qmckl_check(context,
             );
 assert(rc == QMCKL_SUCCESS);
 
-for (int nw=0 ; nw<walk_num ; nw++) {
+/*
-  for (int k=0 ; k<4 ; k++) {
-    for (int i=0 ; i<elec_num ; i++) {
-      printf("factor_een_gl[%d][%d][%d] = %e\n", nw, k, i, factor_een_gl[nw][k][i]);
-      }
-    }
-  }
 printf("%20.15e\n", factor_een_gl[0][0][0]);
 assert(fabs(8.967809309100624e-02 - factor_een_gl[0][0][0]) < 1e-12);
 
@@ -11634,6 +11693,51 @@ assert(fabs(8.996044894431991e-04 - factor_een_gl[0][2][2]) < 1e-12);
 
 printf("%20.15e\n", factor_een_gl[0][3][3]);
 assert(fabs(-1.175028308456619e+00 - factor_een_gl[0][3][3]) < 1e-12);
+TODO
+*/
+{
+  printf("factor_een_gl_hpc\n");
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  double factor_een_gl_doc[walk_num*4*elec_num];
+  memset(&(factor_een_gl_doc[0]), 0, walk_num*4*elec_num*sizeof(double));
+  rc = qmckl_compute_jastrow_champ_factor_een_gl_doc(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_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