diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index 9b0370c..d4eaaa8 100644 --- a/org/qmckl_jastrow.org +++ b/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 - double precision :: x, spin_fact, power_ser + integer*8 :: i, a, p, nw + 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 ) { - - int ipar; - double x, x1, spin_fact, power_ser; - + double x, x1, 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; - N = nucl_num*(cord_num + 1); - K = elec_num; - - LDA = sizeof(een_rescaled_e)/sizeof(double); - LDB = sizeof(een_rescaled_n)/sizeof(double); - LDC = sizeof(tmp_c)/sizeof(double); + if (walk_num <= 0) { + return QMCKL_INVALID_ARG_5; + } - for (int nw=0; nw < walk_num; ++nw) { - for (int i=0; i