Merge branch 'gpu' of github.com:TREX-CoE/qmckl into gpu

org/qmckl_jastrow.org

@@ -1507,7 +1507,6 @@ qmckl_exit_code qmckl_compute_asymp_jasb (
 
 #   #+CALL: generate_c_header(table=qmckl_asymp_jasb_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
 
-    #+RESULTS:
     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
     qmckl_exit_code qmckl_compute_asymp_jasb (
           const qmckl_context context,
@@ -1856,9 +1855,8 @@ qmckl_exit_code qmckl_compute_factor_ee (
 }
 #+end_src
 
-    #+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+#   #+CALL: generate_c_header(table=qmckl_factor_ee_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
     
-    #+RESULTS:
     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
     qmckl_exit_code qmckl_compute_factor_ee (
           const qmckl_context context,
@@ -2169,7 +2167,7 @@ end function qmckl_compute_factor_ee_deriv_e_f
 
 #   #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
 
-    #+RESULTS:
+
     #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
     qmckl_exit_code qmckl_compute_factor_ee_deriv_e (
           const qmckl_context context,
@@ -3242,7 +3240,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)
@@ -3261,7 +3259,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
@@ -3290,6 +3287,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
@@ -3298,6 +3296,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)
@@ -3306,6 +3305,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
@@ -3326,7 +3326,7 @@ 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"))
@@ -3342,13 +3342,14 @@ end function qmckl_compute_een_rescaled_e_f
           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) &
+	(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, &
-        bind(C) result(info)
+        ee_distance, een_rescaled_e) &
+    bind(C) result(info)
 
       use, intrinsic :: iso_c_binding
       implicit none
@@ -3361,13 +3362,186 @@ 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
+      integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_doc_f
-      info = qmckl_compute_een_rescaled_e_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)
+	     (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")
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+    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
@@ -3442,7 +3616,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
@@ -5118,7 +5291,7 @@ qmckl_exit_code qmckl_compute_dim_cord_vect (
     | ~cord_vect_full~   | ~double[dim_cord_vect][nucl_num]~      | out    | Full list of coefficients    |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_cord_vect_full_f( &
+integer function qmckl_compute_cord_vect_full_doc_f( &
      context, nucl_num, dim_cord_vect, type_nucl_num,  &
      type_nucl_vector, cord_vector, cord_vect_full) &
      result(info)
@@ -5161,30 +5334,16 @@ integer function qmckl_compute_cord_vect_full_f( &
     cord_vect_full(a,1:dim_cord_vect) = cord_vector(type_nucl_vector(a),1:dim_cord_vect)
   end do
 
-end function qmckl_compute_cord_vect_full_f
+end function qmckl_compute_cord_vect_full_doc_f
     #+end_src
 
-#   #+CALL: generate_c_header(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+    #+CALL: generate_c_interface(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname="qmckl_compute_cord_vect_full_doc")
-
-    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
-    qmckl_exit_code qmckl_compute_cord_vect_full (
-          const qmckl_context context,
-          const int64_t nucl_num,
-          const int64_t dim_cord_vect,
-          const int64_t type_nucl_num,
-          const int64_t* type_nucl_vector,
-          const double* cord_vector,
-          double* const cord_vect_full ); 
-    #+end_src
-
-
-    #+CALL: generate_c_interface(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
 
     #+RESULTS:
     #+begin_src f90 :tangle (eval f) :comments org :exports none
-    integer(c_int32_t) function qmckl_compute_cord_vect_full &
+    integer(c_int32_t) function qmckl_compute_cord_vect_full_doc &
-        (context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full) &
+	(context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full) &
-        bind(C) result(info)
+	bind(C) result(info)
 
       use, intrinsic :: iso_c_binding
       implicit none
@@ -5197,13 +5356,105 @@ end function qmckl_compute_cord_vect_full_f
       real    (c_double ) , intent(in)          :: cord_vector(type_nucl_num,dim_cord_vect)
       real    (c_double ) , intent(out)         :: cord_vect_full(nucl_num,dim_cord_vect)
 
-      integer(c_int32_t), external :: qmckl_compute_cord_vect_full_f
+      integer(c_int32_t), external :: qmckl_compute_cord_vect_full_doc_f
-      info = qmckl_compute_cord_vect_full_f &
+      info = qmckl_compute_cord_vect_full_doc_f &
-             (context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full)
+	     (context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full)
 
-    end function qmckl_compute_cord_vect_full
+    end function qmckl_compute_cord_vect_full_doc
     #+end_src
 
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_cord_vect_full_hpc (
+	  const qmckl_context context,
+	  const int64_t nucl_num,
+	  const int64_t dim_cord_vect,
+	  const int64_t type_nucl_num,
+	  const int64_t* type_nucl_vector,
+	  const double* cord_vector,
+	  double* const cord_vect_full ) {
+
+  if (context == QMCKL_NULL_CONTEXT) {
+     return QMCKL_INVALID_CONTEXT;
+  }
+
+  if (nucl_num <= 0) {
+     return QMCKL_INVALID_ARG_2;
+  }
+
+  if (type_nucl_num <= 0) {
+     return QMCKL_INVALID_ARG_4;
+  }
+
+  if (dim_cord_vect <= 0) {
+     return QMCKL_INVALID_ARG_5;
+  }
+
+  for (int i=0; i < dim_cord_vect; ++i) {
+    for (int a=0; a < nucl_num; ++a){
+      cord_vect_full[a + i*nucl_num] = cord_vector[(type_nucl_vector[a]-1)+i*type_nucl_num];
+    }
+  }
+
+  return QMCKL_SUCCESS;
+  }
+    #+end_src
+
+
+#   #+CALL: generate_c_header(table=qmckl_factor_cord_vect_full_args,rettyp=get_value("CRetType"),fname="qmckl_compute_cord_vect_full_doc")
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+    qmckl_exit_code qmckl_compute_cord_vect_full (
+	  const qmckl_context context,
+	  const int64_t nucl_num,
+	  const int64_t dim_cord_vect,
+	  const int64_t type_nucl_num,
+	  const int64_t* type_nucl_vector,
+	  const double* cord_vector,
+	  double* const cord_vect_full );
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_cord_vect_full_doc (
+	  const qmckl_context context,
+	  const int64_t nucl_num,
+	  const int64_t dim_cord_vect,
+	  const int64_t type_nucl_num,
+	  const int64_t* type_nucl_vector,
+	  const double* cord_vector,
+	  double* const cord_vect_full );
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_cord_vect_full_hpc (
+	  const qmckl_context context,
+	  const int64_t nucl_num,
+	  const int64_t dim_cord_vect,
+	  const int64_t type_nucl_num,
+	  const int64_t* type_nucl_vector,
+	  const double* cord_vector,
+	  double* const cord_vect_full );
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_cord_vect_full (
+	  const qmckl_context context,
+	  const int64_t nucl_num,
+	  const int64_t dim_cord_vect,
+	  const int64_t type_nucl_num,
+	  const int64_t* type_nucl_vector,
+	  const double* cord_vector,
+	  double* const cord_vect_full ) {
+
+    #ifdef HAVE_HPC
+      return qmckl_compute_cord_vect_full_hpc(context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full);
+    #else
+      return qmckl_compute_cord_vect_full_doc(context, nucl_num, dim_cord_vect, type_nucl_num, type_nucl_vector, cord_vector, cord_vect_full);
+    #endif
+    }
+    #+end_src
+
+
+
 *** Compute lkpm_combined_index
    :PROPERTIES:
    :Name:     qmckl_compute_lkpm_combined_index
@@ -6470,6 +6721,7 @@ assert(fabs(tmp_c[0][0][1][0][0] - 2.7083473948352403) < 1e-12);
 
 printf("%e\n%e\n", tmp_c[0][1][0][0][0],0.237440520852232); 
 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}\)
@@ -6968,6 +7220,7 @@ double factor_een[walk_num];
 rc = qmckl_get_jastrow_factor_een(context, &(factor_een[0]),walk_num);
 
 assert(fabs(factor_een[0] + 0.37407972141304213) < 1e-12);
+return QMCKL_SUCCESS;
      #+end_src
 
 ** Electron-electron-nucleus Jastrow \(f_{een}\) derivative