From 1f3a08fa3027268bad45d0a1ff9113803a3d0cd2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 Aug 2023 14:24:33 +0200 Subject: [PATCH] Avoid FPE in Jastrow GL --- org/qmckl_ao.org | 4 +++- org/qmckl_distance.org | 42 +++++++++++++++++++++++++------------ org/qmckl_jastrow_champ.org | 11 +++++++++- 3 files changed, 42 insertions(+), 15 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index a177779..4f419a6 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -5483,7 +5483,9 @@ integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) if (L(3,j) > 1) then w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2) end if - if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return + if (w /= 0.d0) then + if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return + endif end do test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS diff --git a/org/qmckl_distance.org b/org/qmckl_distance.org index df0e15b..a943cad 100644 --- a/org/qmckl_distance.org +++ b/org/qmckl_distance.org @@ -86,7 +86,7 @@ int main() { const double* B, const int64_t ldb, double* const C, - const int64_t ldc ); + const int64_t ldc ); #+end_src #+begin_src f90 :tangle (eval f) @@ -485,7 +485,7 @@ end function test_qmckl_distance_sq const double* B, const int64_t ldb, double* const C, - const int64_t ldc ); + const int64_t ldc ); #+end_src *** Source @@ -909,7 +909,7 @@ end function test_qmckl_dist const int64_t ldb, double* const C, const int64_t ldc, - const double rescale_factor_kappa ); + const double rescale_factor_kappa ); #+end_src *** Source @@ -1215,7 +1215,7 @@ integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C) integer*8, parameter :: elec_num = 10_8 integer*8, parameter :: nucl_num = 2_8 - + double precision :: nucl_coord(nucl_num,3) = reshape( (/ & 0.0d0, 0.0d0 , & 0.0d0, 0.0d0 , & @@ -1263,19 +1263,19 @@ integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C) 1.20919677d0, 1.3111856d0, 1.26122875d0, 0.22122563d0, 0.55669168d0 /), shape(ref_en) ) double precision, allocatable :: distance_ee(:,:), distance_en(:,:) - + allocate( distance_ee(elec_num,elec_num), distance_en(elec_num,nucl_num) ) - + print *, 'ee' test_qmckl_dist_rescaled = & qmckl_distance_rescaled(context, 'N', 'N', elec_num, elec_num, elec_coord, & size(elec_coord,1)*1_8, elec_coord, size(elec_coord,1)*1_8, & distance_ee, size(distance_ee,1)*1_8, kappa) - + if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return - + test_qmckl_dist_rescaled = QMCKL_FAILURE - + do j=1,elec_num do i=1,elec_num print *, i,j,real(distance_ee(i,j)), real(ref_ee(i,j)) @@ -1294,7 +1294,7 @@ integer(qmckl_exit_code) function test_qmckl_dist_rescaled(context) bind(C) if (test_qmckl_dist_rescaled /= QMCKL_SUCCESS) return test_qmckl_dist_rescaled = QMCKL_FAILURE - + do j=1,nucl_num do i=1,elec_num print *, i,j,real(distance_en(i,j)), real(ref_en(i,j)) @@ -1385,7 +1385,7 @@ end function test_qmckl_dist_rescaled - ~A~ is allocated with at least $3 \times m \times 8$ bytes - ~B~ is allocated with at least $3 \times n \times 8$ bytes - ~C~ is allocated with at least $4 \times m \times n \times 8$ bytes - + #+CALL: generate_c_header(table=qmckl_distance_rescaled_gl_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: @@ -1402,7 +1402,7 @@ end function test_qmckl_dist_rescaled const int64_t ldb, double* const C, const int64_t ldc, - const double rescale_factor_kappa ); + const double rescale_factor_kappa ); #+end_src #+begin_src f90 :tangle (eval f) @@ -1504,6 +1504,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, & y = A(2,i) - B(2,j) z = A(3,i) - B(3,j) dist = dsqrt(x*x + y*y + z*z) + ! Avoid floating-point exception + if (dist == 0.d0) then + dist = 69.d0/rescale_factor_kappa + endif dist_inv = 1.0d0/dist rij = dexp(-rescale_factor_kappa * dist) C(1,i,j) = x * dist_inv * rij @@ -1521,6 +1525,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, & y = A(i,2) - B(2,j) z = A(i,3) - B(3,j) dist = dsqrt(x*x + y*y + z*z) + ! Avoid floating-point exception + if (dist == 0.d0) then + dist = 69.d0/rescale_factor_kappa + endif dist_inv = 1.0d0/dist rij = dexp(-rescale_factor_kappa * dist) C(1,i,j) = x * dist_inv * rij @@ -1538,6 +1546,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, & y = A(2,i) - B(j,2) z = A(3,i) - B(j,3) dist = dsqrt(x*x + y*y + z*z) + ! Avoid floating-point exception + if (dist == 0.d0) then + dist = 69.d0/rescale_factor_kappa + endif dist_inv = 1.0d0/dist rij = dexp(-rescale_factor_kappa * dist) C(1,i,j) = x * dist_inv * rij @@ -1555,6 +1567,10 @@ integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, & y = A(i,2) - B(j,2) z = A(i,3) - B(j,3) dist = dsqrt(x*x + y*y + z*z) + ! Avoid floating-point exception + if (dist == 0.d0) then + dist = 69.d0/rescale_factor_kappa + endif dist_inv = 1.0d0/dist rij = dexp(-rescale_factor_kappa * dist) C(1,i,j) = x * dist_inv * rij @@ -1570,7 +1586,7 @@ end function qmckl_distance_rescaled_gl_f #+end_src This function is more efficient when ~A~ and ~B~ are transposed. - + #+CALL: generate_c_interface(table=qmckl_distance_rescaled_gl_args,fname=get_value("Name")) #+RESULTS: diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org index 9a34500..3484f1d 100644 --- a/org/qmckl_jastrow_champ.org +++ b/org/qmckl_jastrow_champ.org @@ -6095,14 +6095,23 @@ integer function qmckl_compute_jastrow_champ_factor_een_rescaled_e_gl_f( & een_rescaled_e_gl = 0.0d0 do nw = 1, walk_num do j = 1, elec_num - do i = 1, elec_num + do i = 1, j-1 rij_inv = 1.0d0 / ee_distance(i, j, nw) do ii = 1, 3 elec_dist_gl(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv end do elec_dist_gl(4, i, 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(ii, i, j) = (coord_ee(i, ii, nw) - coord_ee(j, ii, nw)) * rij_inv + end do + elec_dist_gl(4, i, j) = 2.0d0 * rij_inv + end do end do ! prepare the actual een table