Moved C version of Jastrow into HPC

org/qmckl_jastrow.org

@@ -1727,7 +1727,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num,
   double precision      , intent(out) :: factor_ee(walk_num)
 
   integer*8 :: i, j, p, ipar, nw
-  double precision   :: pow_ser, x, spin_fact, power_ser
+  double precision   :: x, power_ser, spin_fact
 
   info = QMCKL_SUCCESS
 
@@ -1766,7 +1766,7 @@ integer function qmckl_compute_factor_ee_f(context, walk_num, elec_num, up_num,
           power_ser = power_ser + bord_vector(p + 1) * x
         end do
 
-        if(j .LE. up_num .OR. i .GT. up_num) then
+        if(j <= up_num .OR. i > up_num) then
           spin_fact = 0.5d0
           ipar = 2
         endif
@@ -1785,7 +1785,7 @@ end function qmckl_compute_factor_ee_f
     #+end_src
 
 #+begin_src c   :comments org :tangle (eval c) :noweb yes
-    qmckl_exit_code qmckl_compute_factor_ee (
+qmckl_exit_code qmckl_compute_factor_ee (
 	  const qmckl_context context,
 	  const int64_t walk_num,
 	  const int64_t elec_num,
@@ -1797,7 +1797,7 @@ end function qmckl_compute_factor_ee_f
 	  double* const factor_ee ) {
 
   int ipar; // can we use a smaller integer?
-  double pow_ser, x, x1, spin_fact, power_ser;
+  double x, x1, spin_fact, power_ser;
 
   if (context == QMCKL_NULL_CONTEXT) { 
      return QMCKL_INVALID_CONTEXT;
@@ -2493,8 +2493,8 @@ integer function qmckl_compute_factor_en_f(context, walk_num, elec_num, nucl_num
   double precision      , intent(in)  :: en_distance_rescaled(elec_num, nucl_num, walk_num)
   double precision      , intent(out) :: factor_en(walk_num)
 
-  integer*8 :: i, a, p, ipar, nw
+  integer*8 :: i, a, p, nw
-  double precision   :: x, spin_fact, power_ser
+  double precision   :: x, power_ser
 
   info = QMCKL_SUCCESS
 
@@ -2563,10 +2563,7 @@ qmckl_exit_code qmckl_compute_factor_en (
 	  const double* en_distance_rescaled,
 	  double* const factor_en ) {
 
-
+  double  x, x1, power_ser;
-  int     ipar;
-  double  x, x1, spin_fact, power_ser;
-
 
   if (context == QMCKL_NULL_CONTEXT) {
      return QMCKL_INVALID_CONTEXT;
@@ -2584,10 +2581,30 @@ qmckl_exit_code qmckl_compute_factor_en (
      return QMCKL_INVALID_ARG_4;
   }
 
+  if (type_nucl_num <= 0) {
+    return QMCKL_INVALID_ARG_5;
+  }
+
+  if (type_nucl_vector == NULL) {
+    return QMCKL_INVALID_ARG_6;
+  }
+
   if (aord_num <= 0) {
      return QMCKL_INVALID_ARG_7;
   }
 
+  if (aord_vector == NULL) {
+    return QMCKL_INVALID_ARG_8;
+  }
+
+  if (en_distance_rescaled == NULL) {
+    return QMCKL_INVALID_ARG_9;
+  }
+
+  if (factor_en == NULL) {
+    return QMCKL_INVALID_ARG_10;
+  }
+
 
   for (int nw = 0; nw < walk_num; ++nw ) {
     // init array
@@ -2826,7 +2843,7 @@ integer function qmckl_compute_factor_en_deriv_e_f(context, walk_num, elec_num,
   double precision      , intent(out) :: factor_en_deriv_e(elec_num,4,walk_num)
 
   integer*8 :: i, a, p, ipar, nw, ii
-  double precision   :: x, spin_fact, den, invden, invden2, invden3, xinv
+  double precision   :: x, den, invden, invden2, invden3, xinv
   double precision   :: y, lap1, lap2, lap3, third
   double precision, dimension(3) :: power_ser_g
   double precision, dimension(4) :: dx
@@ -5264,7 +5281,7 @@ qmckl_exit_code qmckl_compute_lkpm_combined_index (
     | ~tmp_c~          | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out    | vector of non-zero coefficients   |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
+integer function qmckl_compute_tmp_c_doc_f(context, cord_num, elec_num, nucl_num, &
      walk_num, een_rescaled_e, een_rescaled_n, tmp_c) &
      result(info)
   use qmckl
@@ -5319,7 +5336,7 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
 
   do nw=1, walk_num
   do i=0, cord_num-1
-  info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     &
+  info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,     &
                      een_rescaled_e(1,1,i,nw),LDA*1_8,                     &
                      een_rescaled_n(1,1,0,nw),LDB*1_8,                     &
                      beta,                                       &
@@ -5327,11 +5344,39 @@ integer function qmckl_compute_tmp_c_f(context, cord_num, elec_num, nucl_num, &
   end do
   end do
 
-end function qmckl_compute_tmp_c_f
+end function qmckl_compute_tmp_c_doc_f
     #+end_src
 
+#+CALL: generate_c_interface(table=qmckl_factor_tmp_c_args,rettyp=get_value("FRetType"),fname="qmckl_compute_tmp_c_doc")
+
+#+RESULTS:
+#+begin_src f90 :tangle (eval f) :comments org :exports none
+integer(c_int32_t) function qmckl_compute_tmp_c_doc &
+    (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c) &
+    bind(C) result(info)
+
+  use, intrinsic :: iso_c_binding
+  implicit none
+
+  integer (c_int64_t) , intent(in)  , value :: context
+  integer (c_int64_t) , intent(in)  , value :: cord_num
+  integer (c_int64_t) , intent(in)  , value :: elec_num
+  integer (c_int64_t) , intent(in)  , value :: nucl_num
+  integer (c_int64_t) , intent(in)  , value :: walk_num
+  real    (c_double ) , intent(in)          :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num)
+  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
+  real    (c_double ) , intent(out)         :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
+
+  integer(c_int32_t), external :: qmckl_compute_tmp_c_doc_f
+  info = qmckl_compute_tmp_c_doc_f &
+         (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c)
+
+end function qmckl_compute_tmp_c_doc
+#+end_src
+
+
     #+begin_src c :comments org :tangle (eval c) :noweb yes
-qmckl_exit_code qmckl_compute_tmp_c (
+qmckl_exit_code qmckl_compute_tmp_c_hpc (
 	  const qmckl_context context,
 	  const int64_t cord_num,
 	  const int64_t elec_num,
@@ -5341,17 +5386,6 @@ qmckl_exit_code qmckl_compute_tmp_c (
 	  const double* een_rescaled_n,
 	  double* const tmp_c ) {
 
-  qmckl_exit_code info;
-  int          i, j, a, l, kk, p, lmax, nw;
-  char         TransA, TransB;
-  double       alpha, beta;
-  int          M, N, K, LDA, LDB, LDC;
-
-  TransA = 'N';
-  TransB = 'N';
-  alpha  = 1.0;
-  beta   = 0.0;
-
   if (context == QMCKL_NULL_CONTEXT) {
     return QMCKL_INVALID_CONTEXT;
   }
@@ -5368,26 +5402,37 @@ qmckl_exit_code qmckl_compute_tmp_c (
     return QMCKL_INVALID_ARG_4;
   } 
 
-  M = elec_num;
+  if (walk_num <= 0) {
-  N = nucl_num*(cord_num + 1);
+    return QMCKL_INVALID_ARG_5;
-  K = elec_num;
+  } 
 
-  LDA = sizeof(een_rescaled_e)/sizeof(double);
+  qmckl_exit_code info = QMCKL_SUCCESS;
-  LDB = sizeof(een_rescaled_n)/sizeof(double);
-  LDC = sizeof(tmp_c)/sizeof(double);
 
-  for (int nw=0; nw < walk_num; ++nw) {
+  const char  TransA = 'N';
-    for (int i=0; i<cord_num; ++i){
+  const char  TransB = 'N';
-    info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     \
+  const double alpha = 1.0;
+  const double beta  = 0.0;
+
+  const int64_t M = elec_num;
+  const int64_t N = nucl_num*(cord_num + 1);
+  const int64_t K = elec_num;
+
+  const int64_t LDA = elec_num;
+  const int64_t LDB = elec_num;
+  const int64_t LDC = elec_num;
+
+  for (int64_t nw=0; nw < walk_num; ++nw) {
+    for (int64_t i=0; i<cord_num; ++i){
+    info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,     \
                  //    &een_rescaled_e[0+0*elec_num+i*elec_num*elec_num+nw*elec_num*elec_num*(cord_num+1)], 
-                       &een_rescaled_e[             i*elec_num*elec_num+nw*elec_num*elec_num*(cord_num+1)], \
+                       &(een_rescaled_e[             i*elec_num*elec_num+nw*elec_num*elec_num*(cord_num+1)]), \
                        LDA,                                        \
                  //    &een_rescaled_n[0+0*elec_num+0*elec_num*nucl_num+nw*elec_num*nucl_num*(cord_num+1)], 
-                       &een_rescaled_n[                                 nw*elec_num*nucl_num*(cord_num+1)], \
+                       &(een_rescaled_n[                                 nw*elec_num*nucl_num*(cord_num+1)]), \
                        LDB,                                        \
                        beta,                                       \
                  //    &tmp_c[0+0*elec_num+0*elec_num*nucl_num+i*elec_num*nucl_num*(cord_num+1)+nw*elec_num*nucl_num*(cord_num+1)*cord_num], 
-                       &tmp_c[                                 i*elec_num*nucl_num*(cord_num+1)+nw*elec_num*nucl_num*(cord_num+1)*cord_num],        \
+                       &(tmp_c[                                 i*elec_num*nucl_num*(cord_num+1)+nw*elec_num*nucl_num*(cord_num+1)*cord_num]), \
                        LDC);
     }
   }
@@ -5396,19 +5441,67 @@ qmckl_exit_code qmckl_compute_tmp_c (
 }
     #+end_src
 
-    #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+    #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname="qmckl_compute_tmp_c")
 
     #+RESULTS:
     #+begin_src c :tangle (eval h_func) :comments org
-    qmckl_exit_code qmckl_compute_tmp_c (
+qmckl_exit_code qmckl_compute_tmp_c (
-	  const qmckl_context context,
+          const qmckl_context context,
-	  const int64_t cord_num,
+          const int64_t cord_num,
-	  const int64_t elec_num,
+          const int64_t elec_num,
-	  const int64_t nucl_num,
+          const int64_t nucl_num,
-	  const int64_t walk_num,
+          const int64_t walk_num,
-	  const double* een_rescaled_e,
+          const double* een_rescaled_e,
-	  const double* een_rescaled_n,
+          const double* een_rescaled_n,
-	  double* const tmp_c );
+          double* const tmp_c ); 
+    #+end_src
+
+#   #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname="qmckl_compute_tmp_c_doc")
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_private_func) :comments org
+qmckl_exit_code qmckl_compute_tmp_c_doc (
+          const qmckl_context context,
+          const int64_t cord_num,
+          const int64_t elec_num,
+          const int64_t nucl_num,
+          const int64_t walk_num,
+          const double* een_rescaled_e,
+          const double* een_rescaled_n,
+          double* const tmp_c ); 
+    #+end_src
+
+#   #+CALL: generate_c_header(table=qmckl_factor_tmp_c_args,rettyp=get_value("CRetType"),fname="qmckl_compute_tmp_c_hpc")
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_private_func) :comments org
+qmckl_exit_code qmckl_compute_tmp_c_hpc (const qmckl_context context,
+                                         const int64_t cord_num,
+                                         const int64_t elec_num,
+                                         const int64_t nucl_num,
+                                         const int64_t walk_num,
+                                         const double* een_rescaled_e,
+                                         const double* een_rescaled_n,
+                                         double* const tmp_c ); 
+    #+end_src
+
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_tmp_c (const qmckl_context context,
+                                     const int64_t cord_num,
+                                     const int64_t elec_num,
+                                     const int64_t nucl_num,
+                                     const int64_t walk_num,
+                                     const double* een_rescaled_e,
+                                     const double* een_rescaled_n,
+                                     double* const tmp_c )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_tmp_c_hpc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c);
+#else
+  return qmckl_compute_tmp_c_doc(context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e, een_rescaled_n, tmp_c);
+#endif
+}
     #+end_src
   
  TODO: FIX dtmp_c dimension in the table.
@@ -5432,7 +5525,7 @@ qmckl_exit_code qmckl_compute_tmp_c (
     | ~dtmp_c~                 | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | out    | vector of non-zero coefficients               |
 
     #+begin_src f90 :comments org :tangle (eval f) :noweb yes
-integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, &
+integer function qmckl_compute_dtmp_c_doc_f(context, cord_num, elec_num, nucl_num, &
      walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c) &
      result(info)
   use qmckl
@@ -5486,38 +5579,56 @@ integer function qmckl_compute_dtmp_c_f(context, cord_num, elec_num, nucl_num, &
   LDC = 4*size(dtmp_c,1)
 
   do nw=1, walk_num
-  do i=0, cord_num-1
+     do i=0, cord_num-1
-  info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,     &
+        info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,  &
-                     een_rescaled_e_deriv_e(1,1,1,i,nw),LDA*1_8,                     &
+             een_rescaled_e_deriv_e(1,1,1,i,nw),LDA*1_8,            &
-                     een_rescaled_n(1,1,0,nw),LDB*1_8,                     &
+             een_rescaled_n(1,1,0,nw),LDB*1_8,                      &
-                     beta,                                       &
+             beta,                                                  &
-                     dtmp_c(1,1,1,0,i,nw),LDC)
+             dtmp_c(1,1,1,0,i,nw),LDC)
-  end do
+     end do
   end do
   
-end function qmckl_compute_dtmp_c_f
+end function qmckl_compute_dtmp_c_doc_f
     #+end_src
 
+#+CALL: generate_c_interface(table=qmckl_factor_dtmp_c_args,rettyp=get_value("FRetType"),fname="qmckl_compute_dtmp_c_doc")
+
+#+RESULTS:
+#+begin_src f90 :tangle (eval f) :comments org :exports none
+integer(c_int32_t) function qmckl_compute_dtmp_c_doc &
+    (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c) &
+    bind(C) result(info)
+
+  use, intrinsic :: iso_c_binding
+  implicit none
+
+  integer (c_int64_t) , intent(in)  , value :: context
+  integer (c_int64_t) , intent(in)  , value :: cord_num
+  integer (c_int64_t) , intent(in)  , value :: elec_num
+  integer (c_int64_t) , intent(in)  , value :: nucl_num
+  integer (c_int64_t) , intent(in)  , value :: walk_num
+  real    (c_double ) , intent(in)          :: een_rescaled_e_deriv_e(elec_num,4,elec_num,0:cord_num,walk_num)
+  real    (c_double ) , intent(in)          :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num)
+  real    (c_double ) , intent(out)         :: dtmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num)
+
+  integer(c_int32_t), external :: qmckl_compute_dtmp_c_doc_f
+  info = qmckl_compute_dtmp_c_doc_f &
+         (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e, een_rescaled_n, dtmp_c)
+
+end function qmckl_compute_dtmp_c_doc
+#+end_src
+
+    
     #+begin_src c :comments org :tangle (eval c) :noweb yes
-qmckl_exit_code qmckl_compute_dtmp_c (
+qmckl_exit_code qmckl_compute_dtmp_c_hpc (const qmckl_context context,
-	  const qmckl_context context,
+                                          const int64_t cord_num,
-	  const int64_t cord_num,
+                                          const int64_t elec_num,
-	  const int64_t elec_num,
+                                          const int64_t nucl_num,
-	  const int64_t nucl_num,
+                                          const int64_t walk_num,
-	  const int64_t walk_num,
+                                          const double* een_rescaled_e_deriv_e,
-	  const double* een_rescaled_e_deriv_e,
+                                          const double* een_rescaled_n,
-	  const double* een_rescaled_n,
+                                          double* const dtmp_c )
-	  double* const dtmp_c ) {
+{
-
-  qmckl_exit_code  info;
-  char             TransA, TransB;
-  double           alpha, beta;
-  int              M, N, K, LDA, LDB, LDC;
-
-  TransA = 'N';
-  TransB = 'N';
-  alpha = 1.0;
-  beta  = 0.0;
 
   if (context == QMCKL_NULL_CONTEXT) {
      return QMCKL_INVALID_CONTEXT;
@@ -5535,26 +5646,37 @@ qmckl_exit_code qmckl_compute_dtmp_c (
      return QMCKL_INVALID_ARG_4;
   } 
 
-  M = 4*elec_num;
+  if (walk_num <= 0) {
-  N = nucl_num*(cord_num + 1);
+     return QMCKL_INVALID_ARG_5;
-  K = elec_num;
+  } 
 
-  LDA = 4*sizeof(een_rescaled_e_deriv_e)/sizeof(double);
+  qmckl_exit_code  info = QMCKL_SUCCESS;
-  LDB = sizeof(een_rescaled_n)/sizeof(double);
-  LDC = 4*sizeof(dtmp_c)/sizeof(double);
 
-  for (int nw=0; nw < walk_num; ++nw) {
+  const char  TransA = 'N';
-    for (int i=0; nw < cord_num; ++i) {
+  const char  TransB = 'N';
-      info = qmckl_dgemm(context,TransA, TransB, M, N, K, alpha,        \
+  const double alpha = 1.0;
+  const double beta  = 0.0;
+
+  const int64_t M = 4*elec_num;
+  const int64_t N = nucl_num*(cord_num + 1);
+  const int64_t K = elec_num;
+
+  const int64_t LDA = 4*elec_num;
+  const int64_t LDB = elec_num;
+  const int64_t LDC = 4*elec_num;
+
+  for (int64_t nw=0; nw < walk_num; ++nw) {
+    for (int64_t i=0; nw < cord_num; ++i) {
+      info = qmckl_dgemm(context, TransA, TransB, M, N, K, alpha,        \
               //&een_rescaled_e_deriv_e[0+0*elec_num+0*elec_num*4+i*elec_num*4*elec_num+nw*elec_num*4*elec_num*(cord_num+1)], 
-                &een_rescaled_e_deriv_e[i*elec_num*4*elec_num+nw*elec_num*4*elec_num*(cord_num+1)],                \
+                &(een_rescaled_e_deriv_e[i*elec_num*4*elec_num+nw*elec_num*4*elec_num*(cord_num+1)]), \
                 LDA,                                               \
               //&een_rescaled_n[0+0*elec_num+0*elec_num*nucl_num+nw*elec_num*nucl_num*(cord_num+1)],  
-                &een_rescaled_n[nw*elec_num*nucl_num*(cord_num+1)],   \
+                &(een_rescaled_n[nw*elec_num*nucl_num*(cord_num+1)]),   \
                 LDB,                                               \
                 beta,                                              \
               //&dtmp_c[0+0*elec_num+0*elec_num*4+0*elec_num*4*nucl_num+i*elec_num*4*nucl_num*(cord_num+1)+nw*elec_num*4*nucl_num*(cord_num+1)*cord_num], 
-                &dtmp_c[i*elec_num*4*nucl_num*(cord_num+1)+nw*elec_num*4*nucl_num*(cord_num+1)*cord_num], \
+                &(dtmp_c[i*elec_num*4*nucl_num*(cord_num+1)+nw*elec_num*4*nucl_num*(cord_num+1)*cord_num]), \
                 LDC);
     }
   }
@@ -5563,11 +5685,41 @@ qmckl_exit_code qmckl_compute_dtmp_c (
 }
     #+end_src
 
-    #+CALL: generate_c_header(table=qmckl_factor_dtmp_c_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
 
-    #+RESULTS:
+    #+begin_src c :comments org :tangle (eval c) :noweb yes
+qmckl_exit_code qmckl_compute_dtmp_c (const qmckl_context context,
+                                      const int64_t cord_num,
+                                      const int64_t elec_num,
+                                      const int64_t nucl_num,
+                                      const int64_t walk_num,
+                                      const double* een_rescaled_e_deriv_e,
+                                      const double* een_rescaled_n,
+                                      double* const dtmp_c )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_dtmp_c_hpc (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e,
+                                   een_rescaled_n, dtmp_c );
+#else
+  return qmckl_compute_dtmp_c_doc (context, cord_num, elec_num, nucl_num, walk_num, een_rescaled_e_deriv_e,
+                                   een_rescaled_n, dtmp_c );
+#endif
+}
+    #+end_src
+
+
     #+begin_src c :tangle (eval h_func) :comments org
-    qmckl_exit_code qmckl_compute_dtmp_c (
+qmckl_exit_code qmckl_compute_dtmp_c (const qmckl_context context,
+                                      const int64_t cord_num,
+                                      const int64_t elec_num,
+                                      const int64_t nucl_num,
+                                      const int64_t walk_num,
+                                      const double* een_rescaled_e_deriv_e,
+                                      const double* een_rescaled_n,
+                                      double* const dtmp_c );
+    #+end_src
+
+    #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_dtmp_c_doc (
 	  const qmckl_context context,
 	  const int64_t cord_num,
 	  const int64_t elec_num,
@@ -5578,6 +5730,17 @@ qmckl_exit_code qmckl_compute_dtmp_c (
 	  double* const dtmp_c );
     #+end_src
 
+    #+begin_src c :tangle (eval h_private_func) :comments org
+    qmckl_exit_code qmckl_compute_dtmp_c_hpc (
+	  const qmckl_context context,
+	  const int64_t cord_num,
+	  const int64_t elec_num,
+	  const int64_t nucl_num,
+	  const int64_t walk_num,
+	  const double* een_rescaled_e_deriv_e,
+	  const double* een_rescaled_n,
+	  double* const dtmp_c );
+    #+end_src
 
 *** Test