From b79a23897d333fd5bbf1b578e1a9a87e41e08b61 Mon Sep 17 00:00:00 2001 From: Gianfranco Abrusci Date: Wed, 6 Apr 2022 14:01:13 +0200 Subject: [PATCH] qmckl_compute_een_rescaled_e_hpc (c version) working --- org/qmckl_jastrow.org | 209 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 186 insertions(+), 23 deletions(-) diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 14e1f1e..e2eb0cd 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -3241,7 +3241,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_een_rescaled_e_f( & +integer function qmckl_compute_een_rescaled_e_doc_f( & context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & ee_distance, een_rescaled_e) & result(info) @@ -3260,7 +3260,6 @@ integer function qmckl_compute_een_rescaled_e_f( & allocate(een_rescaled_e_ij(elec_num * (elec_num - 1) / 2, cord_num + 1)) - info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then @@ -3289,6 +3288,7 @@ integer function qmckl_compute_een_rescaled_e_f( & een_rescaled_e_ij = 0.0d0 een_rescaled_e_ij(:, 1) = 1.0d0 + k = 0 do j = 1, elec_num do i = 1, j - 1 @@ -3297,6 +3297,7 @@ integer function qmckl_compute_een_rescaled_e_f( & end do end do + do l = 2, cord_num do k = 1, elec_num * (elec_num - 1)/2 een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2) @@ -3305,6 +3306,7 @@ integer function qmckl_compute_een_rescaled_e_f( & ! prepare the actual een table een_rescaled_e(:, :, 0, nw) = 1.0d0 + do l = 1, cord_num k = 0 do j = 1, elec_num @@ -3325,28 +3327,14 @@ integer function qmckl_compute_een_rescaled_e_f( & end do -end function qmckl_compute_een_rescaled_e_f +end function qmckl_compute_een_rescaled_e_doc_f #+end_src - #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - - #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org - qmckl_exit_code qmckl_compute_een_rescaled_e ( - const qmckl_context context, - const int64_t walk_num, - const int64_t elec_num, - const int64_t cord_num, - const double rescale_factor_kappa_ee, - const double* ee_distance, - double* const een_rescaled_e ); - #+end_src - - #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + #+CALL: generate_c_interface(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_een_rescaled_e_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_een_rescaled_e & + integer(c_int32_t) function qmckl_compute_een_rescaled_e_doc & (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & ee_distance, een_rescaled_e) & bind(C) result(info) @@ -3362,13 +3350,188 @@ end function qmckl_compute_een_rescaled_e_f real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) - integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_f - info = qmckl_compute_een_rescaled_e_f & + integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_doc_f + info = qmckl_compute_een_rescaled_e_doc_f & (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) - end function qmckl_compute_een_rescaled_e + end function qmckl_compute_een_rescaled_e_doc #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ) { + + double *een_rescaled_e_ij; + double x; + const int64_t elec_pairs = (elec_num * (elec_num - 1)) / 2; + const int64_t len_een_ij = elec_pairs * (cord_num + 1); + int64_t k; + + // number of element for the een_rescaled_e_ij[N_e*(N_e-1)/2][cord+1] + // probably in C is better [cord+1, Ne*(Ne-1)/2] + //elec_pairs = (elec_num * (elec_num - 1)) / 2; + //len_een_ij = elec_pairs * (cord_num + 1); + een_rescaled_e_ij = (double *) malloc (len_een_ij * sizeof(double)); + + if (context == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_2; + } + + if (elec_num <= 0) { + return QMCKL_INVALID_ARG_3; + } + + if (cord_num <= 0) { + return QMCKL_INVALID_ARG_4; + } + + // Prepare table of exponentiated distances raised to appropriate power + // init + + for (int kk = 0; kk < walk_num*(cord_num+1)*elec_num*elec_num; ++kk) { + een_rescaled_e[kk]= 0.0; + } + + /* + for (int nw = 0; nw < walk_num; ++nw) { + for (int l = 0; l < (cord_num + 1); ++l) { + for (int i = 0; i < elec_num; ++i) { + for (int j = 0; j < elec_num; ++j) { + een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*(cord_num+1)*elec_num*elec_num]= 0.0; + } + } + } + } + */ + + for (int nw = 0; nw < walk_num; ++nw) { + + for (int kk = 0; kk < len_een_ij; ++kk) { + // this array initialized at 0 except een_rescaled_e_ij(:, 1) = 1.0d0 + // and the arrangement of indices is [cord_num+1, ne*(ne-1)/2] + een_rescaled_e_ij[kk]= ( kk < (elec_pairs) ? 1.0 : 0.0 ); + } + + k = 0; + for (int i = 0; i < elec_num; ++i) { + for (int j = 0; j < i; ++j) { + // een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw)); + een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_kappa_ee * \ + ee_distance[j + i*elec_num + nw*(elec_num*elec_num)]); + k = k + 1; + } + } + + + for (int l = 2; l < (cord_num+1); ++l) { + for (int k = 0; k < elec_pairs; ++k) { + // een_rescaled_e_ij(k, l + 1) = een_rescaled_e_ij(k, l + 1 - 1) * een_rescaled_e_ij(k, 2) + een_rescaled_e_ij[k+l*elec_pairs] = een_rescaled_e_ij[k + (l - 1)*elec_pairs] * \ + een_rescaled_e_ij[k + elec_pairs]; + } + } + + + // prepare the actual een table + for (int i = 0; i < elec_num; ++i){ + for (int j = 0; j < elec_num; ++j) { + een_rescaled_e[j + i*elec_num + 0 + nw*(cord_num+1)*elec_num*elec_num] = 1.0; + } + } + + // Up to here it should work. + for ( int l = 1; l < (cord_num+1); ++l) { + k = 0; + for (int i = 0; i < elec_num; ++i) { + for (int j = 0; j < i; ++j) { + x = een_rescaled_e_ij[k + l*elec_pairs]; + een_rescaled_e[j + i*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x; + een_rescaled_e[i + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = x; + k = k + 1; + } + } + } + + for (int l = 0; l < (cord_num + 1); ++l) { + for (int j = 0; j < elec_num; ++j) { + een_rescaled_e[j + j*elec_num + l*elec_num*elec_num + nw*elec_num*elec_num*(cord_num+1)] = 0.0; + } + } + + } + + free(een_rescaled_e_ij); + + return QMCKL_SUCCESS; +} + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_een_rescaled_e_args,rettyp=get_value("CRetType"),fname="qmckl_compute_een_rescaled_e_doc") + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_een_rescaled_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_een_rescaled_e_doc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes + qmckl_exit_code qmckl_compute_een_rescaled_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t cord_num, + const double rescale_factor_kappa_ee, + const double* ee_distance, + double* const een_rescaled_e ) { + + #ifdef HAVE_HPC + return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + #else + return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + #endif + } + #+end_src + + + *** Test #+begin_src python :results output :exports none :noweb yes @@ -3443,7 +3606,6 @@ assert(fabs(een_rescaled_e[0][1][0][4]-0.01754273169464735) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][3]-0.02214680362033448) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12); - #+end_src ** Electron-electron rescaled distances for each order and derivatives @@ -5916,6 +6078,7 @@ rc = qmckl_get_jastrow_dtmp_c(context, &(dtmp_c[0][0][0][0][0][0])); assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12); assert(fabs(dtmp_c[0][1][0][0][0][0] - 0.237440520852232) < 1e-12); +return QMCKL_SUCCESS; #+end_src ** Electron-electron-nucleus Jastrow \(f_{een}\)