>
}
qmckl_exit_code qmckl_set_mo_basis_coefficient(qmckl_context context, const double* coefficient) {
- <>
-
+
int32_t mask = 1 << 1;
+ <>
+
if (ctx->mo_basis.coefficient != NULL) {
qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient);
if (rc != QMCKL_SUCCESS) {
@@ -383,7 +395,7 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
NULL);
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
@@ -421,7 +433,464 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
* Computation
-** Computation of MOs
+** Computation of MOs: values only
+
+*** Get
+
+ #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value(qmckl_context context,
+ double* const mo_value,
+ const int64_t size_max);
+ #+end_src
+
+ #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value(qmckl_context context,
+ double* const mo_value,
+ const int64_t size_max)
+{
+
+ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+ return QMCKL_NULL_CONTEXT;
+ }
+
+ qmckl_exit_code rc;
+
+ rc = qmckl_provide_ao_basis_ao_value(context);
+ if (rc != QMCKL_SUCCESS) return rc;
+
+ rc = qmckl_provide_mo_value(context);
+ if (rc != QMCKL_SUCCESS) return rc;
+
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+ assert (ctx != NULL);
+
+ const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num;
+ if (size_max < sze) {
+ return qmckl_failwith( context,
+ QMCKL_INVALID_ARG_3,
+ "qmckl_get_mo_basis_mo_value",
+ "input array too small");
+ }
+ memcpy(mo_value, ctx->mo_basis.mo_value, sze * sizeof(double));
+
+ return QMCKL_SUCCESS;
+}
+ #+end_src
+
+ #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+ interface
+ integer(c_int32_t) function qmckl_get_mo_basis_mo_value (context, &
+ mo_value, size_max) bind(C)
+ use, intrinsic :: iso_c_binding
+ import
+ implicit none
+
+ integer (c_int64_t) , intent(in) , value :: context
+ double precision, intent(out) :: mo_value(*)
+ integer (c_int64_t) , intent(in) , value :: size_max
+ end function qmckl_get_mo_basis_mo_value
+ end interface
+ #+end_src
+
+ Uses the given array to compute the values.
+
+ #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
+ double* const mo_value,
+ const int64_t size_max);
+ #+end_src
+
+ #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
+ double* const mo_value,
+ const int64_t size_max)
+{
+
+ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+ return qmckl_failwith( context,
+ QMCKL_INVALID_CONTEXT,
+ "qmckl_get_mo_basis_mo_value",
+ NULL);
+ }
+
+ qmckl_exit_code rc;
+
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+ assert (ctx != NULL);
+
+ const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num;
+ if (size_max < sze) {
+ return qmckl_failwith( context,
+ QMCKL_INVALID_ARG_3,
+ "qmckl_get_mo_basis_mo_value",
+ "input array too small");
+ }
+
+ rc = qmckl_context_touch(context);
+ if (rc != QMCKL_SUCCESS) return rc;
+
+ double* old_array = ctx->mo_basis.mo_value;
+
+ ctx->mo_basis.mo_value = mo_value;
+
+ rc = qmckl_provide_mo_value(context);
+ if (rc != QMCKL_SUCCESS) return rc;
+
+ ctx->mo_basis.mo_value = old_array;
+
+ return QMCKL_SUCCESS;
+}
+ #+end_src
+
+ #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+ interface
+ integer(c_int32_t) function qmckl_get_mo_basis_mo_value_inplace (context, &
+ mo_value, size_max) bind(C)
+ use, intrinsic :: iso_c_binding
+ import
+ implicit none
+ integer (c_int64_t) , intent(in) , value :: context
+ double precision, intent(out) :: mo_value(*)
+ integer (c_int64_t) , intent(in) , value :: size_max
+ end function qmckl_get_mo_basis_mo_value_inplace
+ end interface
+ #+end_src
+
+*** Provide
+
+ #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_mo_value(qmckl_context context);
+ #+end_src
+
+ #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
+{
+
+ qmckl_exit_code rc;
+ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+ return QMCKL_NULL_CONTEXT;
+ }
+
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+ assert (ctx != NULL);
+
+ if (!ctx->ao_basis.provided) {
+ return qmckl_failwith( context,
+ QMCKL_NOT_PROVIDED,
+ "qmckl_ao_basis",
+ NULL);
+ }
+
+ rc = qmckl_provide_ao_basis_ao_value(context);
+ if (rc != QMCKL_SUCCESS) {
+ return qmckl_failwith( context,
+ QMCKL_NOT_PROVIDED,
+ "qmckl_ao_value",
+ NULL);
+ }
+
+ if (!ctx->mo_basis.provided) {
+ return qmckl_failwith( context,
+ QMCKL_NOT_PROVIDED,
+ "qmckl_mo_basis",
+ NULL);
+ }
+
+ /* Compute if necessary */
+ if (ctx->point.date > ctx->mo_basis.mo_value_date) {
+
+ /* Allocate array */
+ if (ctx->mo_basis.mo_value == NULL) {
+
+ qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+ mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
+ double* mo_value = (double*) qmckl_malloc(context, mem_info);
+
+ if (mo_value == NULL) {
+ return qmckl_failwith( context,
+ QMCKL_ALLOCATION_FAILED,
+ "qmckl_mo_basis_mo_value",
+ NULL);
+ }
+ ctx->mo_basis.mo_value = mo_value;
+ }
+
+ if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {
+
+ // mo_vgl has been computed at this step: Just copy the data.
+
+ double * v = &(ctx->mo_basis.mo_value[0]);
+ double * vgl = &(ctx->mo_basis.mo_vgl[0]);
+ for (int i=0 ; ipoint.num ; ++i) {
+ for (int k=0 ; kmo_basis.mo_num ; ++k) {
+ v[k] = vgl[k];
+ }
+ v += ctx->mo_basis.mo_num;
+ vgl += ctx->mo_basis.mo_num * 5;
+ }
+
+ } else {
+
+ rc = qmckl_compute_mo_basis_mo_value(context,
+ ctx->ao_basis.ao_num,
+ ctx->mo_basis.mo_num,
+ ctx->point.num,
+ ctx->mo_basis.coefficient_t,
+ ctx->ao_basis.ao_value,
+ ctx->mo_basis.mo_value);
+
+ if (rc != QMCKL_SUCCESS) {
+ return rc;
+ }
+
+ }
+
+ ctx->mo_basis.mo_value_date = ctx->date;
+ }
+
+ return QMCKL_SUCCESS;
+}
+ #+end_src
+
+*** Compute
+ :PROPERTIES:
+ :Name: qmckl_compute_mo_basis_mo_value
+ :CRetType: qmckl_exit_code
+ :FRetType: qmckl_exit_code
+ :END:
+
+ #+NAME: qmckl_mo_basis_mo_value_args
+ | Variable | Type | In/Out | Description |
+ |-----------------+-----------------------------+--------+-------------------------------------------------|
+ | ~context~ | ~qmckl_context~ | in | Global state |
+ | ~ao_num~ | ~int64_t~ | in | Number of AOs |
+ | ~mo_num~ | ~int64_t~ | in | Number of MOs |
+ | ~point_num~ | ~int64_t~ | in | Number of points |
+ | ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
+ | ~ao_value~ | ~double[point_num][ao_num]~ | in | Value of the AOs |
+ | ~mo_value~ | ~double[point_num][mo_num]~ | out | Value of the MOs |
+
+
+ The matrix of AO values is very sparse, so we use a sparse-dense
+ matrix multiplication instead of a dgemm, as exposed in
+ https://dx.doi.org/10.1007/978-3-642-38718-0_14.
+
+
+
+ #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
+ ao_num, mo_num, point_num, &
+ coefficient_t, ao_value, mo_value) &
+ result(info)
+ use qmckl
+ implicit none
+ integer(qmckl_context), intent(in) :: context
+ integer*8 , intent(in) :: ao_num, mo_num
+ integer*8 , intent(in) :: point_num
+ double precision , intent(in) :: ao_value(ao_num,point_num)
+ double precision , intent(in) :: coefficient_t(mo_num,ao_num)
+ double precision , intent(out) :: mo_value(mo_num,point_num)
+ integer*8 :: i,j,k
+ double precision :: c1, c2, c3, c4, c5
+
+ integer*8 :: LDA, LDB, LDC
+
+ info = QMCKL_SUCCESS
+ if (.True.) then ! fast algorithm
+ do j=1,point_num
+ mo_value(:,j) = 0.d0
+ do k=1,ao_num
+ if (ao_value(k,j) /= 0.d0) then
+ c1 = ao_value(k,j)
+ do i=1,mo_num
+ mo_value(i,j) = mo_value(i,j) + coefficient_t(i,k) * c1
+ end do
+ end if
+ end do
+ end do
+
+ else ! dgemm
+
+ LDA = size(coefficient_t,1)
+ LDB = size(ao_value,1)
+ LDC = size(mo_value,1)
+
+ info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, &
+ coefficient_t, LDA, ao_value, LDB, &
+ 0.d0, mo_value, LDC)
+
+ end if
+
+end function qmckl_compute_mo_basis_mo_value_doc_f
+ #+end_src
+
+ #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value"))
+
+ #+RESULTS:
+ #+begin_src c :tangle (eval h_func) :comments org
+ qmckl_exit_code qmckl_compute_mo_basis_mo_value (
+ const qmckl_context context,
+ const int64_t ao_num,
+ const int64_t mo_num,
+ const int64_t point_num,
+ const double* coefficient_t,
+ const double* ao_value,
+ double* const mo_value );
+ #+end_src
+
+ #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
+
+ #+RESULTS:
+ #+begin_src c :tangle (eval h_func) :comments org
+ qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
+ const qmckl_context context,
+ const int64_t ao_num,
+ const int64_t mo_num,
+ const int64_t point_num,
+ const double* coefficient_t,
+ const double* ao_value,
+ double* const mo_value );
+ #+end_src
+
+ #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
+
+ #+RESULTS:
+ #+begin_src f90 :tangle (eval f) :comments org :exports none
+ integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc &
+ (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value) &
+ 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 :: ao_num
+ integer (c_int64_t) , intent(in) , value :: mo_num
+ integer (c_int64_t) , intent(in) , value :: point_num
+ real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
+ real (c_double ) , intent(in) :: ao_value(ao_num,point_num)
+ real (c_double ) , intent(out) :: mo_value(mo_num,point_num)
+
+ integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f
+ info = qmckl_compute_mo_basis_mo_value_doc_f &
+ (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value)
+
+ end function qmckl_compute_mo_basis_mo_value_doc
+ #+end_src
+
+ #+begin_src c :tangle (eval c) :comments org
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value (const qmckl_context context,
+ const int64_t ao_num,
+ const int64_t mo_num,
+ const int64_t point_num,
+ const double* coefficient_t,
+ const double* ao_value,
+ double* const mo_value )
+{
+#ifdef HAVE_HPC
+ return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
+#else
+ return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
+#endif
+}
+ #+end_src
+
+*** HPC version
+
+
+ #+begin_src c :tangle (eval h_func) :comments org
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
+ const int64_t ao_num,
+ const int64_t mo_num,
+ const int64_t point_num,
+ const double* coefficient_t,
+ const double* ao_value,
+ double* const mo_value );
+#endif
+ #+end_src
+
+ #+begin_src c :tangle (eval c) :comments org
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
+ const int64_t ao_num,
+ const int64_t mo_num,
+ const int64_t point_num,
+ const double* restrict coefficient_t,
+ const double* restrict ao_value,
+ double* restrict const mo_value )
+{
+ assert (context != QMCKL_NULL_CONTEXT);
+
+#ifdef HAVE_OPENMP
+ #pragma omp parallel for
+#endif
+ for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
+ double* restrict const vgl1 = &(mo_value[ipoint*mo_num]);
+ const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);
+
+ for (int64_t i=0 ; ipoint.num * 5 * ctx->mo_basis.mo_num;
@@ -507,7 +976,7 @@ qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
qmckl_exit_code rc;
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->mo_basis.mo_num * 5 * ctx->point.num;
@@ -563,7 +1032,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
return QMCKL_NULL_CONTEXT;
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->ao_basis.provided) {
@@ -639,7 +1108,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
| ~ao_num~ | ~int64_t~ | in | Number of AOs |
| ~mo_num~ | ~int64_t~ | in | Number of MOs |
| ~point_num~ | ~int64_t~ | in | Number of points |
- | ~coef_normalized_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
+ | ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix |
| ~ao_vgl~ | ~double[point_num][5][ao_num]~ | in | Value, gradients and Laplacian of the AOs |
| ~mo_vgl~ | ~double[point_num][5][mo_num]~ | out | Value, gradients and Laplacian of the MOs |
@@ -653,7 +1122,7 @@ qmckl_exit_code qmckl_provide_mo_vgl(qmckl_context context)
#+begin_src f90 :comments org :tangle (eval f) :noweb yes
integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
ao_num, mo_num, point_num, &
- coef_normalized_t, ao_vgl, mo_vgl) &
+ coefficient_t, ao_vgl, mo_vgl) &
result(info)
use qmckl
implicit none
@@ -661,7 +1130,7 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
integer*8 , intent(in) :: ao_num, mo_num
integer*8 , intent(in) :: point_num
double precision , intent(in) :: ao_vgl(ao_num,5,point_num)
- double precision , intent(in) :: coef_normalized_t(mo_num,ao_num)
+ double precision , intent(in) :: coefficient_t(mo_num,ao_num)
double precision , intent(out) :: mo_vgl(mo_num,5,point_num)
integer*8 :: i,j,k
double precision :: c1, c2, c3, c4, c5
@@ -676,15 +1145,21 @@ integer function qmckl_compute_mo_basis_mo_vgl_doc_f(context, &
c4 = ao_vgl(k,4,j)
c5 = ao_vgl(k,5,j)
do i=1,mo_num
- mo_vgl(i,1,j) = mo_vgl(i,1,j) + coef_normalized_t(i,k) * c1
- mo_vgl(i,2,j) = mo_vgl(i,2,j) + coef_normalized_t(i,k) * c2
- mo_vgl(i,3,j) = mo_vgl(i,3,j) + coef_normalized_t(i,k) * c3
- mo_vgl(i,4,j) = mo_vgl(i,4,j) + coef_normalized_t(i,k) * c4
- mo_vgl(i,5,j) = mo_vgl(i,5,j) + coef_normalized_t(i,k) * c5
+ mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1
+ mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2
+ mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3
+ mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
+ mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
end do
end if
end do
end do
+ info = QMCKL_SUCCESS
+
+! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, &
+! coefficient_t, int(size(coefficient_t,1),8), &
+! ao_vgl, int(size(ao_vgl,1),8), 0.d0, &
+! mo_vgl, int(size(mo_vgl,1),8))
end function qmckl_compute_mo_basis_mo_vgl_doc_f
#+end_src
@@ -698,7 +1173,7 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
- const double* coef_normalized_t,
+ const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#+end_src
@@ -712,7 +1187,7 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
- const double* coef_normalized_t,
+ const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#+end_src
@@ -722,7 +1197,7 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_compute_mo_basis_mo_vgl_doc &
- (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl) &
+ (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
@@ -732,13 +1207,13 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f
integer (c_int64_t) , intent(in) , value :: ao_num
integer (c_int64_t) , intent(in) , value :: mo_num
integer (c_int64_t) , intent(in) , value :: point_num
- real (c_double ) , intent(in) :: coef_normalized_t(ao_num,mo_num)
+ real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num)
real (c_double ) , intent(in) :: ao_vgl(ao_num,5,point_num)
real (c_double ) , intent(out) :: mo_vgl(mo_num,5,point_num)
integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_vgl_doc_f
info = qmckl_compute_mo_basis_mo_vgl_doc_f &
- (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl)
+ (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl)
end function qmckl_compute_mo_basis_mo_vgl_doc
#+end_src
@@ -749,19 +1224,18 @@ qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
- const double* coef_normalized_t,
+ const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl )
{
#ifdef HAVE_HPC
- return qmckl_compute_mo_basis_mo_vgl_hpc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl);
+ return qmckl_compute_mo_basis_mo_vgl_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl);
#else
- return qmckl_compute_mo_basis_mo_vgl_doc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl);
+ return qmckl_compute_mo_basis_mo_vgl_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl);
#endif
}
#+end_src
-
*** HPC version
@@ -772,7 +1246,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
- const double* coef_normalized_t,
+ const double* coefficient_t,
const double* ao_vgl,
double* const mo_vgl );
#endif
@@ -785,20 +1259,23 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
const int64_t ao_num,
const int64_t mo_num,
const int64_t point_num,
- const double* restrict coef_normalized_t,
+ const double* restrict coefficient_t,
const double* restrict ao_vgl,
double* restrict const mo_vgl )
{
+ assert (context != QMCKL_NULL_CONTEXT);
+
#ifdef HAVE_OPENMP
#pragma omp parallel for
#endif
for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]);
- const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
double* restrict const vgl2 = vgl1 + mo_num;
double* restrict const vgl3 = vgl1 + (mo_num << 1);
double* restrict const vgl4 = vgl1 + (mo_num << 1) + mo_num;
double* restrict const vgl5 = vgl1 + (mo_num << 2);
+
+ const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
const double* restrict avgl2 = avgl1 + ao_num;
const double* restrict avgl3 = avgl1 + (ao_num << 1);
const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
@@ -820,7 +1297,6 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
double av4[ao_num];
double av5[ao_num];
for (int64_t k=0 ; knucleus.uninitialized = (1 << 3) - 1;
@@ -167,7 +167,7 @@ qmckl_get_nucleus_num (const qmckl_context context, int64_t* const num) {
"num is a null pointer");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 0;
@@ -226,7 +226,7 @@ qmckl_get_nucleus_charge (const qmckl_context context,
"charge is a null pointer");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 1;
@@ -293,7 +293,7 @@ qmckl_get_nucleus_rescale_factor (const qmckl_context context,
"rescale_factor_kappa is a null pointer");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->nucleus.rescale_factor_kappa > 0.0);
@@ -351,7 +351,7 @@ qmckl_get_nucleus_coord (const qmckl_context context,
"coord is a null pointer");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int32_t mask = 1 << 2;
@@ -410,7 +410,7 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
return false;
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
return ctx->nucleus.provided;
@@ -422,10 +422,20 @@ bool qmckl_nucleus_provided(const qmckl_context context) {
#+NAME:pre2
#+begin_src c :exports none
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
- return QMCKL_NULL_CONTEXT;
- }
+ return qmckl_failwith( context,
+ QMCKL_NULL_CONTEXT,
+ "qmckl_set_nucleus_*",
+ NULL);
+}
-qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+
+if (mask != 0 && !(ctx->nucleus.uninitialized & mask)) {
+ return qmckl_failwith( context,
+ QMCKL_ALREADY_SET,
+ "qmckl_set_nucleus_*",
+ NULL);
+}
#+end_src
#+NAME:post2
@@ -452,6 +462,8 @@ qmckl_exit_code
qmckl_set_nucleus_num(qmckl_context context,
const int64_t num)
{
+ int32_t mask = 1 << 0;
+
<>
if (num <= 0) {
@@ -461,8 +473,6 @@ qmckl_set_nucleus_num(qmckl_context context,
"num <= 0");
}
- int32_t mask = 1 << 0;
-
ctx->nucleus.num = num;
<>
@@ -498,6 +508,8 @@ qmckl_set_nucleus_charge(qmckl_context context,
const double* charge,
const int64_t size_max)
{
+ int32_t mask = 1 << 1;
+
<>
if (charge == NULL) {
@@ -510,8 +522,6 @@ qmckl_set_nucleus_charge(qmckl_context context,
int64_t num;
qmckl_exit_code rc;
- int32_t mask = 1 << 1;
-
rc = qmckl_get_nucleus_num(context, &num);
if (rc != QMCKL_SUCCESS) return rc;
@@ -569,12 +579,12 @@ qmckl_set_nucleus_coord(qmckl_context context,
const double* coord,
const int64_t size_max)
{
+ int32_t mask = 1 << 2;
+
<>
qmckl_exit_code rc;
- int32_t mask = 1 << 2;
-
const int64_t nucl_num = (int64_t) ctx->nucleus.num;
if (ctx->nucleus.coord.data != NULL) {
@@ -641,6 +651,8 @@ qmckl_exit_code
qmckl_set_nucleus_rescale_factor(qmckl_context context,
const double rescale_factor_kappa)
{
+ int32_t mask = 0; // Can be updated
+
<>
if (rescale_factor_kappa <= 0.0) {
@@ -672,7 +684,6 @@ end interface
** Test
#+begin_src c :tangle (eval c_test)
-const int64_t nucl_num = chbrclf_nucl_num;
const double* nucl_charge = chbrclf_charge;
const double* nucl_coord = &(chbrclf_nucl_coord[0][0]);
const double nucl_rescale_factor_kappa = 2.0;
@@ -688,13 +699,13 @@ rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_NOT_PROVIDED);
-rc = qmckl_set_nucleus_num (context, nucl_num);
+rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
rc = qmckl_get_nucleus_num (context, &n);
assert(rc == QMCKL_SUCCESS);
-assert(n == nucl_num);
+assert(n == chbrclf_nucl_num);
double k;
rc = qmckl_get_nucleus_rescale_factor (context, &k);
@@ -709,41 +720,41 @@ rc = qmckl_get_nucleus_rescale_factor (context, &k);
assert(rc == QMCKL_SUCCESS);
assert(k == nucl_rescale_factor_kappa);
-double nucl_coord2[3*nucl_num];
+double nucl_coord2[3*chbrclf_nucl_num];
-rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num);
+rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_NOT_PROVIDED);
-rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num);
+rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
assert(!qmckl_nucleus_provided(context));
-rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*nucl_num);
+rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num);
assert(rc == QMCKL_SUCCESS);
for (size_t k=0 ; k<3 ; ++k) {
- for (int64_t i=0 ; inucleus.num * ctx->nucleus.num;
@@ -828,7 +839,7 @@ qmckl_exit_code qmckl_provide_nn_distance(qmckl_context context)
return (char) 0;
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
@@ -940,10 +951,10 @@ qmckl_exit_code qmckl_compute_nn_distance (
assert(qmckl_nucleus_provided(context));
-double distance[nucl_num*nucl_num];
-rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num);
+double distance[chbrclf_nucl_num*chbrclf_nucl_num];
+rc = qmckl_get_nucleus_nn_distance(context, distance, chbrclf_nucl_num*chbrclf_nucl_num);
assert(distance[0] == 0.);
-assert(distance[1] == distance[nucl_num]);
+assert(distance[1] == distance[chbrclf_nucl_num]);
assert(fabs(distance[1]-2.070304721365169) < 1.e-12);
#+end_src
@@ -973,7 +984,7 @@ qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context,
qmckl_exit_code rc = qmckl_provide_nn_distance_rescaled(context);
if (rc != QMCKL_SUCCESS) return rc;
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
const int64_t sze = ctx->nucleus.num * ctx->nucleus.num;
@@ -1019,7 +1030,7 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context)
return (char) 0;
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED;
@@ -1167,7 +1178,7 @@ qmckl_exit_code qmckl_get_nucleus_repulsion(qmckl_context context, double* const
qmckl_exit_code rc = qmckl_provide_nucleus_repulsion(context);
if (rc != QMCKL_SUCCESS) return rc;
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
*energy = ctx->nucleus.repulsion;
@@ -1203,7 +1214,7 @@ qmckl_exit_code qmckl_provide_nucleus_repulsion(qmckl_context context)
return (char) 0;
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
diff --git a/org/qmckl_numprec.org b/org/qmckl_numprec.org
index 7c5f26f..3174464 100644
--- a/org/qmckl_numprec.org
+++ b/org/qmckl_numprec.org
@@ -141,7 +141,7 @@ qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const i
"precision > 53");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
/* This should be always true because the context is valid */
assert (ctx != NULL);
@@ -185,7 +185,7 @@ int qmckl_get_numprec_precision(const qmckl_context context) {
"");
}
- const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
return ctx->numprec.precision;
}
#+end_src
@@ -232,7 +232,7 @@ qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int r
"range > 11");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
/* This should be always true because the context is valid */
assert (ctx != NULL);
@@ -275,7 +275,7 @@ int qmckl_get_numprec_range(const qmckl_context context) {
"");
}
- const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
return ctx->numprec.range;
}
#+end_src
diff --git a/org/qmckl_point.org b/org/qmckl_point.org
index 07b0863..6344975 100644
--- a/org/qmckl_point.org
+++ b/org/qmckl_point.org
@@ -77,11 +77,13 @@ int main() {
The following data stored in the context:
- | Variable | Type | Description |
- |----------+----------------+-------------------------------------------|
- | ~num~ | ~int64_t~ | Total number of points |
- | ~date~ | ~uint64_t~ | Last modification date of the coordinates |
- | ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix3 |
+ | Variable | Type | Description |
+ |--------------+----------------+-------------------------------------------|
+ | ~num~ | ~int64_t~ | Total number of points |
+ | ~alloc_num~ | ~int64_t~ | Numer of allocated number of points |
+ | ~date~ | ~uint64_t~ | Last modification date of the coordinates |
+ | ~alloc_date~ | ~uint64_t~ | Last modification date of the allocation |
+ | ~coord~ | ~qmckl_matrix~ | ~num~ \times 3 matrix |
We consider that the matrix is stored 'transposed' and 'normal'
corresponds to the 3 \times ~num~ matrix.
@@ -91,7 +93,9 @@ int main() {
#+begin_src c :comments org :tangle (eval h_private_type)
typedef struct qmckl_point_struct {
int64_t num;
+ int64_t alloc_num;
uint64_t date;
+ uint64_t alloc_date;
qmckl_matrix coord;
} qmckl_point_struct;
@@ -108,7 +112,7 @@ qmckl_exit_code qmckl_init_point(qmckl_context context) {
return false;
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
memset(&(ctx->point), 0, sizeof(qmckl_point_struct));
@@ -148,7 +152,7 @@ qmckl_get_point_num (const qmckl_context context, int64_t* const num) {
"num is a null pointer");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
assert (ctx->point.num >= (int64_t) 0);
@@ -202,7 +206,7 @@ qmckl_get_point(const qmckl_context context,
"coord is a null pointer");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
int64_t point_num = ctx->point.num;
@@ -263,8 +267,9 @@ end interface
#+begin_src c :comments org :tangle (eval h_func)
qmckl_exit_code qmckl_set_point (qmckl_context context,
const char transp,
+ const int64_t num,
const double* coord,
- const int64_t num);
+ const int64_t size_max);
#+end_src
Copy a sequence of ~num~ points $(x,y,z)$ into the context.
@@ -273,14 +278,22 @@ qmckl_exit_code qmckl_set_point (qmckl_context context,
qmckl_exit_code
qmckl_set_point (qmckl_context context,
const char transp,
+ const int64_t num,
const double* coord,
- const int64_t num)
+ const int64_t size_max)
{
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
+ if (size_max < 3*num) {
+ return qmckl_failwith( context,
+ QMCKL_INVALID_ARG_4,
+ "qmckl_set_point",
+ "Array too small");
+ }
+
if (transp != 'N' && transp != 'T') {
return qmckl_failwith( context,
QMCKL_INVALID_ARG_2,
@@ -295,11 +308,11 @@ qmckl_set_point (qmckl_context context,
"coord is a NULL pointer");
}
- qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
assert (ctx != NULL);
qmckl_exit_code rc;
- if (ctx->point.num < num) {
+ if (num > ctx->point.alloc_num) {
if (ctx->point.coord.data != NULL) {
rc = qmckl_matrix_free(context, &(ctx->point.coord));
@@ -313,7 +326,6 @@ qmckl_set_point (qmckl_context context,
"qmckl_set_point",
NULL);
}
-
};
ctx->point.num = num;
@@ -341,6 +353,11 @@ qmckl_set_point (qmckl_context context,
rc = qmckl_context_touch(context);
assert (rc == QMCKL_SUCCESS);
+ if (num > ctx->point.alloc_num) {
+ ctx->point.alloc_num = num;
+ ctx->point.alloc_date = ctx->point.date;
+ };
+
return QMCKL_SUCCESS;
}
@@ -349,15 +366,16 @@ qmckl_set_point (qmckl_context context,
#+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes
interface
integer(c_int32_t) function qmckl_set_point(context, &
- transp, coord, num) bind(C)
+ transp, num, coord, size_max) bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
character(c_char) , intent(in) , value :: transp
- real (c_double ) , intent(in) :: coord(*)
integer (c_int64_t) , intent(in) , value :: num
+ real (c_double ) , intent(in) :: coord(*)
+ integer (c_int64_t) , intent(in) , value :: size_max
end function
end interface
#+end_src
@@ -380,7 +398,7 @@ double coord3[point_num*3];
rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_NOT_PROVIDED);
-rc = qmckl_set_point (context, 'N', coord, point_num);
+rc = qmckl_set_point (context, 'N', point_num, coord, (point_num*3));
assert(rc == QMCKL_SUCCESS);
int64_t n;
@@ -404,7 +422,7 @@ for (int64_t i=0 ; i=0 ; --i) {
const int k = tmp_array[i];
if (k < 0 || k >= nucleus_num) {
@@ -1086,7 +1079,7 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6
qmckl_exit_code rc;
char file_name_new[size_max+1];
- strncpy(file_name_new, file_name, size_max+1);
+ strncpy(file_name_new, file_name, size_max);
file_name_new[size_max] = '\0';
#ifdef HAVE_TREXIO
diff --git a/org/table_of_contents b/org/table_of_contents
index f7cece8..6cc895a 100644
--- a/org/table_of_contents
+++ b/org/table_of_contents
@@ -18,3 +18,4 @@ qmckl_utils.org
qmckl_trexio.org
qmckl_tests.org
qmckl_verificarlo.org
+qmckl_examples.org
diff --git a/python/src/qmckl.i b/python/src/qmckl.i
index e19e1cc..fe64be1 100644
--- a/python/src/qmckl.i
+++ b/python/src/qmckl.i
@@ -49,6 +49,10 @@ import_array();
/* Include typemaps generated by the process_header.py script */
%include "qmckl_include.i"
+/* Handle properly get_point */
+
+
+
/* exception.i is a generic (language-independent) module */
%include "exception.i"
diff --git a/tools/lib.org b/tools/lib.org
index fd198a9..5e3d78d 100644
--- a/tools/lib.org
+++ b/tools/lib.org
@@ -4,10 +4,10 @@
** Defines the name of the current file
#+NAME: filename
- #+begin_src elisp :tangle no
+ #+begin_src elisp :tangle no
(file-name-nondirectory (substring buffer-file-name 0 -4))
#+end_src
-
+
** Function to get the value of a property.
#+NAME: get_value
#+begin_src elisp :var key="Type"
@@ -15,7 +15,6 @@
(org-entry-get nil key t))
#+end_src
-
** Table of function arguments
#+NAME: test
@@ -32,7 +31,7 @@
| ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ |
| ~C~ | ~double[n][ldc]~ | out | Array containing the $m \times n$ matrix $C$ |
| ~ldc~ | ~int64_t~ | in | Leading dimension of array ~C~ |
-
+
*** Fortran-C type conversions
@@ -124,7 +123,7 @@ for d in parse_table(table):
const = "const "
else:
const = ""
-
+
results += [ f" {const}{c_type} {name}" ]
results=',\n'.join(results)
@@ -146,10 +145,9 @@ return template
const double* B,
const int64_t ldb,
double* const C,
- const int64_t ldc );
+ const int64_t ldc );
#+end_src
-
*** Generates a C interface to the Fortran function
#+NAME: generate_c_interface
@@ -258,4 +256,161 @@ return results
#+END_SRC
-
+
+** Creating provide functions
+
+ #+NAME: write_provider_header
+ #+BEGIN_SRC python :var group="GROUP" :var data="DATA" :results drawer :noweb yes :wrap "src c :comments org :tangle (eval h_private_func) :noweb yes :export none"
+template = "qmckl_exit_code qmckl_provide_{{ group }}_{{ data }}(qmckl_context context);"
+
+msg = template.replace("{{ group }}", group) \
+ .replace("{{ data }}", data)
+return msg
+ #+END_SRC
+
+ #+RESULTS: write_provider_header
+ #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :export none
+ qmckl_exit_code qmckl_provide_GROUP_DATA(qmckl_context context);
+ #+end_src
+
+ #+NAME: write_provider_pre
+ #+BEGIN_SRC python :var group="GROUP" :var data="DATA" :var dimension="DIMENSION" :results drawer :noweb yes :wrap "src c :comments org :tangle (eval c) :noweb yes :export none"
+template = """qmckl_exit_code qmckl_provide_{{ group }}_{{ data }}(qmckl_context context)
+{
+
+ qmckl_exit_code rc;
+
+ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+ return qmckl_failwith( context,
+ QMCKL_INVALID_CONTEXT,
+ "qmckl_provide_{{ group }}_{{ data }}",
+ NULL);
+ }
+
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+ assert (ctx != NULL);
+
+ if (!ctx->{{ group }}.provided) {
+ return qmckl_failwith( context,
+ QMCKL_NOT_PROVIDED,
+ "qmckl_provide_{{ group }}_{{ data }}",
+ NULL);
+ }
+
+ /* Compute if necessary */
+ if (ctx->point.date > ctx->{{ group }}.{{ data }}_date) {
+
+ if (ctx->point.alloc_date > ctx->{{ group }}.{{ data }}_date) {
+ if (ctx->{{ group }}.{{ data }} != NULL) {
+ rc = qmckl_free(context, ctx->{{ group }}.{{ data }});
+ assert (rc == QMCKL_SUCCESS);
+ ctx->{{ group }}.{{ data }} = NULL;
+ }
+ }
+
+ /* Allocate array */
+ if (ctx->{{ group }}.{{ data }} == NULL) {
+
+ qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+ mem_info.size = {{ dimension }} * sizeof(double);
+ double* {{ data }} = (double*) qmckl_malloc(context, mem_info);
+
+ if ({{ data }} == NULL) {
+ return qmckl_failwith( context,
+ QMCKL_ALLOCATION_FAILED,
+ "qmckl_{{ group }}_{{ data }}",
+ NULL);
+ }
+ ctx->{{ group }}.{{ data }} = {{ data }};
+ }
+
+"""
+
+msg = template.replace("{{ group }}", group) \
+ .replace("{{ data }}", data) \
+ .replace("{{ dimension }}", dimension)
+return msg
+ #+END_SRC
+
+ #+RESULTS: write_provider_pre
+ #+begin_src c :comments org :tangle (eval c) :noweb yes :export none
+ qmckl_exit_code qmckl_provide_GROUP_DATA(qmckl_context context)
+ {
+
+ qmckl_exit_code rc;
+
+ if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+ return qmckl_failwith( context,
+ QMCKL_INVALID_CONTEXT,
+ "qmckl_provide_GROUP_DATA",
+ NULL);
+ }
+
+ qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+ assert (ctx != NULL);
+
+ if (!ctx->GROUP.provided) {
+ return qmckl_failwith( context,
+ QMCKL_NOT_PROVIDED,
+ "qmckl_provide_GROUP_DATA",
+ NULL);
+ }
+
+ /* Compute if necessary */
+ if (ctx->point.date > ctx->GROUP.DATA_date) {
+
+ if (ctx->point.alloc_date > ctx->GROUP.DATA_date) {
+ rc = qmckl_free(context, ctx->GROUP.DATA);
+ assert (rc == QMCKL_SUCCESS);
+ ctx->GROUP.DATA = NULL;
+ }
+
+ /* Allocate array */
+ if (ctx->GROUP.DATA == NULL) {
+
+ qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+ mem_info.size = DIMENSION * sizeof(double);
+ double* DATA = (double*) qmckl_malloc(context, mem_info);
+
+ if (DATA == NULL) {
+ return qmckl_failwith( context,
+ QMCKL_ALLOCATION_FAILED,
+ "qmckl_GROUP_DATA",
+ NULL);
+ }
+ ctx->GROUP.DATA = DATA;
+ }
+
+ #+end_src
+
+ #+NAME: write_provider_post
+ #+BEGIN_SRC python :var group="BASIS" :var data="DATA" :results drawer :noweb yes :wrap "src c :comments org :tangle (eval c) :noweb yes :export none"
+template = """ if (rc != QMCKL_SUCCESS) {
+ return rc;
+ }
+
+ ctx->{{ group }}.{{ data }}_date = ctx->date;
+ }
+
+ return QMCKL_SUCCESS;
+}
+"""
+
+msg = template.replace("{{ group }}", group) \
+ .replace("{{ data }}", data)
+
+return msg
+ #+END_SRC
+
+ #+RESULTS: write_provider_post
+ #+begin_src c :comments org :tangle (eval c) :noweb yes :export none
+ if (rc != QMCKL_SUCCESS) {
+ return rc;
+ }
+
+ ctx->BASIS.DATA_date = ctx->date;
+ }
+
+ return QMCKL_SUCCESS;
+ }
+ #+end_src
+
The domain of quantum chemistry needs a library in which the main
diff --git a/autogen.sh b/autogen.sh
index 588976a..018116f 100755
--- a/autogen.sh
+++ b/autogen.sh
@@ -1,5 +1,5 @@
#!/bin/bash
export srcdir="."
-python ${srcdir}/tools/build_makefile.py
+python3 ${srcdir}/tools/build_makefile.py
autoreconf -i -Wall --no-recursive
diff --git a/configure.ac b/configure.ac
index 31322e1..e8e09c6 100644
--- a/configure.ac
+++ b/configure.ac
@@ -35,7 +35,7 @@
AC_PREREQ([2.69])
-AC_INIT([qmckl],[0.1.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html])
+AC_INIT([qmckl],[0.2.1],[https://github.com/TREX-CoE/qmckl/issues],[],[https://trex-coe.github.io/qmckl/index.html])
AC_CONFIG_AUX_DIR([tools])
AM_INIT_AUTOMAKE([subdir-objects color-tests parallel-tests silent-rules 1.11])
@@ -93,6 +93,7 @@ AC_PROG_F77
m4_version_prereq([2.70],[], [AC_PROG_CC_C99])
AS_IF([test "$ac_cv_prog_cc_c99" = "no"], [AC_MSG_ERROR([The compiler does not support C99])])
AC_PROG_CC_C_O
+AM_PROG_CC_C_O
AC_PROG_FC
AC_PROG_FC_C_O
AC_FC_PP_DEFINE
@@ -137,10 +138,10 @@ case "$with_chameleon" in
[PKG_CFLAGS="$PKG_CFLAGS $LIBCHAMELEON_CFLAGS"
PKG_LIBS="$PKG_LIBS $LIBCHAMELEON_LIBS"]
,[
-
+
## something went wrong.
## try to find the package without pkg-config
-
+
## check that the library is actually new enough.
## by testing for a 1.0.0+ function which we use
AC_CHECK_LIB(chameleon,CHAMELEON_finalize,[LIBCHAMELEON_LIBS="-lchameleon"])
@@ -212,6 +213,9 @@ esac
case $CC in
+ *gcc*)
+ CFLAGS="$CFLAGS -fPIC"
+ ;;
*nvc*)
CFLAGS="$CFLAGS -fPIC"
;;
@@ -224,6 +228,131 @@ AS_IF([test "$HAVE_HPC" = "yes"], [
AC_DEFINE([HAVE_HPC], [1], [If defined, activate HPC routines])
])
+# Enable Verificarlo tests
+AC_ARG_ENABLE([vfc_ci],
+[ --enable-vfc_ci Build the library with vfc_ci support],
+[case "${enableval}" in
+ yes) vfc_ci=true && FCFLAGS="-D VFC_CI $FCFLAGS" && CFLAGS="-D VFC_CI $CFLAGS";;
+ no) vfc_ci=false ;;
+ *) AC_MSG_ERROR([bad value ${enableval} for --enable_vfc_ci]) ;;
+esac],[vfc_ci=false])
+AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue])
+
+if test "$FC" = "verificarlo-f"; then
+ AC_MSG_NOTICE(verificarlo-f detected)
+ # Arguments order is important here
+ FCFLAGS="-Mpreprocess $FCFLAGS"
+fi
+
+## Enable GPU offloading
+
+# GPU offloading
+AC_ARG_ENABLE(gpu, [AS_HELP_STRING([--enable-gpu],[openmp|openacc : Use GPU-offloaded functions])], enable_gpu=$enableval, enable_gpu=no)
+AS_IF([test "$enable_gpu" = "yes"], [enable_gpu="openmp"])
+
+# OpenMP offloading
+HAVE_OPENMP_OFFLOAD="no"
+AS_IF([test "$enable_gpu" = "openmp"], [
+ AC_DEFINE([HAVE_OPENMP_OFFLOAD], [1], [If defined, activate OpenMP-offloaded routines])
+ HAVE_OPENMP_OFFLOAD="yes"
+ case $CC in
+
+ *gcc*)
+ CFLAGS="$CFLAGS -fopenmp"
+ ;;
+ *nvc*)
+ CFLAGS="$CFLAGS -mp=gpu"
+ ;;
+ esac
+
+ case $FC in
+
+ *gfortran*)
+ FCFLAGS="$FCFLAGS -fopenmp"
+ ;;
+ *nvfortran*)
+ FCFLAGS="$FCFLAGS -mp=gpu"
+ ;;
+ esac]
+)
+
+# OpenMP offloading
+HAVE_OPENACC_OFFLOAD="no"
+AS_IF([test "$enable_gpu" = "openacc"], [
+ AC_DEFINE([HAVE_OPENACC_OFFLOAD], [1], [If defined, activate OpenACC-offloaded routines])
+ HAVE_OPENACC_OFFLOAD="yes"
+ case $CC in
+
+ *gcc*)
+ CFLAGS="$CFLAGS -fopenacc"
+ ;;
+ *nvc*)
+ CFLAGS="$CFLAGS -acc=gpu"
+ ;;
+ esac
+
+ case $FC in
+
+ *gfortran*)
+ FCFLAGS="$FCFLAGS -fopenacc"
+ ;;
+ *nvfortran*)
+ FCFLAGS="$FCFLAGS -acc=gpu"
+ ;;
+ esac
+
+])
+
+# cuBLAS offloading
+AC_ARG_WITH(cublas, [AS_HELP_STRING([--with-cublas],[Use cuBLAS-offloaded functions])], HAVE_CUBLAS_OFFLOAD=$withval, HAVE_CUBLAS_OFFLOAD=no)
+AS_IF([test "$HAVE_CUBLAS_OFFLOAD" = "yes"], [
+ AC_DEFINE([HAVE_CUBLAS_OFFLOAD], [1], [If defined, activate cuBLAS-offloaded routines])
+ HAVE_OPENACC_OFFLOAD="yes"
+ case $CC in
+
+ *gcc*)
+ CFLAGS="$CFLAGS -fopenmp"
+ LDFLAGS="-lcublas"
+ ;;
+ *nvc*)
+ CFLAGS="$CFLAGS -mp=gpu -cudalib=cublas"
+ ;;
+ esac
+
+ case $FC in
+
+ *gfortran*)
+ FCFLAGS="$FCFLAGS -fopenmp"
+ ;;
+ *nvfortran*)
+ FCFLAGS="$FCFLAGS -mp=gpu -cudalib=cublas"
+ ;;
+ esac
+])
+
+AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no)
+if test "$ok" = "yes"; then
+ AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free])
+ ARGS="${ARGS} malloc-trace"
+fi
+
+AC_ARG_ENABLE(prof, [AS_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no)
+if test "$ok" = "yes"; then
+ CFLAGS="${CFLAGS} -pg"
+ AC_DEFINE(ENABLE_PROF,1,[Define when using the profiler tool])
+ ARGS="${ARGS} prof"
+fi
+
+AC_ARG_WITH(efence, [AS_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no)
+if test "$ok" = "yes"; then
+ AC_CHECK_LIB([efence], [malloc])
+ ARGS="${ARGS} efence"
+fi
+
+
+
+##
+
AC_ARG_ENABLE(debug, [AS_HELP_STRING([--enable-debug],[compile for debugging])], ok=$enableval, ok=no)
if test "$ok" = "yes"; then
if test "$GCC" = "yes"; then
@@ -247,26 +376,6 @@ if test "$ok" = "yes"; then
ARGS="${ARGS} debug"
fi
-AC_ARG_ENABLE(malloc-trace, [AS_HELP_STRING([--enable-malloc-trace],[use debug malloc/free])], ok=$enableval, ok=no)
-if test "$ok" = "yes"; then
- AC_DEFINE(MALLOC_TRACE,"malloc_trace.dat",[Define to use debugging malloc/free])
- ARGS="${ARGS} malloc-trace"
-fi
-
-AC_ARG_ENABLE(prof, [AS_HELP_STRING([--enable-prof],[compile for profiling])], ok=$enableval, ok=no)
-if test "$ok" = "yes"; then
- CFLAGS="${CFLAGS} -pg"
- AC_DEFINE(ENABLE_PROF,1,[Define when using the profiler tool])
- ARGS="${ARGS} prof"
-fi
-
-AC_ARG_WITH(efence, [AS_HELP_STRING([--with-efence],[use ElectricFence library])], ok=$withval, ok=no)
-if test "$ok" = "yes"; then
- AC_CHECK_LIB(efence, malloc)
- ARGS="${ARGS} efence"
-fi
-
-
# Checks for header files.
## qmckl
@@ -321,21 +430,6 @@ if test "x${QMCKL_DEVEL}" != "x"; then
fi
-# Enable Verificarlo tests
-AC_ARG_ENABLE([vfc_ci],
-[ --enable-vfc_ci Build the library with vfc_ci support],
-[case "${enableval}" in
- yes) vfc_ci=true && FCFLAGS="-D VFC_CI $FCFLAGS" && CFLAGS="-D VFC_CI $CFLAGS";;
- no) vfc_ci=false ;;
- *) AC_MSG_ERROR([bad value ${enableval} for --enable_vfc_ci]) ;;
-esac],[vfc_ci=false])
-AM_CONDITIONAL([VFC_CI], [test x$vfc_ci = xtrue])
-
-if test "$FC" = "verificarlo-f"; then
- AC_MSG_NOTICE(verificarlo-f detected)
- # Arguments order is important here
- FCFLAGS="-Mpreprocess $FCFLAGS"
-fi
#PKG-CONFIG
#mkl-dynamic-lp64-seq
@@ -371,6 +465,9 @@ LDFLAGS:........: ${LDFLAGS}
LIBS............: ${LIBS}
USE CHAMELEON...: ${with_chameleon}
HPC version.....: ${HAVE_HPC}
+OpenMP offload..: ${HAVE_OPENMP_OFFLOAD}
+OpenACC offload.: ${HAVE_OPENACC_OFFLOAD}
+cuBLAS offload..: ${HAVE_CUBLAS_OFFLOAD}
Package features:
${ARGS}
diff --git a/org/examples.org b/org/examples.org
deleted file mode 100644
index 58dc366..0000000
--- a/org/examples.org
+++ /dev/null
@@ -1,199 +0,0 @@
-#+TITLE: Code examples
-#+SETUPFILE: ../tools/theme.setup
-#+INCLUDE: ../tools/lib.org
-
-In this section, we present examples of usage of QMCkl.
-For simplicity, we assume that the wave function parameters are stores
-in a [[https://github.com/TREX-CoE/trexio][TREXIO]] file.
-
-* Checking errors
-
- All QMCkl functions return an error code. A convenient way to handle
- errors is to write an error-checking function that displays the
- error in text format and exits the program.
-
- #+NAME: qmckl_check_error
- #+begin_src f90
-subroutine qmckl_check_error(rc, message)
- use qmckl
- implicit none
- integer(qmckl_exit_code), intent(in) :: rc
- character(len=*) , intent(in) :: message
- character(len=128) :: str_buffer
- if (rc /= QMCKL_SUCCESS) then
- print *, message
- call qmckl_string_of_error(rc, str_buffer)
- print *, str_buffer
- call exit(rc)
- end if
-end subroutine qmckl_check_error
- #+end_src
-
-* Computing an atomic orbital on a grid
- :PROPERTIES:
- :header-args: :tangle ao_grid.f90
- :END:
-
- The following program, in Fortran, computes the values of an atomic
- orbital on a regular 3-dimensional grid. The 100^3 grid points are
- automatically defined, such that the molecule fits in a box with 5
- atomic units in the borders.
-
- This program uses the ~qmckl_check_error~ function defined above.
-
- To use this program, run
-
- #+begin_src bash :tangle no
-$ ao_grid