diff --git a/src/qmckl.org b/src/qmckl.org index 3a568c7..a871b73 100644 --- a/src/qmckl.org +++ b/src/qmckl.org @@ -59,7 +59,6 @@ Both files are located in the =include/= directory. Moreover, within the Emacs text editor the source code blocks can be executed interactively, in the same spirit as Jupyter notebooks. - ** Source code editing For a tutorial on literate programming with org-mode, follow [[http://www.howardism.org/Technical/Emacs/literate-programming-tutorial.html][this link]]. @@ -79,7 +78,6 @@ Both files are located in the =include/= directory. Note that pandoc can be used to convert multiple markdown formats into org-mode. - ** Choice of the programming language Most of the codes of the [[https://trex-coe.eu][TREX CoE]] are written in Fortran with some scripts in @@ -111,31 +109,21 @@ Both files are located in the =include/= directory. ~iso_c_binding~ module. The name of the Fortran source files should end with =_f.f90= to be properly handled by the =Makefile=. The names of the functions defined in Fortran should be the same as those exposed in the API suffixed by - =_f=. Fortran interfaces should also be written in the =qmckl_f.f90= file. + =_f=. For more guidelines on using Fortran to generate a C interface, see [[http://fortranwiki.org/fortran/show/Generating+C+Interfaces][this link]]. ** Coding rules + + The authors should follow the recommendations of the C99 + [[https://wiki.sei.cmu.edu/confluence/display/c/SEI+CERT+C+Coding+Standard][SEI+CERT C Coding Standard]]. - The authors should follow the recommendations of the - [[https://wiki.sei.cmu.edu/confluence/display/c/SEI+CERT+C+Coding+Standard][SEI+CERT+C+Coding+Standard]]. - - - Store a new value in pointers immediately after the memory is - freed - - Free dynamically allocated memory when no longer needed - - # # TODO: decide on a coding style - - # To improve readability, we maintain a consistent coding style in - # the library. - - # - For C source files, we will use __(decide on a coding style)__ - # - For Fortran source files, we will use __(decide on a coding - # style)__ - - # Coding style can be automatically checked with [[https://clang.llvm.org/docs/ClangFormat.html][clang-format]]. + Compliance can be checked with =cppcheck= as: + #+begin_src bash +cppcheck --addon=cert --enable=all *.c &> cppcheck.out + #+end_src ** Design of the library @@ -199,14 +187,43 @@ Both files are located in the =include/= directory. The internal structure of the context is not specified, to give a maximum of freedom to the different implementations. Modifying the state is done by setters and getters, prefixed by - =qmckl_context_set_= an =qmckl_context_get_=. When a context - variable is modified by a setter, a copy of the old data structure - is made and updated, and the pointer to the new data structure is - returned, such that the old contexts can still be accessed. It is - also possible to modify the state in an impure fashion, using the - =qmckl_context_update_= functions. The context and its old - versions can be destroyed with =qmckl_context_destroy=. + =qmckl_set_= an =qmckl_get_=. +** Headers + + A single =qmckl.h= header to be distributed by the library + is built by concatenating some of the produced header files. + To facilitate building the =qmckl.h= file, we separate types from + function declarations in headers. Types should be defined in header + files suffixed by =_type.h=, and functions in files suffixed by + =_func.h=. + As these files will be concatenated in a single file, they should + not be guarded by ~#ifndef *_H~, and they should not include other + produced headers. + + Some particular types that are not exported need to be known by the + context, and there are some functions to update instances of these + types contained inside the context. For example, a + ~qmckl_numprec_struct~ is present in the context, and the function + ~qmckl_set_numprec_range~ takes a context as a parameter, and set a + value in the ~qmckl_numprec_struct~ contained in the context. + Because of these intricate dependencies, a private header is + created, containing the ~qmckl_numprec_struct~. This header is + included in the private header which defines the type of the + context. Headers for private types are suffixed by =_private_type.h= + and headers for private functions, =_private_func.h=. + Fortran interfaces should also be written in the =*_f_func.f90= file, + and the types definitions should be written in the =*_f_type.f90= file. + + | File | Scope | Description | + |--------------------+---------+------------------------------| + | =*_type.h= | Public | Type definitions | + | =*_func.h= | Public | Function definitions | + | =*_private_type.h= | Private | Type definitions | + | =*_private_func.h= | Private | Function definitions | + | =*fh_type.f90= | Public | Fortran type definitions | + | =*fh_func.f90= | Public | Fortran function definitions | + ** Low-level functions Low-level functions are very simple functions which are leaves of diff --git a/src/qmckl_ao.org b/src/qmckl_ao.org deleted file mode 100644 index e26b2d1..0000000 --- a/src/qmckl_ao.org +++ /dev/null @@ -1,809 +0,0 @@ -#+TITLE: Atomic Orbitals -#+SETUPFILE: ../docs/theme.setup - -The atomic basis set is defined as a list of shells. Each shell $s$ is -centered on a nucleus $A$, possesses a given angular momentum $l$ and a -radial function $R_s$. The radial function is a linear combination of -\emph{primitive} functions that can be of type Slater ($p=1$) or -Gaussian ($p=2$): - -\[ - R_s(\mathbf{r}) = \mathcal{N}_s |\mathbf{r}-\mathbf{R}_A|^{n_s} - \sum_{k=1}^{N_{\text{prim}}} a_{ks} - \exp \left( - \gamma_{ks} | \mathbf{r}-\mathbf{R}_A | ^p \right). -\] - -In the case of Gaussian functions, $n_s$ is always zero. -The normalization factor $\mathcal{N}_s$ ensures that all the functions -of the shell are normalized to unity. As this normalization requires -the ability to compute overlap integrals, it should be written in the -file to ensure that the file is self-contained and does not require -the client program to have the ability to compute such integrals. - -Atomic orbitals (AOs) are defined as - -\[ -\chi_i (\mathbf{r}) = P_{\eta(i)}(\mathbf{r})\, R_{\theta(i)} (\mathbf{r}) -\] - -where $\theta(i)$ returns the shell on which the AO is expanded, -and $\eta(i)$ denotes which angular function is chosen. - -In this section we describe the kernels used to compute the values, -gradients and Laplacian of the atomic basis functions. - -* Headers :noexport: - - #+NAME: filename - #+begin_src elisp tangle: no -(file-name-nondirectory (substring buffer-file-name 0 -4)) - #+end_src - - - #+begin_src c :tangle (eval c_test) :noweb yes -#include "qmckl.h" -#include "munit.h" -MunitResult test_<>() { - qmckl_context context; - context = qmckl_context_create(); - #+end_src - -* Polynomial part - -** Powers of $x-X_i$ - - The ~qmckl_ao_power~ function computes all the powers of the ~n~ - input data up to the given maximum value given in input for each of - the $n$ points: - - \[ P_{ik} = X_i^k \] - - | ~context~ | input | Global state | - | ~n~ | input | Number of values | - | ~X(n)~ | input | Array containing the input values | - | ~LMAX(n)~ | input | Array containing the maximum power for each value | - | ~P(LDP,n)~ | output | Array containing all the powers of ~X~ | - | ~LDP~ | input | Leading dimension of array ~P~ | - - Requirements: - - - ~context~ is not ~QMCKL_NULL_CONTEXT~ - - ~n~ > 0 - - ~X~ is allocated with at least $n \times 8$ bytes - - ~LMAX~ is allocated with at least $n \times 4$ bytes - - ~P~ is allocated with at least $n \times \max_i \text{LMAX}_i \times 8$ bytes - - ~LDP~ >= $\max_i$ ~LMAX[i]~ - - #+begin_src c :tangle (eval h) -qmckl_exit_code -qmckl_ao_power(const qmckl_context context, - const int64_t n, - const double *X, - const int32_t *LMAX, - const double *P, - const int64_t LDP); - #+end_src - - #+begin_src f90 :tangle (eval f) -integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info) - use qmckl - implicit none - integer*8 , intent(in) :: context - integer*8 , intent(in) :: n - real*8 , intent(in) :: X(n) - integer , intent(in) :: LMAX(n) - real*8 , intent(out) :: P(ldp,n) - integer*8 , intent(in) :: ldp - - integer*8 :: i,k - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (n <= ldp) then - info = QMCKL_INVALID_ARG_2 - return - endif - - k = MAXVAL(LMAX) - if (LDP < k) then - info = QMCKL_INVALID_ARG_6 - return - endif - - if (k <= 0) then - info = QMCKL_INVALID_ARG_4 - return - endif - - do i=1,n - P(1,i) = X(i) - do k=2,LMAX(i) - P(k,i) = P(k-1,i) * X(i) - end do - end do - -end function qmckl_ao_power_f - #+end_src - - #+begin_src f90 :tangle (eval f) :exports none -integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) & - 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 :: n - real (c_double) , intent(in) :: X(n) - integer (c_int32_t) , intent(in) :: LMAX(n) - real (c_double) , intent(out) :: P(ldp,n) - integer (c_int64_t) , intent(in) , value :: ldp - - integer, external :: qmckl_ao_power_f - info = qmckl_ao_power_f(context, n, X, LMAX, P, ldp) -end function qmckl_ao_power - #+end_src - - #+begin_src f90 :tangle (eval fh) :exports none - interface - integer(c_int32_t) function qmckl_ao_power(context, n, X, LMAX, P, ldp) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: n - integer (c_int64_t) , intent(in) , value :: ldp - real (c_double) , intent(in) :: X(n) - integer (c_int32_t) , intent(in) :: LMAX(n) - real (c_double) , intent(out) :: P(ldp,n) - end function qmckl_ao_power - end interface - #+end_src - - # Test - #+begin_src f90 :tangle (eval f_test) -integer(c_int32_t) function test_qmckl_ao_power(context) bind(C) - use qmckl - implicit none - - integer(c_int64_t), intent(in), value :: context - - integer*8 :: n, LDP - integer, allocatable :: LMAX(:) - double precision, allocatable :: X(:), P(:,:) - integer*8 :: i,j - double precision :: epsilon - - epsilon = qmckl_context_get_epsilon(context) - - n = 100; - LDP = 10; - - allocate(X(n), P(LDP,n), LMAX(n)) - - do j=1,n - X(j) = -5.d0 + 0.1d0 * dble(j) - LMAX(j) = 1 + int(mod(j, 5),4) - end do - - test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP) - if (test_qmckl_ao_power /= 0) return - - test_qmckl_ao_power = -1 - - do j=1,n - do i=1,LMAX(j) - if ( X(j)**i == 0.d0 ) then - if ( P(i,j) /= 0.d0) return - else - if ( dabs(1.d0 - P(i,j) / (X(j)**i)) > epsilon ) return - end if - end do - end do - - test_qmckl_ao_power = 0 - deallocate(X,P,LMAX) -end function test_qmckl_ao_power - #+end_src - - #+begin_src c :tangle (eval c_test) :exports none -int test_qmckl_ao_power(qmckl_context context); -munit_assert_int(0, ==, test_qmckl_ao_power(context)); - #+end_src - -** Value, Gradient and Laplacian of a polynomial - - A polynomial is centered on a nucleus $\mathbf{R}_i$ - - \[ - P_l(\mathbf{r},\mathbf{R}_i) = (x-X_i)^a (y-Y_i)^b (z-Z_i)^c - \] - - The gradients with respect to electron coordinates are - - \begin{eqnarray*} - \frac{\partial }{\partial x} P_l\left(\mathbf{r},\mathbf{R}_i \right) & - = & a (x-X_i)^{a-1} (y-Y_i)^b (z-Z_i)^c \\ - \frac{\partial }{\partial y} P_l\left(\mathbf{r},\mathbf{R}_i \right) & - = & b (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c \\ - \frac{\partial }{\partial z} P_l\left(\mathbf{r},\mathbf{R}_i \right) & - = & c (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1} \\ - \end{eqnarray*} - - and the Laplacian is - - \begin{eqnarray*} - \left( \frac{\partial }{\partial x^2} + - \frac{\partial }{\partial y^2} + - \frac{\partial }{\partial z^2} \right) P_l - \left(\mathbf{r},\mathbf{R}_i \right) & = & - a(a-1) (x-X_i)^{a-2} (y-Y_i)^b (z-Z_i)^c + \\ - && b(b-1) (x-X_i)^a (y-Y_i)^{b-1} (z-Z_i)^c + \\ - && c(c-1) (x-X_i)^a (y-Y_i)^b (z-Z_i)^{c-1}. - \end{eqnarray*} - - ~qmckl_ao_polynomial_vgl~ computes the values, gradients and - Laplacians at a given point in space, of all polynomials with an - angular momentum up to ~lmax~. - - | ~context~ | input | Global state | - | ~X(3)~ | input | Array containing the coordinates of the points | - | ~R(3)~ | input | Array containing the x,y,z coordinates of the center | - | ~lmax~ | input | Maximum angular momentum | - | ~n~ | output | Number of computed polynomials | - | ~L(ldl,n)~ | output | Contains a,b,c for all ~n~ results | - | ~ldl~ | input | Leading dimension of ~L~ | - | ~VGL(ldv,n)~ | output | Value, gradients and Laplacian of the polynomials | - | ~ldv~ | input | Leading dimension of array ~VGL~ | - - Requirements: - - - ~context~ is not ~QMCKL_NULL_CONTEXT~ - - ~n~ > 0 - - ~lmax~ >= 0 - - ~ldl~ >= 3 - - ~ldv~ >= 5 - - ~X~ is allocated with at least $3 \times 8$ bytes - - ~R~ is allocated with at least $3 \times 8$ bytes - - ~n~ >= ~(lmax+1)(lmax+2)(lmax+3)/6~ - - ~L~ is allocated with at least $3 \times n \times 4$ bytes - - ~VGL~ is allocated with at least $5 \times n \times 8$ bytes - - On output, ~n~ should be equal to ~(lmax+1)(lmax+2)(lmax+3)/6~ - - On output, the powers are given in the following order (l=a+b+c): - - Increasing values of ~l~ - - Within a given value of ~l~, alphabetical order of the - string made by a*"x" + b*"y" + c*"z" (in Python notation). - For example, with a=0, b=2 and c=1 the string is "yyz" - - # Header - #+begin_src c :tangle (eval h) -qmckl_exit_code -qmckl_ao_polynomial_vgl(const qmckl_context context, - const double *X, - const double *R, - const int32_t lmax, - const int64_t *n, - const int32_t *L, - const int64_t ldl, - const double *VGL, - const int64_t ldv); - #+end_src - - # Source - #+begin_src f90 :tangle (eval f) -integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info) - use qmckl - implicit none - integer*8 , intent(in) :: context - real*8 , intent(in) :: X(3), R(3) - integer , intent(in) :: lmax - integer*8 , intent(out) :: n - integer , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6) - integer*8 , intent(in) :: ldl - real*8 , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6) - integer*8 , intent(in) :: ldv - - integer*8 :: i,j - integer :: a,b,c,d - real*8 :: Y(3) - integer :: lmax_array(3) - real*8 :: pows(-2:lmax,3) - integer, external :: qmckl_ao_power_f - double precision :: xy, yz, xz - double precision :: da, db, dc, dd - - info = 0 - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (lmax < 0) then - info = QMCKL_INVALID_ARG_4 - return - endif - - if (ldl < 3) then - info = QMCKL_INVALID_ARG_7 - return - endif - - if (ldv < 5) then - info = QMCKL_INVALID_ARG_9 - return - endif - - - do i=1,3 - Y(i) = X(i) - R(i) - end do - - lmax_array(1:3) = lmax - if (lmax == 0) then - VGL(1,1) = 1.d0 - vgL(2:5,1) = 0.d0 - l(1:3,1) = 0 - n=1 - else if (lmax > 0) then - pows(-2:0,1:3) = 1.d0 - do i=1,lmax - pows(i,1) = pows(i-1,1) * Y(1) - pows(i,2) = pows(i-1,2) * Y(2) - pows(i,3) = pows(i-1,3) * Y(3) - end do - - VGL(1:5,1:4) = 0.d0 - l (1:3,1:4) = 0 - - VGL(1 ,1 ) = 1.d0 - vgl(1:5,2:4) = 0.d0 - - l (1,2) = 1 - vgl(1,2) = pows(1,1) - vgL(2,2) = 1.d0 - - l (2,3) = 1 - vgl(1,3) = pows(1,2) - vgL(3,3) = 1.d0 - - l (3,4) = 1 - vgl(1,4) = pows(1,3) - vgL(4,4) = 1.d0 - - n=4 - endif - - ! l>=2 - dd = 2.d0 - do d=2,lmax - da = dd - do a=d,0,-1 - db = dd-da - do b=d-a,0,-1 - c = d - a - b - dc = dd - da - db - n = n+1 - - l(1,n) = a - l(2,n) = b - l(3,n) = c - - xy = pows(a,1) * pows(b,2) - yz = pows(b,2) * pows(c,3) - xz = pows(a,1) * pows(c,3) - - vgl(1,n) = xy * pows(c,3) - - xy = dc * xy - xz = db * xz - yz = da * yz - - vgl(2,n) = pows(a-1,1) * yz - vgl(3,n) = pows(b-1,2) * xz - vgl(4,n) = pows(c-1,3) * xy - - vgl(5,n) = & - (da-1.d0) * pows(a-2,1) * yz + & - (db-1.d0) * pows(b-2,2) * xz + & - (dc-1.d0) * pows(c-2,3) * xy - - db = db - 1.d0 - end do - da = da - 1.d0 - end do - dd = dd + 1.d0 - end do - - info = QMCKL_SUCCESS - -end function qmckl_ao_polynomial_vgl_f - #+end_src - - - #+begin_src f90 :tangle (eval f) :exports none -integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) & - bind(C) result(info) - use, intrinsic :: iso_c_binding - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(in) :: X(3), R(3) - integer (c_int32_t) , intent(in) , value :: lmax - integer (c_int64_t) , intent(out) :: n - integer (c_int32_t) , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6) - integer (c_int64_t) , intent(in) , value :: ldl - real (c_double) , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6) - integer (c_int64_t) , intent(in) , value :: ldv - - integer, external :: qmckl_ao_polynomial_vgl_f - info = qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) -end function qmckl_ao_polynomial_vgl - #+end_src - - #+begin_src f90 :tangle (eval fh) :exports none - interface - integer(c_int32_t) function qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) & - bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t) , intent(in) , value :: context - integer (c_int32_t) , intent(in) , value :: lmax - integer (c_int64_t) , intent(in) , value :: ldl - integer (c_int64_t) , intent(in) , value :: ldv - real (c_double) , intent(in) :: X(3), R(3) - integer (c_int64_t) , intent(out) :: n - integer (c_int32_t) , intent(out) :: L(ldl,(lmax+1)*(lmax+2)*(lmax+3)/6) - real (c_double) , intent(out) :: VGL(ldv,(lmax+1)*(lmax+2)*(lmax+3)/6) - end function qmckl_ao_polynomial_vgl - end interface - #+end_src - - #+begin_src f90 :tangle (eval f_test) -integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C) - use qmckl - implicit none - - integer(c_int64_t), intent(in), value :: context - - integer :: lmax, d, i - integer, allocatable :: L(:,:) - integer*8 :: n, ldl, ldv, j - double precision :: X(3), R(3), Y(3) - double precision, allocatable :: VGL(:,:) - double precision :: w - double precision :: epsilon - - epsilon = qmckl_context_get_epsilon(context) - - X = (/ 1.1 , 2.2 , 3.3 /) - R = (/ 0.1 , 1.2 , -2.3 /) - Y(:) = X(:) - R(:) - - lmax = 4; - ldl = 3; - ldv = 100; - - d = (lmax+1)*(lmax+2)*(lmax+3)/6 - - allocate (L(ldl,d), VGL(ldv,d)) - - test_qmckl_ao_polynomial_vgl = & - qmckl_ao_polynomial_vgl(context, X, R, lmax, n, L, ldl, VGL, ldv) - - if (test_qmckl_ao_polynomial_vgl /= QMCKL_SUCCESS) return - if (n /= d) return - - do j=1,n - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - do i=1,3 - if (L(i,j) < 0) return - end do - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (dabs(1.d0 - VGL(1,j) / (& - Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**L(3,j) & - )) > epsilon ) return - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (L(1,j) < 1) then - if (VGL(2,j) /= 0.d0) return - else - if (dabs(1.d0 - VGL(2,j) / (& - L(1,j) * Y(1)**(L(1,j)-1) * Y(2)**L(2,j) * Y(3)**L(3,j) & - )) > epsilon ) return - end if - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (L(2,j) < 1) then - if (VGL(3,j) /= 0.d0) return - else - if (dabs(1.d0 - VGL(3,j) / (& - L(2,j) * Y(1)**L(1,j) * Y(2)**(L(2,j)-1) * Y(3)**L(3,j) & - )) > epsilon ) return - end if - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - if (L(3,j) < 1) then - if (VGL(4,j) /= 0.d0) return - else - if (dabs(1.d0 - VGL(4,j) / (& - L(3,j) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-1) & - )) > epsilon ) return - end if - - test_qmckl_ao_polynomial_vgl = QMCKL_FAILURE - w = 0.d0 - if (L(1,j) > 1) then - w = w + L(1,j) * (L(1,j)-1) * Y(1)**(L(1,j)-2) * Y(2)**L(2,j) * Y(3)**L(3,j) - end if - if (L(2,j) > 1) then - w = w + L(2,j) * (L(2,j)-1) * Y(1)**L(1,j) * Y(2)**(L(2,j)-2) * Y(3)**L(3,j) - end if - if (L(3,j) > 1) then - w = w + L(3,j) * (L(3,j)-1) * Y(1)**L(1,j) * Y(2)**L(2,j) * Y(3)**(L(3,j)-2) - end if - if (dabs(1.d0 - VGL(5,j) / w) > epsilon ) return - end do - - test_qmckl_ao_polynomial_vgl = QMCKL_SUCCESS - - deallocate(L,VGL) -end function test_qmckl_ao_polynomial_vgl - #+end_src - - #+begin_src c :tangle (eval c_test) -int test_qmckl_ao_polynomial_vgl(qmckl_context context); -munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context)); - #+end_src - -* Gaussian basis functions - - ~qmckl_ao_gaussian_vgl~ computes the values, gradients and - Laplacians at a given point of ~n~ Gaussian functions centered at - the same point: - - \[ v_i = \exp(-a_i |X-R|^2) \] - \[ \nabla_x v_i = -2 a_i (X_x - R_x) v_i \] - \[ \nabla_y v_i = -2 a_i (X_y - R_y) v_i \] - \[ \nabla_z v_i = -2 a_i (X_z - R_z) v_i \] - \[ \Delta v_i = a_i (4 |X-R|^2 a_i - 6) v_i \] - - | ~context~ | input | Global state | - | ~X(3)~ | input | Array containing the coordinates of the points | - | ~R(3)~ | input | Array containing the x,y,z coordinates of the center | - | ~n~ | input | Number of computed Gaussians | - | ~A(n)~ | input | Exponents of the Gaussians | - | ~VGL(ldv,5)~ | output | Value, gradients and Laplacian of the Gaussians | - | ~ldv~ | input | Leading dimension of array ~VGL~ | - - Requirements : - - - ~context~ is not 0 - - ~n~ > 0 - - ~ldv~ >= 5 - - ~A(i)~ > 0 for all ~i~ - - ~X~ is allocated with at least $3 \times 8$ bytes - - ~R~ is allocated with at least $3 \times 8$ bytes - - ~A~ is allocated with at least $n \times 8$ bytes - - ~VGL~ is allocated with at least $n \times 5 \times 8$ bytes - - #+begin_src c :tangle (eval h) -qmckl_exit_code -qmckl_ao_gaussian_vgl(const qmckl_context context, - const double *X, - const double *R, - const int64_t *n, - const int64_t *A, - const double *VGL, - const int64_t ldv); - #+end_src - - #+begin_src f90 :tangle (eval f) -integer function qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) result(info) - use qmckl - implicit none - integer*8 , intent(in) :: context - real*8 , intent(in) :: X(3), R(3) - integer*8 , intent(in) :: n - real*8 , intent(in) :: A(n) - real*8 , intent(out) :: VGL(ldv,5) - integer*8 , intent(in) :: ldv - - integer*8 :: i,j - real*8 :: Y(3), r2, t, u, v - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (n <= 0) then - info = QMCKL_INVALID_ARG_4 - return - endif - - if (ldv < n) then - info = QMCKL_INVALID_ARG_7 - return - endif - - - do i=1,3 - Y(i) = X(i) - R(i) - end do - r2 = Y(1)*Y(1) + Y(2)*Y(2) + Y(3)*Y(3) - - do i=1,n - VGL(i,1) = dexp(-A(i) * r2) - end do - - do i=1,n - VGL(i,5) = A(i) * VGL(i,1) - end do - - t = -2.d0 * ( X(1) - R(1) ) - u = -2.d0 * ( X(2) - R(2) ) - v = -2.d0 * ( X(3) - R(3) ) - - do i=1,n - VGL(i,2) = t * VGL(i,5) - VGL(i,3) = u * VGL(i,5) - VGL(i,4) = v * VGL(i,5) - end do - - t = 4.d0 * r2 - do i=1,n - VGL(i,5) = (t * A(i) - 6.d0) * VGL(i,5) - end do - -end function qmckl_ao_gaussian_vgl_f - #+end_src - - #+begin_src f90 :tangle (eval f) :exports none -integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) & - bind(C) result(info) - use, intrinsic :: iso_c_binding - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(in) :: X(3), R(3) - integer (c_int64_t) , intent(in) , value :: n - real (c_double) , intent(in) :: A(n) - real (c_double) , intent(out) :: VGL(ldv,5) - integer (c_int64_t) , intent(in) , value :: ldv - - integer, external :: qmckl_ao_gaussian_vgl_f - info = qmckl_ao_gaussian_vgl_f(context, X, R, n, A, VGL, ldv) -end function qmckl_ao_gaussian_vgl - #+end_src - - #+begin_src f90 :tangle (eval fh) :exports none - interface - integer(c_int32_t) function qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) & - bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: ldv - integer (c_int64_t) , intent(in) , value :: n - real (c_double) , intent(in) :: X(3), R(3), A(n) - real (c_double) , intent(out) :: VGL(ldv,5) - end function qmckl_ao_gaussian_vgl - end interface - #+end_src - - # Test - #+begin_src f90 :tangle (eval f_test) -integer(c_int32_t) function test_qmckl_ao_gaussian_vgl(context) bind(C) - use qmckl - implicit none - - integer(c_int64_t), intent(in), value :: context - - integer*8 :: n, ldv, j, i - double precision :: X(3), R(3), Y(3), r2 - double precision, allocatable :: VGL(:,:), A(:) - double precision :: epsilon - - epsilon = qmckl_context_get_epsilon(context) - - X = (/ 1.1 , 2.2 , 3.3 /) - R = (/ 0.1 , 1.2 , -2.3 /) - Y(:) = X(:) - R(:) - r2 = Y(1)**2 + Y(2)**2 + Y(3)**2 - - n = 10; - ldv = 100; - - allocate (A(n), VGL(ldv,5)) - do i=1,n - A(i) = 0.0013 * dble(ishft(1,i)) - end do - - - test_qmckl_ao_gaussian_vgl = & - qmckl_ao_gaussian_vgl(context, X, R, n, A, VGL, ldv) - if (test_qmckl_ao_gaussian_vgl /= 0) return - - test_qmckl_ao_gaussian_vgl = -1 - - do i=1,n - test_qmckl_ao_gaussian_vgl = -11 - if (dabs(1.d0 - VGL(i,1) / (& - dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -12 - if (dabs(1.d0 - VGL(i,2) / (& - -2.d0 * A(i) * Y(1) * dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -13 - if (dabs(1.d0 - VGL(i,3) / (& - -2.d0 * A(i) * Y(2) * dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -14 - if (dabs(1.d0 - VGL(i,4) / (& - -2.d0 * A(i) * Y(3) * dexp(-A(i) * r2) & - )) > epsilon ) return - - test_qmckl_ao_gaussian_vgl = -15 - if (dabs(1.d0 - VGL(i,5) / (& - A(i) * (4.d0*r2*A(i) - 6.d0) * dexp(-A(i) * r2) & - )) > epsilon ) return - end do - - test_qmckl_ao_gaussian_vgl = 0 - - deallocate(VGL) -end function test_qmckl_ao_gaussian_vgl - #+end_src - - #+begin_src c :tangle (eval c_test) :exports none -int test_qmckl_ao_gaussian_vgl(qmckl_context context); -munit_assert_int(0, ==, test_qmckl_ao_gaussian_vgl(context)); - #+end_src - -* TODO Slater basis functions - -* End of files :noexport: - -*** Test - #+begin_src c :tangle (eval c_test) - if (qmckl_context_destroy(context) != QMCKL_SUCCESS) - return QMCKL_FAILURE; - return MUNIT_OK; -} - #+end_src - -**✸ Compute file names - #+begin_src emacs-lisp -; The following is required to compute the file names - -(setq pwd (file-name-directory buffer-file-name)) -(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) -(setq f (concat pwd name "_f.f90")) -(setq fh (concat pwd name "_fh.f90")) -(setq c (concat pwd name ".c")) -(setq h (concat name ".h")) -(setq h_private (concat name "_private.h")) -(setq c_test (concat pwd "test_" name ".c")) -(setq f_test (concat pwd "test_" name "_f.f90")) - -; Minted -(require 'ox-latex) -(setq org-latex-listings 'minted) -(add-to-list 'org-latex-packages-alist '("" "listings")) -(add-to-list 'org-latex-packages-alist '("" "color")) - - #+end_src - - #+RESULTS: - | | color | - | | listings | - - -# -*- mode: org -*- -# vim: syntax=c - - diff --git a/src/qmckl_context.org b/src/qmckl_context.org index e195336..df721d5 100644 --- a/src/qmckl_context.org +++ b/src/qmckl_context.org @@ -1,24 +1,6 @@ #+TITLE: Context #+SETUPFILE: ../docs/theme.setup - - The context variable is a handle for the state of the library, - and is stored in a data structure which can't be seen outside of - the library. To simplify compatibility with other languages, the - pointer to the internal data structure is converted into a 64-bit - signed integer, defined in the ~qmckl_context~ type. - A value of ~QMCKL_NULL_CONTEXT~ for the context is equivalent to a - ~NULL~ pointer. - - #+begin_src c :comments org :tangle (eval h) -typedef int64_t qmckl_context ; -#define QMCKL_NULL_CONTEXT (qmckl_context) 0 - #+end_src - - #+begin_src f90 :comments org :tangle (eval fh) :exports none - integer*8, parameter :: QMCKL_NULL_CONTEXT = 0 - #+end_src - * Headers :noexport: #+NAME: filename @@ -33,14 +15,15 @@ typedef int64_t qmckl_context ; MunitResult test_<>() { #+end_src - #+begin_src c :tangle (eval h_private) -#ifndef __QMCKL_CONTEXT__ -#define __QMCKL_CONTEXT__ + #+begin_src c :tangle (eval h_private_type) :noweb yes +#ifndef QMCKL_CONTEXT_HPT +#define QMCKL_CONTEXT_HPT #include #include -#include "qmckl_error.h" +#include "qmckl_error_private_type.h" +#include "qmckl_numprec_private_type.h" #+end_src #+begin_src c :tangle (eval c) @@ -48,90 +31,97 @@ MunitResult test_<>() { #include #include #include -#include - -#include "qmckl_error.h" -#include "qmckl_context.h" -#include "qmckl_context_private.h" -#include "qmckl_memory.h" - #include +#include +#include + +#include "qmckl_error_type.h" +#include "qmckl_context_private_type.h" +#include "qmckl_context_type.h" + +#include "qmckl_memory_func.h" +#include "qmckl_context_func.h" + +static void init_lock(pthread_mutex_t* mutex); #+end_src * Context handling - The context appears as an immutable data structure: modifying a - context returns a new context with the modifications. Therefore, it - is necessary to store a pointer to the old version of context so - that it can be freed when necessary. - Note that we also provide a possibility to mutate the context, but - this should be done with caution, only when it is justified. + The context variable is a handle for the state of the library, + and is stored in a data structure which can't be seen outside of + the library. To simplify compatibility with other languages, the + pointer to the internal data structure is converted into a 64-bit + signed integer, defined in the ~qmckl_context~ type. + A value of ~QMCKL_NULL_CONTEXT~ for the context is equivalent to a + ~NULL~ pointer. + + #+NAME: qmckl_context + #+begin_src c :comments org :tangle (eval h_type) +typedef int64_t qmckl_context ; +#define QMCKL_NULL_CONTEXT (qmckl_context) 0 + #+end_src + + #+begin_src f90 :comments org :tangle (eval fh_type) :exports none + integer , parameter :: qmckl_context = c_int64_t + integer*8, parameter :: QMCKL_NULL_CONTEXT = 0 + #+end_src + + An immutable context would have required to implement a garbage + collector. To keep the library simple, we have chosen to implement + the context as a mutable data structure, so it has to be handled + with care. By convention, in this file ~context~ is a ~qmckl_context~ variable and ~ctx~ is a ~qmckl_context_struct*~ pointer. ** Data structure - The main data structure contains pointers to other data structures, - containing the data specific to each given domain, such that the - modified contexts don't need to duplicate the data but only the - pointers. - - #+NAME: qmckl_context_struct - #+begin_src c :comments org :tangle none :noweb yes + #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes :exports none typedef struct qmckl_context_struct { - - /* Pointer to the previous context, before modification */ - struct qmckl_context_struct * prev; - - /* Molecular system */ - qmckl_ao_basis_struct * ao_basis; - - /* To be implemented: - qmckl_nucleus_struct * nucleus; - qmckl_electron_struct * electron; - qmckl_mo_struct * mo; - qmckl_determinant_struct * det; - ,*/ - - /* Numerical precision */ - qmckl_precision_struct * fp; - - /* Error handling */ - qmckl_error_struct * error; - - /* Memory allocation */ - qmckl_memory_struct * alloc; - - /* Thread lock */ - int lock_count; - pthread_mutex_t mutex; + /* -- State of the library -- */ /* Validity checking */ - uint32_t tag; + uint64_t tag; + + /* Numerical precision */ + qmckl_numprec_struct numprec; + + /* Thread lock */ + int lock_count; + pthread_mutex_t mutex; + + /* Error handling */ + qmckl_error_struct error; + + /* Memory allocation */ + /* + qmckl_memory_struct memory; + ,*/ + + /* -- Molecular system -- */ + /* To be implemented: + qmckl_ao_basis_struct ao_basis; + + qmckl_nucleus_struct nucleus; + qmckl_electron_struct electron; + qmckl_mo_struct mo; + qmckl_determinant_struct det; + ,*/ } qmckl_context_struct; #+end_src - - #+begin_src c :comments org :tangle (eval h_private) :noweb yes :exports none -<> - -<> - -<> - -<> - -<> - #+end_src + When a new element is added to the context, the functions + [[Creation][qmckl_context_create]], [[Destroy][qmckl_context_destroy]] and [[Copy][qmckl_context_copy]] + should be updated inorder to make deep copies. + A tag is used internally to check if the memory domain pointed by a pointer is a valid context. This allows to check that even if the pointer associated with a context is non-null, we can still verify that it points to the expected data structure. - #+begin_src c :comments org :tangle (eval h_private) :noweb yes + #+begin_src c :comments org :tangle (eval h_private_type) :noweb yes #define VALID_TAG 0xBEEFFACE #define INVALID_TAG 0xDEADBEEF #+end_src @@ -140,7 +130,7 @@ typedef struct qmckl_context_struct { the pointer is a valid context. It returns the input ~qmckl_context~ if the context is valid, ~QMCKL_NULL_CONTEXT~ otherwise. - #+begin_src c :comments org :tangle (eval h) :noexport + #+begin_src c :comments org :tangle (eval h_func) :noexport qmckl_context qmckl_context_check(const qmckl_context context) ; #+end_src @@ -150,7 +140,7 @@ qmckl_context qmckl_context_check(const qmckl_context context) { if (context == QMCKL_NULL_CONTEXT) return QMCKL_NULL_CONTEXT; - const qmckl_context_struct* ctx = (qmckl_context_struct*) context; + const qmckl_context_struct* const ctx = (const qmckl_context_struct* const) context; /* Try to access memory */ if (ctx->tag != VALID_TAG) { @@ -168,7 +158,7 @@ qmckl_context qmckl_context_check(const qmckl_context context) { - It returns ~QMCKL_NULL_CONTEXT~ upon failure to allocate the internal data structure # Header - #+begin_src c :comments org :tangle (eval h) :exports none + #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_context qmckl_context_create(); #+end_src @@ -176,8 +166,8 @@ qmckl_context qmckl_context_create(); #+begin_src c :tangle (eval c) qmckl_context qmckl_context_create() { - qmckl_context_struct* ctx = - (qmckl_context_struct*) qmckl_malloc (QMCKL_NULL_CONTEXT, sizeof(qmckl_context_struct)); + qmckl_context_struct* const ctx = + (qmckl_context_struct* const) malloc (sizeof(qmckl_context_struct)); if (ctx == NULL) { return QMCKL_NULL_CONTEXT; @@ -192,7 +182,7 @@ qmckl_context qmckl_context_create() { /* Initialize data */ ctx->tag = VALID_TAG; - const qmckl_context context = (qmckl_context) ctx; + const qmckl_context context = (const qmckl_context) ctx; assert ( qmckl_context_check(context) != QMCKL_NULL_CONTEXT ); return context; @@ -200,10 +190,11 @@ qmckl_context qmckl_context_create() { #+end_src # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none + #+begin_src f90 :tangle (eval fh_func) :exports none interface - integer (c_int64_t) function qmckl_context_create() bind(C) + integer (qmckl_context) function qmckl_context_create() bind(C) use, intrinsic :: iso_c_binding + import end function qmckl_context_create end interface #+end_src @@ -217,46 +208,6 @@ munit_assert_int64( context, !=, QMCKL_NULL_CONTEXT ); munit_assert_int64( qmckl_context_check(context), ==, context ); #+end_src -** Access to the previous context - - ~qmckl_context_previous~ returns the previous context. It returns - ~QMCKL_NULL_CONTEXT~ for the initial context and for the ~NULL~ context. - - # Header - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_context qmckl_context_previous(const qmckl_context context); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -qmckl_context qmckl_context_previous(const qmckl_context context) { - - const qmckl_context checked_context = qmckl_context_check(context); - if (checked_context == (qmckl_context) 0) { - return (qmckl_context) 0; - } - - const qmckl_context_struct* ctx = (qmckl_context_struct*) checked_context; - return qmckl_context_check((qmckl_context) ctx->prev); -} - #+end_src - - # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none - interface - integer (c_int64_t) function qmckl_context_previous(context) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - end function qmckl_context_previous - end interface - #+end_src - - # Test - #+begin_src c :comments link :tangle (eval c_test) :exports none -munit_assert_int64(qmckl_context_previous(context), ==, QMCKL_NULL_CONTEXT); -munit_assert_int64(qmckl_context_previous(QMCKL_NULL_CONTEXT), ==, QMCKL_NULL_CONTEXT); - #+end_src - ** Locking For thread safety, the context may be locked/unlocked. The lock is @@ -265,16 +216,14 @@ munit_assert_int64(qmckl_context_previous(QMCKL_NULL_CONTEXT), ==, QMCKL_NULL_CO ~lock_count~ attribute. # Header - #+begin_src c :comments org :tangle (eval h) :exports none + #+begin_src c :comments org :tangle (eval h_func) :exports none void qmckl_lock (qmckl_context context); void qmckl_unlock(qmckl_context context); - -void init_lock(pthread_mutex_t* mutex); #+end_src # Source #+begin_src c :tangle (eval c) -void init_lock(pthread_mutex_t* mutex) { +static void init_lock(pthread_mutex_t* mutex) { pthread_mutexattr_t attr; int rc; @@ -289,48 +238,49 @@ void init_lock(pthread_mutex_t* mutex) { (void)pthread_mutexattr_destroy(&attr); } -void qmckl_lock(qmckl_context context) { +void qmckl_lock(const qmckl_context context) { if (context == QMCKL_NULL_CONTEXT) return ; - qmckl_context_struct *ctx = (qmckl_context_struct*) context; + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; errno = 0; int rc = pthread_mutex_lock( &(ctx->mutex) ); if (rc != 0) { - fprintf(stderr, "qmckl_lock:%s\n", strerror(rc) ); + fprintf(stderr, "DEBUG qmckl_lock:%s\n", strerror(rc) ); fflush(stderr); } assert (rc == 0); - ctx->lock_count++; + ctx->lock_count += 1; /* printf(" lock : %d\n", ctx->lock_count); */ } -void qmckl_unlock(qmckl_context context) { - qmckl_context_struct *ctx = (qmckl_context_struct*) context; +void qmckl_unlock(const qmckl_context context) { + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; int rc = pthread_mutex_unlock( &(ctx->mutex) ); if (rc != 0) { - fprintf(stderr, "qmckl_unlock:%s\n", strerror(rc) ); + fprintf(stderr, "DEBUG qmckl_unlock:%s\n", strerror(rc) ); fflush(stderr); } assert (rc == 0); - ctx->lock_count--; + ctx->lock_count -= 1; /* printf("unlock : %d\n", ctx->lock_count); */ } #+end_src -** Copy +** TODO Copy ~qmckl_context_copy~ makes a shallow copy of a context. It returns ~QMCKL_NULL_CONTEXT~ upon failure. # Header - #+begin_src c :comments org :tangle (eval h) :exports none + #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_context qmckl_context_copy(const qmckl_context context); #+end_src + # Source #+begin_src c :tangle (eval c) qmckl_context qmckl_context_copy(const qmckl_context context) { @@ -345,11 +295,11 @@ qmckl_context qmckl_context_copy(const qmckl_context context) { } - qmckl_context_struct* old_ctx = - (qmckl_context_struct*) checked_context; + const qmckl_context_struct* const old_ctx = + (qmckl_context_struct* const) checked_context; - qmckl_context_struct* new_ctx = - (qmckl_context_struct*) qmckl_malloc (context, sizeof(qmckl_context_struct)); + qmckl_context_struct* const new_ctx = + (qmckl_context_struct* const) qmckl_malloc (context, sizeof(qmckl_context_struct)); if (new_ctx == NULL) { qmckl_unlock(context); @@ -357,10 +307,10 @@ qmckl_context qmckl_context_copy(const qmckl_context context) { } /* Copy the old context on the new one */ + /* TODO Deep copies should be done here */ memcpy(new_ctx, old_ctx, sizeof(qmckl_context_struct)); - new_ctx->prev = old_ctx; - + /* As the lock was copied, both need to be unlocked */ qmckl_unlock( (qmckl_context) new_ctx ); qmckl_unlock( (qmckl_context) old_ctx ); @@ -370,11 +320,12 @@ qmckl_context qmckl_context_copy(const qmckl_context context) { #+end_src # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none + #+begin_src f90 :tangle (eval fh_func) :exports none interface - integer (c_int64_t) function qmckl_context_copy(context) bind(C) + integer (qmckl_context) function qmckl_context_copy(context) bind(C) use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context + import + integer (qmckl_context), intent(in), value :: context end function qmckl_context_copy end interface #+end_src @@ -385,7 +336,6 @@ qmckl_context new_context = qmckl_context_copy(context); munit_assert_int64(new_context, !=, QMCKL_NULL_CONTEXT); munit_assert_int64(new_context, !=, context); munit_assert_int64(qmckl_context_check(new_context), ==, new_context); -munit_assert_int64(qmckl_context_previous(new_context), ==, context); #+end_src ** Destroy @@ -394,42 +344,28 @@ munit_assert_int64(qmckl_context_previous(new_context), ==, context); It frees the context, and returns the previous context. # Header - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_context qmckl_context_destroy(qmckl_context context); + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code qmckl_context_destroy(const qmckl_context context); #+end_src # Source #+begin_src c :tangle (eval c) -qmckl_context qmckl_context_destroy(const qmckl_context context) { +qmckl_exit_code qmckl_context_destroy(const qmckl_context context) { const qmckl_context checked_context = qmckl_context_check(context); - if (checked_context == QMCKL_NULL_CONTEXT) return QMCKL_NULL_CONTEXT; + if (checked_context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; qmckl_lock(context); - qmckl_context_struct* ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); /* Shouldn't be true because the context is valid */ + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); /* Shouldn't be possible because the context is valid */ - qmckl_unlock(context); - - const qmckl_context prev_context = (qmckl_context) ctx->prev; - if (prev_context == QMCKL_NULL_CONTEXT) { - /* This is the first context, free all memory. */ - struct qmckl_memory_struct* new = NULL; - while (ctx->alloc != NULL) { - new = ctx->alloc->next; - free(ctx->alloc->pointer); - ctx->alloc->pointer = NULL; - free(ctx->alloc); - ctx->alloc = new; - } - } + /* TODO Remove all allocated data */ + /* + qmckl_memory_free_all(context); + ,*/ - qmckl_exit_code rc; - rc = qmckl_context_remove_memory(context,ctx); - assert (rc == QMCKL_SUCCESS); - - ctx->tag = INVALID_TAG; + qmckl_unlock(context); const int rc_destroy = pthread_mutex_destroy( &(ctx->mutex) ); if (rc_destroy != 0) { @@ -437,19 +373,22 @@ qmckl_context qmckl_context_destroy(const qmckl_context context) { abort(); } - rc = qmckl_free(context,ctx); + ctx->tag = INVALID_TAG; + + const qmckl_exit_code rc = qmckl_free(context,ctx); assert (rc == QMCKL_SUCCESS); - return prev_context; + return QMCKL_SUCCESS; } #+end_src # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none + #+begin_src f90 :tangle (eval fh_func) :exports none interface - integer (c_int64_t) function qmckl_context_destroy(context) bind(C) + integer (qmckl_exit_code) function qmckl_context_destroy(context) bind(C) use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context + import + integer (qmckl_context), intent(in), value :: context end function qmckl_context_destroy end interface #+end_src @@ -457,1022 +396,17 @@ qmckl_context qmckl_context_destroy(const qmckl_context context) { # Test #+begin_src c :tangle (eval c_test) :exports none munit_assert_int64(qmckl_context_check(new_context), ==, new_context); -munit_assert_int64(qmckl_context_destroy(new_context), ==, context); +munit_assert_int32(qmckl_context_destroy(new_context), ==, QMCKL_SUCCESS); munit_assert_int64(qmckl_context_check(new_context), !=, new_context); munit_assert_int64(qmckl_context_check(new_context), ==, QMCKL_NULL_CONTEXT); -munit_assert_int64(qmckl_context_destroy(context), ==, QMCKL_NULL_CONTEXT); -munit_assert_int64(qmckl_context_destroy(QMCKL_NULL_CONTEXT), ==, QMCKL_NULL_CONTEXT); +munit_assert_int32(qmckl_context_destroy(context), ==, QMCKL_SUCCESS); +munit_assert_int32(qmckl_context_destroy(QMCKL_NULL_CONTEXT), ==, QMCKL_INVALID_CONTEXT); #+end_src -* Memory allocation handling - -** Data structure - - Pointers to all allocated memory domains are stored in the context, - in a linked list. The size is also stored, to enable the - computation of the amount of currently used memory by the library. - - #+NAME: qmckl_memory_struct - #+begin_src c :comments org :tangle no -typedef struct qmckl_memory_struct { - struct qmckl_memory_struct * next ; - void * pointer ; - size_t size ; -} qmckl_memory_struct; - #+end_src - -** Append memory - - The following function, called in [[./qmckl_memory.html][=qmckl_memory.c=]], appends a new - pair (pointer, size) to the data structure. - It is forbidden to pass the ~NULL~ pointer, or a zero size. - If the context is ~QMCKL_NULL_CONTEXT~, the function returns - immediately with ~QMCKL_SUCCESS~. - - # Header - #+begin_src c :comments org :tangle (eval h_private) :exports none -qmckl_exit_code qmckl_context_append_memory(qmckl_context context, - void* pointer, - const size_t size); - #+end_src - - # Source - #+begin_src c :comments org :tangle (eval c) -qmckl_exit_code qmckl_context_append_memory(qmckl_context context, - void* pointer, - const size_t size) { - assert (pointer != NULL); - assert (size > 0L); - - qmckl_lock(context); - - if ( qmckl_context_check(context) == QMCKL_NULL_CONTEXT ) { - qmckl_unlock(context); - return QMCKL_SUCCESS; - } - - qmckl_context_struct* ctx = (qmckl_context_struct*) context; - - qmckl_memory_struct* new_alloc = (qmckl_memory_struct*) - malloc(sizeof(qmckl_memory_struct)); - - if (new_alloc == NULL) { - qmckl_unlock(context); - return QMCKL_ALLOCATION_FAILED; - } - - new_alloc->next = NULL; - new_alloc->pointer = pointer; - new_alloc->size = size; - - qmckl_memory_struct* alloc = ctx->alloc; - if (alloc == NULL) { - ctx->alloc = new_alloc; - } else { - while (alloc->next != NULL) { - alloc = alloc->next; - } - alloc->next = new_alloc; - } - - qmckl_unlock(context); - - return QMCKL_SUCCESS; - -} - #+end_src - -** Remove memory - - The following function, called in [[./qmckl_memory.html][=qmckl_memory.c=]], removes a - pointer from the data structure. - It is forbidden to pass the ~NULL~ pointer. - If the context is ~QMCKL_NULL_CONTEXT~, the function returns - immediately with ~QMCKL_SUCCESS~. - - # Header - #+begin_src c :comments org :tangle (eval h_private) :exports none -qmckl_exit_code qmckl_context_remove_memory(qmckl_context context, - const void* pointer); - #+end_src - - # Source - #+begin_src c :comments org :tangle (eval c) -qmckl_exit_code qmckl_context_remove_memory(qmckl_context context, - const void* pointer) { - assert (pointer != NULL); - - qmckl_lock(context); - - if ( qmckl_context_check(context) == QMCKL_NULL_CONTEXT ) { - qmckl_unlock(context); - return QMCKL_SUCCESS; - } - - qmckl_context_struct* ctx = (qmckl_context_struct*) context; - - qmckl_memory_struct* alloc = ctx->alloc; - qmckl_memory_struct* prev = ctx->alloc; - - while ( (alloc != NULL) && (alloc->pointer != pointer) ) { - prev = alloc; - alloc = alloc->next; - } - - if (alloc != NULL) { - prev->next = alloc->next; - free(alloc); - } - - qmckl_unlock(context); - - return QMCKL_SUCCESS; -} - #+end_src - - #+RESULTS: - -* Error handling - -** Data structure - - #+NAME: qmckl_error_struct - #+begin_src c :comments org :tangle no -#define QMCKL_MAX_FUN_LEN 256 -#define QMCKL_MAX_MSG_LEN 1024 - -typedef struct qmckl_error_struct { - - qmckl_exit_code exit_code; - char function[QMCKL_MAX_FUN_LEN]; - char message [QMCKL_MAX_MSG_LEN]; - -} qmckl_error_struct; - #+end_src - -** Updating errors - - The error is updated in the context using - ~qmckl_context_update_error~, although it is recommended to use - ~qmckl_context_set_error~ for the immutable variant. - When the error is set in the context, it is mandatory to specify - from which function the error is triggered, and a message - explaining the error. The exit code can't be ~QMCKL_SUCCESS~. - - # Header - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_exit_code -qmckl_context_update_error(qmckl_context context, - const qmckl_exit_code exit_code, - const char* function_name, - const char* message); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -qmckl_exit_code -qmckl_context_update_error(qmckl_context context, - const qmckl_exit_code exit_code, - const char* function_name, - const char* message) -{ - /* Passing a function name and a message is mandatory. */ - assert (function_name != NULL); - assert (message != NULL); - - /* Exit codes are assumed valid. */ - assert (exit_code >= 0); - assert (exit_code != QMCKL_SUCCESS); - assert (exit_code < QMCKL_INVALID_EXIT_CODE); - - qmckl_lock(context); - - /* The context is assumed to exist. */ - assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT); - - qmckl_context_struct* ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); /* Impossible because the context is valid. */ - - if (ctx->error != NULL) { - free(ctx->error); - ctx->error = NULL; - } - - qmckl_error_struct* error = - (qmckl_error_struct*) qmckl_malloc (context, sizeof(qmckl_error_struct)); - error->exit_code = exit_code; - strncpy(error->function, function_name, QMCKL_MAX_FUN_LEN); - strncpy(error->message, message, QMCKL_MAX_MSG_LEN); - - ctx->error = error; - - qmckl_unlock(context); - - return QMCKL_SUCCESS; -} - #+end_src - - The ~qmckl_context_set_error~ function returns a new context with - the error domain updated. - - # Header - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_context -qmckl_context_set_error(qmckl_context context, - const qmckl_exit_code exit_code, - const char* function_name, - const char* message); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -qmckl_context -qmckl_context_set_error(qmckl_context context, - const qmckl_exit_code exit_code, - const char* function_name, - const char* message) -{ - /* Passing a function name and a message is mandatory. */ - assert (function_name != NULL); - assert (message != NULL); - - /* Exit codes are assumed valid. */ - assert (exit_code >= 0); - assert (exit_code != QMCKL_SUCCESS); - assert (exit_code < QMCKL_INVALID_EXIT_CODE); - - /* The context is assumed to be valid */ - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) - return QMCKL_NULL_CONTEXT; - - qmckl_context new_context = qmckl_context_copy(context); - - /* Should be impossible because the context is valid */ - assert (new_context != QMCKL_NULL_CONTEXT); - - if (qmckl_context_update_error(new_context, - exit_code, - function_name, - message) != QMCKL_SUCCESS) { - return context; - } - - return new_context; -} - #+end_src - - - To make a function fail, the ~qmckl_failwith~ function should be - called, such that information about the failure is stored in - the context. The desired exit code is given as an argument, as - well as the name of the function and an error message. The return - code of the function is the desired return code. - - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_exit_code qmckl_failwith(qmckl_context context, - const qmckl_exit_code exit_code, - const char* function, - const char* message) ; - #+end_src - - #+begin_src c :comments org :tangle (eval c) -qmckl_exit_code qmckl_failwith(qmckl_context context, - const qmckl_exit_code exit_code, - const char* function, - const char* message) { - - assert (exit_code > 0); - assert (exit_code < QMCKL_INVALID_EXIT_CODE); - assert (function != NULL); - assert (message != NULL); - assert (strlen(function) < QMCKL_MAX_FUN_LEN); - assert (strlen(message) < QMCKL_MAX_MSG_LEN); - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) - return QMCKL_NULL_CONTEXT; - - const qmckl_exit_code rc = - qmckl_context_update_error(context, exit_code, function, message); - - assert (rc == QMCKL_SUCCESS); - - return exit_code; -} - - #+end_src - - For example, this function can be used as - #+begin_src c :tangle no -if (x < 0) { - return qmckl_failwith(context, - QMCKL_INVALID_ARG_2, - "qmckl_function", - "Expected x >= 0"); - } - #+end_src - - # To decode the error messages, ~qmckl_strerror~ converts an - # error code into a string. - -* Control of the numerical precision - - Controlling numerical precision enables optimizations. Here, the - default parameters determining the target numerical precision and - range are defined. - - #+NAME: table-precision - | ~QMCKL_DEFAULT_PRECISION~ | 53 | - | ~QMCKL_DEFAULT_RANGE~ | 11 | - - # We need to force Emacs not to indent the Python code: - # -*- org-src-preserve-indentation: t - -#+begin_src python :var table=table-precision :results drawer :exports result -""" This script generates the C and Fortran constants from the org-mode table. -""" - -result = [ "#+begin_src c :comments org :tangle (eval h)" ] -for (text, code) in table: - text=text.replace("~","") - result += [ f"#define {text:30s} {code:d}" ] -result += [ "#+end_src" ] - -result += [ "" ] - -result += [ "#+begin_src f90 :comments org :tangle (eval fh) :exports none" ] -for (text, code) in table: - text=text.replace("~","") - result += [ f" integer, parameter :: {text:30s} = {code:d}" ] -result += [ "#+end_src" ] - -return '\n'.join(result) - -#+end_src - -#+RESULTS: -:results: -#+begin_src c :comments org :tangle (eval h) -#define QMCKL_DEFAULT_PRECISION 53 -#define QMCKL_DEFAULT_RANGE 11 -#+end_src - -#+begin_src f90 :comments org :tangle (eval fh) :exports none - integer, parameter :: QMCKL_DEFAULT_PRECISION = 53 - integer, parameter :: QMCKL_DEFAULT_RANGE = 11 -#+end_src -:end: - - #+NAME: qmckl_precision_struct - #+begin_src c :comments org :tangle no -typedef struct qmckl_precision_struct { - int precision; - int range; -} qmckl_precision_struct; - #+end_src - - The following functions set and get the required precision and - range. ~precision~ is an integer between 2 and 53, and ~range~ is an - integer between 2 and 11. - - The setter functions functions return a new context as a 64-bit - integer. The getter functions return the value, as a 32-bit - integer. The update functions return ~QMCKL_SUCCESS~ or - ~QMCKL_FAILURE~. - -** Precision - ~qmckl_context_update_precision~ modifies the parameter for the - numerical precision in a context. If the context doesn't have any - precision set yet, the default values are used. - - # Header - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision) { - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) - return QMCKL_INVALID_CONTEXT; - - if (precision < 2) { - return qmckl_failwith(context, - QMCKL_INVALID_ARG_2, - "qmckl_context_update_precision", - "precision < 2"); - } - - if (precision > 53) { - return qmckl_failwith(context, - QMCKL_INVALID_ARG_2, - "qmckl_context_update_precision", - "precision > 53"); - } - - qmckl_context_struct* ctx = (qmckl_context_struct*) context; - - /* This should be always true */ - assert (ctx != NULL); - - qmckl_lock(context); - - if (ctx->fp == NULL) { - - ctx->fp = (qmckl_precision_struct*) - qmckl_malloc(context, sizeof(qmckl_precision_struct)); - - if (ctx->fp == NULL) { - return qmckl_failwith(context, - QMCKL_ALLOCATION_FAILED, - "qmckl_context_update_precision", - "ctx->fp"); - } - ctx->fp->range = QMCKL_DEFAULT_RANGE; - } - - ctx->fp->precision = precision; - - qmckl_unlock(context); - - return QMCKL_SUCCESS; -} - #+end_src - - # Fortran interface - - #+begin_src f90 :tangle (eval fh) - interface - integer (c_int32_t) function qmckl_context_update_precision(context, precision) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - integer (c_int32_t), intent(in), value :: precision - end function qmckl_context_update_precision - end interface - #+end_src - - ~qmckl_context_set_precision~ returns a copy of the context with a - different precision parameter. - - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision) { - qmckl_context new_context = qmckl_context_copy(context); - if (new_context == 0) return 0; - - if (qmckl_context_update_precision(new_context, precision) == QMCKL_FAILURE) return 0; - - return new_context; -} - #+end_src - - # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none - interface - integer (c_int64_t) function qmckl_context_set_precision(context, precision) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - integer (c_int32_t), intent(in), value :: precision - end function qmckl_context_set_precision - end interface - #+end_src - - ~qmckl_context_get_precision~ returns the value of the numerical precision in the context. - - #+begin_src c :comments org :tangle (eval h) :exports none -int32_t qmckl_context_get_precision(const qmckl_context context); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -int qmckl_context_get_precision(const qmckl_context context) { - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return qmckl_failwith(context, - QMCKL_INVALID_CONTEXT, - "qmckl_context_get_precision", - ""); - } - - const qmckl_context_struct* ctx = (qmckl_context_struct*) context; - if (ctx->fp != NULL) - return ctx->fp->precision; - else - return QMCKL_DEFAULT_PRECISION; -} - #+end_src - - # Fortran interface - #+begin_src f90 :tangle (eval fh) - interface - integer (c_int32_t) function qmckl_context_get_precision(context) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - end function qmckl_context_get_precision - end interface - #+end_src - -** Range - - ~qmckl_context_update_range~ modifies the parameter for the numerical range in a given context. - - # Header - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range) { - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) - return QMCKL_INVALID_CONTEXT; - - if (range < 2) { - return qmckl_failwith(context, - QMCKL_INVALID_ARG_2, - "qmckl_context_update_range", - "range < 2"); - } - - if (range > 11) { - return qmckl_failwith(context, - QMCKL_INVALID_ARG_2, - "qmckl_context_update_range", - "range > 11"); - } - - qmckl_context_struct* ctx = (qmckl_context_struct*) context; - - /* This should be always true */ - assert (ctx != NULL); - - qmckl_lock(context); - - if (ctx->fp == NULL) { - - ctx->fp = (qmckl_precision_struct*) - qmckl_malloc(context, sizeof(qmckl_precision_struct)); - - if (ctx->fp == NULL) { - return qmckl_failwith(context, - QMCKL_ALLOCATION_FAILED, - "qmckl_context_update_range", - "ctx->fp"); - } - - ctx->fp->precision = QMCKL_DEFAULT_PRECISION; - } - - ctx->fp->range = range; - - qmckl_unlock(context); - - return QMCKL_SUCCESS; -} - #+end_src - - # Fortran interface - #+begin_src f90 :tangle (eval fh) - interface - integer (c_int32_t) function qmckl_context_update_range(context, range) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - integer (c_int32_t), intent(in), value :: range - end function qmckl_context_update_range - end interface - #+end_src - - ~qmckl_context_set_range~ returns a copy of the context with a different precision parameter. - - #+begin_src c :comments org :tangle (eval h) :exports none -qmckl_context qmckl_context_set_range(const qmckl_context context, const int range); - #+end_src - - # Source - - #+begin_src c :tangle (eval c) -qmckl_context qmckl_context_set_range(const qmckl_context context, const int range) { - qmckl_context new_context = qmckl_context_copy(context); - if (new_context == 0) return 0; - - if (qmckl_context_update_range(new_context, range) == QMCKL_FAILURE) return 0; - - return new_context; -} - #+end_src - - # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none - interface - integer (c_int64_t) function qmckl_context_set_range(context, range) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - integer (c_int32_t), intent(in), value :: range - end function qmckl_context_set_range - end interface - #+end_src - - ~qmckl_context_get_range~ returns the value of the numerical range in the context. - - #+begin_src c :comments org :tangle (eval h) :exports none -int32_t qmckl_context_get_range(const qmckl_context context); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -int qmckl_context_get_range(const qmckl_context context) { - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return qmckl_failwith(context, - QMCKL_INVALID_CONTEXT, - "qmckl_context_get_range", - ""); - } - - const qmckl_context_struct* ctx = (qmckl_context_struct*) context; - if (ctx->fp != NULL) - return ctx->fp->range; - else - return QMCKL_DEFAULT_RANGE; -} - #+end_src - - # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none - interface - integer (c_int32_t) function qmckl_context_get_range(context) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - end function qmckl_context_get_range - end interface - #+end_src - -** Helper functions - - ~qmckl_context_get_epsilon~ returns $\epsilon = 2^{1-n}$ where ~n~ is the precision. - - #+begin_src c :comments org :tangle (eval h) :exports none -double qmckl_context_get_epsilon(const qmckl_context context); - #+end_src - - # Source - #+begin_src c :tangle (eval c) -double qmckl_context_get_epsilon(const qmckl_context context) { - const int precision = qmckl_context_get_precision(context); - return 1. / (double) (1L << (precision-1)); -} - #+end_src - - # Fortran interface - #+begin_src f90 :tangle (eval fh) :exports none - interface - real (c_double) function qmckl_context_get_epsilon(context) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - end function qmckl_context_get_epsilon - end interface - #+end_src - - -* TODO Basis set - - For H_2 with the following basis set, - - #+BEGIN_EXAMPLE -HYDROGEN -S 5 -1 3.387000E+01 6.068000E-03 -2 5.095000E+00 4.530800E-02 -3 1.159000E+00 2.028220E-01 -4 3.258000E-01 5.039030E-01 -5 1.027000E-01 3.834210E-01 -S 1 -1 3.258000E-01 1.000000E+00 -S 1 -1 1.027000E-01 1.000000E+00 -P 1 -1 1.407000E+00 1.000000E+00 -P 1 -1 3.880000E-01 1.000000E+00 -D 1 -1 1.057000E+00 1.0000000 - #+END_EXAMPLE - - we have: - - #+BEGIN_EXAMPLE -type = 'G' -shell_num = 12 -prim_num = 20 -SHELL_CENTER = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2] -SHELL_ANG_MOM = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D'] -SHELL_PRIM_NUM = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1] -prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20] -EXPONENT = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027, - 1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027, - 0.3258, 0.1027, 1.407, 0.388, 1.057] -COEFFICIENT = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421, - 1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822, - 0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0] - #+END_EXAMPLE - -** Data structure - - #+NAME: qmckl_ao_basis_struct - #+begin_src c :comments org :tangle no -typedef struct qmckl_ao_basis_struct { - - int64_t shell_num; - int64_t prim_num; - int64_t * shell_center; - int32_t * shell_ang_mom; - double * shell_factor; - double * exponent ; - double * coefficient ; - int64_t * shell_prim_num; - char type; - -} qmckl_ao_basis_struct; - #+end_src - -** ~qmckl_context_update_ao_basis~ - - Updates the data describing the AO basis set into the context. - - | ~type~ | Gaussian or Slater | - | ~shell_num~ | Number of shells | - | ~prim_num~ | Total number of primitives | - | ~SHELL_CENTER(shell_num)~ | Id of the nucleus on which the shell is centered | - | ~SHELL_ANG_MOM(shell_num)~ | Id of the nucleus on which the shell is centered | - | ~SHELL_FACTOR(shell_num)~ | Normalization factor for the shell | - | ~SHELL_PRIM_NUM(shell_num)~ | Number of primitives in the shell | - | ~SHELL_PRIM_INDEX(shell_num)~ | Address of the first primitive of the shelll in the ~EXPONENT~ array | - | ~EXPONENT(prim_num)~ | Array of exponents | - | ~COEFFICIENT(prim_num)~ | Array of coefficients | - - #+begin_src c :comments org :tangle (eval h) -qmckl_exit_code -qmckl_context_update_ao_basis(qmckl_context context , - const char type , - const int64_t shell_num , - const int64_t prim_num , - const int64_t * SHELL_CENTER , - const int32_t * SHELL_ANG_MOM , - const double * SHELL_FACTOR , - const int64_t * SHELL_PRIM_NUM , - const int64_t * SHELL_PRIM_INDEX, - const double * EXPONENT , - const double * COEFFICIENT); - #+end_src - -*** Source - #+begin_src c :tangle (eval c) -qmckl_exit_code -qmckl_context_update_ao_basis(qmckl_context context , const char type, - const int64_t shell_num , const int64_t prim_num, - const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM, - const double * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM, - const int64_t * SHELL_PRIM_INDEX, - const double * EXPONENT , const double * COEFFICIENT) -{ - - int64_t i; - - /* Check input */ - - if (type != 'G' && type != 'S') return QMCKL_FAILURE; - if (shell_num <= 0) return QMCKL_FAILURE; - if (prim_num <= 0) return QMCKL_FAILURE; - if (prim_num < shell_num) return QMCKL_FAILURE; - - for (i=0 ; ishell_center == NULL); - basis->shell_center = (int64_t*) qmckl_malloc (context, shell_num * sizeof(int64_t)); - if (basis->shell_center == NULL) { - qmckl_free(context, basis); - basis = NULL; - return QMCKL_FAILURE; - } - - assert (basis->shell_ang_mom == NULL); - basis->shell_ang_mom = (int32_t*) qmckl_malloc (context, shell_num * sizeof(int32_t)); - if (basis->shell_ang_mom == NULL) { - qmckl_free(context, basis->shell_center); - basis->shell_center = NULL; - qmckl_free(context, basis); - basis = NULL; - return QMCKL_FAILURE; - } - - assert (basis->shell_prim_num == NULL); - basis->shell_prim_num= (int64_t*) qmckl_malloc (context, shell_num * sizeof(int64_t)); - if (basis->shell_prim_num == NULL) { - qmckl_free(context, basis->shell_ang_mom); - basis->shell_ang_mom = NULL; - qmckl_free(context, basis->shell_center); - basis->shell_center = NULL; - qmckl_free(context, basis); - basis = NULL; - return QMCKL_FAILURE; - } - - assert (basis->shell_factor == NULL); - basis->shell_factor = (double *) qmckl_malloc (context, shell_num * sizeof(double)); - if (basis->shell_factor == NULL) { - qmckl_free(context, basis->shell_prim_num); - basis->shell_prim_num = NULL; - qmckl_free(context, basis->shell_ang_mom); - basis->shell_ang_mom = NULL; - qmckl_free(context, basis->shell_center); - basis->shell_center = NULL; - qmckl_free(context, basis); - basis = NULL; - return QMCKL_FAILURE; - } - - assert (basis->exponent == NULL); - basis->exponent = (double *) qmckl_malloc (context, prim_num * sizeof(double)); - if (basis->exponent == NULL) { - qmckl_free(context, basis->shell_factor); - basis->shell_factor = NULL; - qmckl_free(context, basis->shell_prim_num); - basis->shell_prim_num = NULL; - qmckl_free(context, basis->shell_ang_mom); - basis->shell_ang_mom = NULL; - qmckl_free(context, basis->shell_center); - basis->shell_center = NULL; - qmckl_free(context, basis); - basis = NULL; - return QMCKL_FAILURE; - } - - assert (basis->coefficient == NULL); - basis->coefficient = (double *) qmckl_malloc (context, prim_num * sizeof(double)); - if (basis->coefficient == NULL) { - qmckl_free(context, basis->exponent); - basis->exponent = NULL; - qmckl_free(context, basis->shell_factor); - basis->shell_factor = NULL; - qmckl_free(context, basis->shell_prim_num); - basis->shell_prim_num = NULL; - qmckl_free(context, basis->shell_ang_mom); - basis->shell_ang_mom = NULL; - qmckl_free(context, basis->shell_center); - basis->shell_center = NULL; - qmckl_free(context, basis); - basis = NULL; - return QMCKL_FAILURE; - } - - - /* Assign data */ - - basis->type = type; - basis->shell_num = shell_num; - basis->prim_num = prim_num; - - for (i=0 ; ishell_center [i] = SHELL_CENTER [i]; - basis->shell_ang_mom [i] = SHELL_ANG_MOM [i]; - basis->shell_prim_num[i] = SHELL_PRIM_NUM[i]; - basis->shell_factor [i] = SHELL_FACTOR [i]; - } - - for (i=0 ; iexponent [i] = EXPONENT[i]; - basis->coefficient[i] = COEFFICIENT[i]; - } - - ctx->ao_basis = basis; - return QMCKL_SUCCESS; -} - #+end_src - -*** Fortran interface - #+begin_src f90 :tangle (eval fh) - interface - integer (c_int32_t) function qmckl_context_update_ao_basis(context, & - typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, & - SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - character(c_char) , intent(in), value :: typ - integer (c_int64_t), intent(in), value :: shell_num - integer (c_int64_t), intent(in), value :: prim_num - integer (c_int64_t), intent(in) :: SHELL_CENTER(shell_num) - integer (c_int32_t), intent(in) :: SHELL_ANG_MOM(shell_num) - double precision , intent(in) :: SHELL_FACTOR(shell_num) - integer (c_int64_t), intent(in) :: SHELL_PRIM_NUM(shell_num) - integer (c_int64_t), intent(in) :: SHELL_PRIM_INDEX(shell_num) - double precision , intent(in) :: EXPONENT(prim_num) - double precision , intent(in) :: COEFFICIENT(prim_num) - end function qmckl_context_update_ao_basis - end interface - #+end_src - -*** TODO Test - -** ~qmckl_context_set_ao_basis~ - - Sets the data describing the AO basis set into the context. - - | ~type~ | Gaussian or Slater | - | ~shell_num~ | Number of shells | - | ~prim_num~ | Total number of primitives | - | ~SHELL_CENTER(shell_num)~ | Id of the nucleus on which the shell is centered | - | ~SHELL_ANG_MOM(shell_num)~ | Id of the nucleus on which the shell is centered | - | ~SHELL_FACTOR(shell_num)~ | Normalization factor for the shell | - | ~SHELL_PRIM_NUM(shell_num)~ | Number of primitives in the shell | - | ~SHELL_PRIM_INDEX(shell_num)~ | Address of the first primitive of the shelll in the ~EXPONENT~ array | - | ~EXPONENT(prim_num)~ | Array of exponents | - | ~COEFFICIENT(prim_num)~ | Array of coefficients | - - #+begin_src c :comments org :tangle (eval h) -qmckl_context -qmckl_context_set_ao_basis(const qmckl_context context , const char type, - const int64_t shell_num , const int64_t prim_num, - const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM, - const double * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM, - const int64_t * SHELL_PRIM_INDEX, - const double * EXPONENT , const double * COEFFICIENT); - #+end_src - -*** Source - #+begin_src c :tangle (eval c) -qmckl_context -qmckl_context_set_ao_basis(const qmckl_context context , const char type, - const int64_t shell_num , const int64_t prim_num, - const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM, - const double * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM, - const int64_t * SHELL_PRIM_INDEX, - const double * EXPONENT , const double * COEFFICIENT) -{ - - qmckl_context new_context = qmckl_context_copy(context); - if (new_context == 0) return 0; - - if (qmckl_context_update_ao_basis(new_context, type, shell_num, prim_num, - SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, - SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, - COEFFICIENT - ) == QMCKL_FAILURE) - return 0; - - return new_context; -} - #+end_src - -*** Fortran interface - #+begin_src f90 :tangle (eval fh) - interface - integer (c_int64_t) function qmckl_context_set_ao_basis(context, & - typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, & - SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C) - use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - character(c_char) , intent(in), value :: typ - integer (c_int64_t), intent(in), value :: shell_num - integer (c_int64_t), intent(in), value :: prim_num - integer (c_int64_t), intent(in) :: SHELL_CENTER(shell_num) - integer (c_int32_t), intent(in) :: SHELL_ANG_MOM(shell_num) - double precision , intent(in) :: SHELL_FACTOR(shell_num) - integer (c_int64_t), intent(in) :: SHELL_PRIM_NUM(shell_num) - integer (c_int64_t), intent(in) :: SHELL_PRIM_INDEX(shell_num) - double precision , intent(in) :: EXPONENT(prim_num) - double precision , intent(in) :: COEFFICIENT(prim_num) - end function qmckl_context_set_ao_basis - end interface - #+end_src - -*** TODO Test - * End of files :noexport: - #+begin_src c :comments link :tangle (eval h_private) + #+begin_src c :comments link :tangle (eval h_private_type) #endif #+end_src @@ -1482,32 +416,4 @@ return MUNIT_OK; } #+end_src -*** Compute file names - #+begin_src emacs-lisp -; The following is required to compute the file names - -(setq pwd (file-name-directory buffer-file-name)) -(setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) -(setq f (concat pwd name "_f.f90")) -(setq fh (concat pwd name "_fh.f90")) -(setq c (concat pwd name ".c")) -(setq h (concat name ".h")) -(setq h_private (concat name "_private.h")) -(setq c_test (concat pwd "test_" name ".c")) -(setq f_test (concat pwd "test_" name "_f.f90")) - -; Minted -(require 'ox-latex) -(setq org-latex-listings 'minted) -(add-to-list 'org-latex-packages-alist '("" "listings")) -(add-to-list 'org-latex-packages-alist '("" "color")) - - #+end_src - - #+RESULTS: - | | color | - | | listings | - -# -*- mode: org -*- -# vim: syntax=c - + diff --git a/src/qmckl_distance.org b/src/qmckl_distance.org index b7e264b..cad3809 100644 --- a/src/qmckl_distance.org +++ b/src/qmckl_distance.org @@ -30,7 +30,7 @@ MunitResult test_<>() { :PROPERTIES: :Name: qmckl_distance_sq :CRetType: qmckl_exit_code - :FRetType: integer + :FRetType: qmckl_exit_code :END: ~qmckl_distance_sq~ computes the matrix of the squared distances @@ -72,7 +72,7 @@ MunitResult test_<>() { #+CALL: generate_c_header(table=qmckl_distance_sq_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: - #+begin_src c :tangle (eval h) :comments org + #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_distance_sq ( const qmckl_context context, const char transa, @@ -227,14 +227,15 @@ end function qmckl_distance_sq_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_distance_sq & + integer (qmckl_exit_code) function qmckl_distance_sq & (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) & bind(C) result(info) use, intrinsic :: iso_c_binding + import implicit none - integer (c_int64_t) , intent(in) :: context + integer (qmckl_context), intent(in) :: context character , intent(in) :: transa character , intent(in) :: transb integer (c_int64_t) , intent(in) :: m @@ -246,7 +247,7 @@ end function qmckl_distance_sq_f real (c_double ) , intent(out) :: C(ldc,n) integer (c_int64_t) , intent(in) :: ldc - integer(c_int32_t), external :: qmckl_distance_sq_f + integer (qmckl_exit_code), external :: qmckl_distance_sq_f info = qmckl_distance_sq_f & (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) @@ -258,13 +259,14 @@ end function qmckl_distance_sq_f #+RESULTS: #+begin_src f90 :tangle (eval fh) :comments org :exports none interface - integer(c_int32_t) function qmckl_distance_sq & + integer (qmckl_exit_code) function qmckl_distance_sq & (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) & bind(C) use, intrinsic :: iso_c_binding + import implicit none - integer (c_int64_t) , intent(in) :: context + integer (qmckl_context), intent(in) :: context character , intent(in) :: transa character , intent(in) :: transb integer (c_int64_t) , intent(in) :: m @@ -282,10 +284,10 @@ end function qmckl_distance_sq_f *** Test :noexport: #+begin_src f90 :tangle (eval f_test) -integer(c_int32_t) function test_qmckl_distance_sq(context) bind(C) +integer(qmckl_exit_code) function test_qmckl_distance_sq(context) bind(C) use qmckl implicit none - integer(c_int64_t), intent(in), value :: context + integer(qmckl_context), intent(in), value :: context double precision, allocatable :: A(:,:), B(:,:), C(:,:) integer*8 :: m, n, LDA, LDB, LDC diff --git a/src/qmckl_error.org b/src/qmckl_error.org index 4f52854..b1fc1dc 100644 --- a/src/qmckl_error.org +++ b/src/qmckl_error.org @@ -8,9 +8,24 @@ (file-name-nondirectory (substring buffer-file-name 0 -4)) #+end_src + #+begin_src c :tangle (eval h_private_type) +#ifndef QMCKL_ERROR_HPT +#define QMCKL_ERROR_HPT + #+end_src + #+begin_src c :tangle (eval c) #include -#include "qmckl_error.h" +#include +#include +#include +#include + +#include "qmckl_error_type.h" +#include "qmckl_context_private_type.h" +#include "qmckl_context_type.h" + +#include "qmckl_context_func.h" +#include "qmckl_error_func.h" #+end_src #+begin_src c :tangle (eval c_test) :noweb yes @@ -19,9 +34,6 @@ MunitResult test_<>() { #+end_src - #+begin_src c :comments org :tangle (eval h) -#include -#include #+end_src * @@ -35,10 +47,15 @@ MunitResult test_<>() { All the functions return with an exit code, defined as #+NAME: type-exit-code - #+begin_src c :comments org :tangle (eval h) + #+begin_src c :comments org :tangle (eval h_type) typedef int32_t qmckl_exit_code; #+end_src + + #+begin_src f90 :comments org :tangle (eval fh_type) :exports none + integer , parameter :: qmckl_exit_code = c_int32_t + #+end_src + The exit code returns the completion status of the function to the calling program. When a function call completed successfully, ~QMCKL_SUCCESS~ is returned. If one of the functions of @@ -76,18 +93,18 @@ typedef int32_t qmckl_exit_code; codes from the org-mode table. """ -result = [ "#+begin_src c :comments org :tangle (eval h) :exports none" ] +result = [ "#+begin_src c :comments org :tangle (eval h_type) :exports none" ] for (text, code,_) in table: text=text.replace("~","") - result += [ f"#define {text:30s} {code:d}" ] + result += [ f"#define {text:30s} ((qmckl_exit_code) {code:d})" ] result += [ "#+end_src" ] result += [ "" ] -result += [ "#+begin_src f90 :comments org :tangle (eval fh) :exports none" ] +result += [ "#+begin_src f90 :comments org :tangle (eval fh_type) :exports none" ] for (text, code,_) in table: text=text.replace("~","") - result += [ f" integer, parameter :: {text:30s} = {code:d}" ] + result += [ f" integer(qmckl_exit_code), parameter :: {text:30s} = {code:d}" ] result += [ "#+end_src" ] return '\n'.join(result) @@ -96,44 +113,44 @@ return '\n'.join(result) #+RESULTS: :results: - #+begin_src c :comments org :tangle (eval h) :exports none - #define QMCKL_SUCCESS 0 - #define QMCKL_INVALID_ARG_1 1 - #define QMCKL_INVALID_ARG_2 2 - #define QMCKL_INVALID_ARG_3 3 - #define QMCKL_INVALID_ARG_4 4 - #define QMCKL_INVALID_ARG_5 5 - #define QMCKL_INVALID_ARG_6 6 - #define QMCKL_INVALID_ARG_7 7 - #define QMCKL_INVALID_ARG_8 8 - #define QMCKL_INVALID_ARG_9 9 - #define QMCKL_INVALID_ARG_10 10 - #define QMCKL_FAILURE 101 - #define QMCKL_ERRNO 102 - #define QMCKL_INVALID_CONTEXT 103 - #define QMCKL_ALLOCATION_FAILED 104 - #define QMCKL_DEALLOCATION_FAILED 105 - #define QMCKL_INVALID_EXIT_CODE 106 + #+begin_src c :comments org :tangle (eval h_type) :exports none + #define QMCKL_SUCCESS ((qmckl_exit_code) 0) + #define QMCKL_INVALID_ARG_1 ((qmckl_exit_code) 1) + #define QMCKL_INVALID_ARG_2 ((qmckl_exit_code) 2) + #define QMCKL_INVALID_ARG_3 ((qmckl_exit_code) 3) + #define QMCKL_INVALID_ARG_4 ((qmckl_exit_code) 4) + #define QMCKL_INVALID_ARG_5 ((qmckl_exit_code) 5) + #define QMCKL_INVALID_ARG_6 ((qmckl_exit_code) 6) + #define QMCKL_INVALID_ARG_7 ((qmckl_exit_code) 7) + #define QMCKL_INVALID_ARG_8 ((qmckl_exit_code) 8) + #define QMCKL_INVALID_ARG_9 ((qmckl_exit_code) 9) + #define QMCKL_INVALID_ARG_10 ((qmckl_exit_code) 10) + #define QMCKL_FAILURE ((qmckl_exit_code) 101) + #define QMCKL_ERRNO ((qmckl_exit_code) 102) + #define QMCKL_INVALID_CONTEXT ((qmckl_exit_code) 103) + #define QMCKL_ALLOCATION_FAILED ((qmckl_exit_code) 104) + #define QMCKL_DEALLOCATION_FAILED ((qmckl_exit_code) 105) + #define QMCKL_INVALID_EXIT_CODE ((qmckl_exit_code) 106) #+end_src - #+begin_src f90 :comments org :tangle (eval fh) :exports none - integer, parameter :: QMCKL_SUCCESS = 0 - integer, parameter :: QMCKL_INVALID_ARG_1 = 1 - integer, parameter :: QMCKL_INVALID_ARG_2 = 2 - integer, parameter :: QMCKL_INVALID_ARG_3 = 3 - integer, parameter :: QMCKL_INVALID_ARG_4 = 4 - integer, parameter :: QMCKL_INVALID_ARG_5 = 5 - integer, parameter :: QMCKL_INVALID_ARG_6 = 6 - integer, parameter :: QMCKL_INVALID_ARG_7 = 7 - integer, parameter :: QMCKL_INVALID_ARG_8 = 8 - integer, parameter :: QMCKL_INVALID_ARG_9 = 9 - integer, parameter :: QMCKL_INVALID_ARG_10 = 10 - integer, parameter :: QMCKL_FAILURE = 101 - integer, parameter :: QMCKL_ERRNO = 102 - integer, parameter :: QMCKL_INVALID_CONTEXT = 103 - integer, parameter :: QMCKL_ALLOCATION_FAILED = 104 - integer, parameter :: QMCKL_DEALLOCATION_FAILED = 105 - integer, parameter :: QMCKL_INVALID_EXIT_CODE = 106 + #+begin_src f90 :comments org :tangle (eval fh_type) :exports none + integer(qmckl_exit_code), parameter :: QMCKL_SUCCESS = 0 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_1 = 1 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_2 = 2 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_3 = 3 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_4 = 4 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_5 = 5 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_6 = 6 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_7 = 7 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_8 = 8 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_9 = 9 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_ARG_10 = 10 + integer(qmckl_exit_code), parameter :: QMCKL_FAILURE = 101 + integer(qmckl_exit_code), parameter :: QMCKL_ERRNO = 102 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_CONTEXT = 103 + integer(qmckl_exit_code), parameter :: QMCKL_ALLOCATION_FAILED = 104 + integer(qmckl_exit_code), parameter :: QMCKL_DEALLOCATION_FAILED = 105 + integer(qmckl_exit_code), parameter :: QMCKL_INVALID_EXIT_CODE = 106 #+end_src :end: @@ -144,7 +161,7 @@ return '\n'.join(result) #+NAME: MAX_STRING_LENGTH : 128 - #+begin_src c :comments org :tangle (eval h) :exports none :noweb yes + #+begin_src c :comments org :tangle (eval h_func) :exports none :noweb yes const char* qmckl_string_of_error(const qmckl_exit_code error); void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<>]); #+end_src @@ -183,25 +200,157 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<>) end subroutine qmckl_string_of_error end interface #+end_src +* Data structure in context + + The strings are declared with a maximum fixed size to avoid + dynamic memory allocation. + + #+begin_src c :comments org :tangle (eval h_private_type) +#define QMCKL_MAX_FUN_LEN 256 +#define QMCKL_MAX_MSG_LEN 1024 + +typedef struct qmckl_error_struct { + + qmckl_exit_code exit_code; + char function[QMCKL_MAX_FUN_LEN]; + char message [QMCKL_MAX_MSG_LEN]; + +} qmckl_error_struct; + #+end_src + +* Updating errors in the context + + The error is updated in the context using + ~qmckl_set_error~. + When the error is set in the context, it is mandatory to specify + from which function the error is triggered, and a message + explaining the error. The exit code can't be ~QMCKL_SUCCESS~. + + # Header + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code +qmckl_set_error(qmckl_context context, + const qmckl_exit_code exit_code, + const char* function_name, + const char* message); + #+end_src + + # Source + #+begin_src c :tangle (eval c) +qmckl_exit_code +qmckl_set_error(qmckl_context context, + const qmckl_exit_code exit_code, + const char* function_name, + const char* message) +{ + /* Passing a function name and a message is mandatory. */ + assert (function_name != NULL); + assert (message != NULL); + + /* Exit codes are assumed valid. */ + assert (exit_code >= 0); + assert (exit_code != QMCKL_SUCCESS); + assert (exit_code < QMCKL_INVALID_EXIT_CODE); + + /* The context is assumed to exist. */ + assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT); + + qmckl_lock(context); + { + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); /* Impossible because the context is valid. */ + + ctx->error.exit_code = exit_code; + strncpy(ctx->error.function, function_name, QMCKL_MAX_FUN_LEN); + strncpy(ctx->error.message, message, QMCKL_MAX_MSG_LEN); + } + qmckl_unlock(context); + + return QMCKL_SUCCESS; +} + #+end_src + +* Failing + + To make a function fail, the ~qmckl_failwith~ function should be + called, such that information about the failure is stored in + the context. The desired exit code is given as an argument, as + well as the name of the function and an error message. The return + code of the function is the desired return code. + + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code qmckl_failwith(qmckl_context context, + const qmckl_exit_code exit_code, + const char* function, + const char* message) ; + #+end_src + + #+begin_src c :comments org :tangle (eval c) +qmckl_exit_code qmckl_failwith(qmckl_context context, + const qmckl_exit_code exit_code, + const char* function, + const char* message) { + + assert (exit_code > 0); + assert (exit_code < QMCKL_INVALID_EXIT_CODE); + assert (function != NULL); + assert (message != NULL); + assert (strlen(function) < QMCKL_MAX_FUN_LEN); + assert (strlen(message) < QMCKL_MAX_MSG_LEN); + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) + return QMCKL_NULL_CONTEXT; + + const qmckl_exit_code rc = + qmckl_set_error(context, exit_code, function, message); + + assert (rc == QMCKL_SUCCESS); + + return exit_code; +} + + #+end_src + + For example, this function can be used as + #+begin_src c :tangle no +if (x < 0) { + return qmckl_failwith(context, + QMCKL_INVALID_ARG_2, + "qmckl_function", + "Expected x >= 0"); + } + #+end_src + +* TODO Decoding errors + + To decode the error messages, ~qmckl_strerror~ converts an + error code into a string. + * End of files :noexport: + #+begin_src c :comments link :tangle (eval h_private_type) +#endif + #+end_src + + ** Test - #+begin_src c :comments link :tangle (eval c_test) + #+begin_src c :comments link :tangle (eval c_test) return MUNIT_OK; } - #+end_src + #+end_src -# -*- mode: org -*- -# vim: syntax=c + # -*- mode: org -*- + # vim: syntax=c diff --git a/src/qmckl_memory.org b/src/qmckl_memory.org index 3f3a4e3..5016212 100644 --- a/src/qmckl_memory.org +++ b/src/qmckl_memory.org @@ -17,10 +17,13 @@ optimized libraries to fine-tune the memory allocation. #include #include -#include "qmckl_error.h" -#include "qmckl_context.h" -#include "qmckl_context_private.h" -#include "qmckl_memory.h" +#include "qmckl_error_type.h" +#include "qmckl_context_type.h" +#include "qmckl_context_private_type.h" + +#include "qmckl_memory_func.h" +#include "qmckl_context_func.h" +#include "qmckl_error_func.h" #+end_src #+begin_src c :tangle (eval c_test) :noweb yes @@ -44,7 +47,7 @@ MunitResult test_<>() { If the allocation failed, the ~NULL~ pointer is returned. # Header - #+begin_src c :tangle (eval h) :noexport + #+begin_src c :tangle (eval h_func) :noexport void* qmckl_malloc(qmckl_context context, const size_t size); #+end_src @@ -56,25 +59,29 @@ void* qmckl_malloc(qmckl_context context, #+begin_src c :tangle (eval c) void* qmckl_malloc(qmckl_context context, const size_t size) { + assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT); void * pointer = calloc(size, (size_t) 1); + /* if (qmckl_context_check(context) != QMCKL_NULL_CONTEXT) { qmckl_exit_code rc; rc = qmckl_context_append_memory(context, pointer, size); assert (rc == QMCKL_SUCCESS); } + */ return pointer; } #+end_src # Fortran interface - #+begin_src f90 :tangle (eval fh) :noexport + #+begin_src f90 :tangle (eval fh_func) :noexport interface type (c_ptr) function qmckl_malloc (context, size) bind(C) use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context - integer (c_int64_t), intent(in), value :: size + import + integer (qmckl_context), intent(in), value :: context + integer (c_int64_t) , intent(in), value :: size end function qmckl_malloc end interface #+end_src @@ -95,16 +102,17 @@ a[2] = 3; munit_assert_int(a[2], ==, 3); case some important information has been stored related to memory allocation and needs to be updated. - #+begin_src c :tangle (eval h) + #+begin_src c :tangle (eval h_func) qmckl_exit_code qmckl_free(qmckl_context context, void *ptr); #+end_src - #+begin_src f90 :tangle (eval fh) + #+begin_src f90 :tangle (eval fh_func) interface - integer (c_int32_t) function qmckl_free (context, ptr) bind(C) + integer (qmckl_exit_code) function qmckl_free (context, ptr) bind(C) use, intrinsic :: iso_c_binding - integer (c_int64_t), intent(in), value :: context + import + integer (qmckl_context), intent(in), value :: context type (c_ptr), intent(in), value :: ptr end function qmckl_free end interface @@ -122,10 +130,12 @@ qmckl_exit_code qmckl_free(qmckl_context context, void *ptr) { "NULL pointer"); } + /* qmckl_exit_code rc; rc = qmckl_context_remove_memory(context, ptr); assert (rc == QMCKL_SUCCESS); + */ } free(ptr); return QMCKL_SUCCESS; diff --git a/src/qmckl_numprec.org b/src/qmckl_numprec.org new file mode 100644 index 0000000..8349529 --- /dev/null +++ b/src/qmckl_numprec.org @@ -0,0 +1,328 @@ +#+TITLE: Numerical precision +#+SETUPFILE: ../docs/theme.setup + +* Headers :noexport: + + #+NAME: filename + #+begin_src elisp tangle: no +(file-name-nondirectory (substring buffer-file-name 0 -4)) + #+end_src + + + #+begin_src c :tangle (eval c_test) :noweb yes +#include "qmckl.h" +#include "munit.h" +MunitResult test_<>() { + #+end_src + + #+begin_src c :tangle (eval h_private_type) +#ifndef QMCKL_NUMPREC_HPT +#define QMCKL_NUMPREC_HPT + +#include + #+end_src + + #+begin_src c :tangle (eval c) +#include +#include +#include +#include +#include + +#include "qmckl_error_type.h" +#include "qmckl_context_type.h" +#include "qmckl_context_private_type.h" +#include "qmckl_numprec_type.h" + +#include "qmckl_numprec_func.h" +#include "qmckl_error_func.h" +#include "qmckl_context_func.h" + + #+end_src + +* Control of the numerical precision + + Controlling numerical precision enables optimizations. Here, the + default parameters determining the target numerical precision and + range are defined. Following the IEEE Standard for Floating-Point + Arithmetic (IEEE 754), + /precision/ refers to the number of significand bits and /range/ + refers to the number of exponent bits. + + #+NAME: table-precision + | ~QMCKL_DEFAULT_PRECISION~ | 53 | + | ~QMCKL_DEFAULT_RANGE~ | 11 | + + # We need to force Emacs not to indent the Python code: + # -*- org-src-preserve-indentation: t + +#+begin_src python :var table=table-precision :results drawer :exports results +""" This script generates the C and Fortran constants from the org-mode table. +""" + +result = [ "#+begin_src c :comments org :tangle (eval h_type)" ] +for (text, code) in table: + text=text.replace("~","") + result += [ f"#define {text:30s} {code:d}" ] +result += [ "#+end_src" ] + +result += [ "" ] + +result += [ "#+begin_src f90 :comments org :tangle (eval fh_func) :exports none" ] +for (text, code) in table: + text=text.replace("~","") + result += [ f" integer, parameter :: {text:30s} = {code:d}" ] +result += [ "#+end_src" ] + +return '\n'.join(result) + +#+end_src + +#+RESULTS: +:results: +#+begin_src c :comments org :tangle (eval h_type) +#define QMCKL_DEFAULT_PRECISION 53 +#define QMCKL_DEFAULT_RANGE 11 +#+end_src + +#+begin_src f90 :comments org :tangle (eval fh_func) :exports none + integer, parameter :: QMCKL_DEFAULT_PRECISION = 53 + integer, parameter :: QMCKL_DEFAULT_RANGE = 11 +#+end_src +:end: + + #+begin_src c :comments org :tangle (eval h_private_type) +typedef struct qmckl_numprec_struct { + uint32_t precision; + uint32_t range; +} qmckl_numprec_struct; + #+end_src + + The following functions set and get the required precision and + range. ~precision~ is an integer between 2 and 53, and ~range~ is an + integer between 2 and 11. + + The setter functions functions return a new context as a 64-bit + integer. The getter functions return the value, as a 32-bit + integer. The update functions return ~QMCKL_SUCCESS~ or + ~QMCKL_FAILURE~. + +* Precision + ~qmckl_context_set_numprec_precision~ modifies the parameter for the + numerical precision in the context. + + # Header + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const int precision); + #+end_src + + # Source + #+begin_src c :tangle (eval c) +qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const int precision) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) + return QMCKL_INVALID_CONTEXT; + + if (precision < 2) { + return qmckl_failwith(context, + QMCKL_INVALID_ARG_2, + "qmckl_update_numprec_precision", + "precision < 2"); + } + + if (precision > 53) { + return qmckl_failwith(context, + QMCKL_INVALID_ARG_2, + "qmckl_update_numprec_precision", + "precision > 53"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + + /* This should be always true because the context is valid */ + assert (ctx != NULL); + + qmckl_lock(context); + { + ctx->numprec.precision = (uint32_t) precision; + } + qmckl_unlock(context); + + return QMCKL_SUCCESS; +} + #+end_src + + # Fortran interface + + #+begin_src f90 :tangle (eval fh_func) + interface + integer (qmckl_exit_code) function qmckl_set_numprec_precision(context, precision) bind(C) + use, intrinsic :: iso_c_binding + import + integer (qmckl_context), intent(in), value :: context + integer (c_int32_t), intent(in), value :: precision + end function qmckl_set_numprec_precision + end interface + #+end_src + + ~qmckl_get_numprec_precision~ returns the value of the numerical precision in the context. + + #+begin_src c :comments org :tangle (eval h_func) :exports none +int32_t qmckl_get_numprec_precision(const qmckl_context context); + #+end_src + + # Source + #+begin_src c :tangle (eval c) +int qmckl_get_numprec_precision(const qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith(context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_numprec_precision", + ""); + } + + const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + return ctx->numprec.precision; +} + #+end_src + + # Fortran interface + #+begin_src f90 :tangle (eval fh_func) + interface + integer (qmckl_exit_code) function qmckl_get_numprec_precision(context) bind(C) + use, intrinsic :: iso_c_binding + import + integer (qmckl_context), intent(in), value :: context + end function qmckl_get_numprec_precision + end interface + #+end_src + +* Range + + ~qmckl_set_numprec_range~ modifies the parameter for the numerical + range in a given context. + + # Header + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int range); + #+end_src + + # Source + #+begin_src c :tangle (eval c) +qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int range) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) + return QMCKL_INVALID_CONTEXT; + + if (range < 2) { + return qmckl_failwith(context, + QMCKL_INVALID_ARG_2, + "qmckl_set_numprec_range", + "range < 2"); + } + + if (range > 11) { + return qmckl_failwith(context, + QMCKL_INVALID_ARG_2, + "qmckl_set_numprec_range", + "range > 11"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + + /* This should be always true because the context is valid */ + assert (ctx != NULL); + + qmckl_lock(context); + { + ctx->numprec.range = (uint32_t) range; + } + qmckl_unlock(context); + + return QMCKL_SUCCESS; +} + #+end_src + + # Fortran interface + #+begin_src f90 :tangle (eval fh_func) + interface + integer (qmckl_exit_code) function qmckl_numprec_set_range(context, range) bind(C) + use, intrinsic :: iso_c_binding + import + integer (qmckl_context), intent(in), value :: context + integer (c_int32_t), intent(in), value :: range + end function qmckl_numprec_set_range + end interface + #+end_src + + ~qmckl_get_numprec_range~ returns the value of the numerical range in the context. + + #+begin_src c :comments org :tangle (eval h_func) :exports none +int32_t qmckl_context_get_range(const qmckl_context context); + #+end_src + + # Source + #+begin_src c :tangle (eval c) +int qmckl_get_numprec_range(const qmckl_context context) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith(context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_numprec_range", + ""); + } + + const qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + return ctx->numprec.range; +} + #+end_src + + # Fortran interface + #+begin_src f90 :tangle (eval fh_func) :exports none + interface + integer (qmckl_exit_code) function qmckl_get_numprec_range(context) bind(C) + use, intrinsic :: iso_c_binding + import + integer (qmckl_context), intent(in), value :: context + end function qmckl_get_numprec_range + end interface + #+end_src + +* Helper functions + + ~qmckl_context_get_epsilon~ returns $\epsilon = 2^{1-n}$ where ~n~ is the precision. + + #+begin_src c :comments org :tangle (eval h_func) :exports none +double qmckl_get_numprec_epsilon(const qmckl_context context); + #+end_src + + # Source + #+begin_src c :tangle (eval c) +double qmckl_get_numprec_epsilon(const qmckl_context context) { + const int precision = qmckl_get_numprec_precision(context); + return 1. / (double) (1L << (precision-1)); +} + #+end_src + + # Fortran interface + #+begin_src f90 :tangle (eval fh_func) :exports none + interface + real (c_double) function qmckl_get_numprec_epsilon(context) bind(C) + use, intrinsic :: iso_c_binding + import + integer (qmckl_context), intent(in), value :: context + end function qmckl_get_numprec_epsilon + end interface + #+end_src + +* End of files :noexport: + + #+begin_src c :comments link :tangle (eval h_private_type) +#endif + #+end_src + +*** Test + #+begin_src c :comments link :tangle (eval c_test) +return MUNIT_OK; +} + #+end_src + diff --git a/src/table_of_contents b/src/table_of_contents index a6a9011..03c9622 100644 --- a/src/table_of_contents +++ b/src/table_of_contents @@ -2,6 +2,4 @@ qmckl.org qmckl_error.org qmckl_context.org qmckl_memory.org -qmckl_distance.org -qmckl_ao.org test_qmckl.org diff --git a/tools/Building.org b/tools/Building.org index 32ccf3f..959deb8 100644 --- a/tools/Building.org +++ b/tools/Building.org @@ -353,7 +353,12 @@ EOF HEADERS="" for i in $(cat table_of_contents) do - HEADERS+="${i%.org}.h " + HEADERS+="${i%.org}_type.h " +done + +for i in $(cat table_of_contents) +do + HEADERS+="${i%.org}_func.h " done #+end_src @@ -389,7 +394,8 @@ EOF Generate Fortran interface file from all =qmckl_*_fh.f90= files #+begin_src bash -HEADERS="qmckl_*_fh.f90" +HEADERS_TYPE="qmckl_*_fh_type.f90" +HEADERS="qmckl_*_fh_func.f90" OUTPUT="../include/qmckl_f.f90" cat << EOF > ${OUTPUT} @@ -400,6 +406,11 @@ module qmckl use, intrinsic :: iso_c_binding EOF +for i in ${HEADERS_TYPE} +do + cat $i >> ${OUTPUT} +done + for i in ${HEADERS} do cat $i >> ${OUTPUT} diff --git a/tools/build_qmckl_h.sh b/tools/build_qmckl_h.sh index 78f2665..ff57803 100755 --- a/tools/build_qmckl_h.sh +++ b/tools/build_qmckl_h.sh @@ -19,7 +19,12 @@ HEADERS="" for i in $(cat table_of_contents) do - HEADERS+="${i%.org}.h " + HEADERS+="${i%.org}_type.h " +done + +for i in $(cat table_of_contents) +do + HEADERS+="${i%.org}_func.h " done @@ -96,7 +101,8 @@ EOF # Generate Fortran interface file from all =qmckl_*_fh.f90= files -HEADERS="qmckl_*_fh.f90" +HEADERS_TYPE="qmckl_*_fh_type.f90" +HEADERS="qmckl_*_fh_func.f90" OUTPUT="../include/qmckl_f.f90" cat << EOF > ${OUTPUT} @@ -146,6 +152,11 @@ module qmckl use, intrinsic :: iso_c_binding EOF +for i in ${HEADERS_TYPE} +do + cat $i >> ${OUTPUT} +done + for i in ${HEADERS} do cat $i >> ${OUTPUT} diff --git a/tools/config_tangle.el b/tools/config_tangle.el index 6f1eed1..91c122f 100755 --- a/tools/config_tangle.el +++ b/tools/config_tangle.el @@ -36,10 +36,13 @@ (setq pwd (file-name-directory buffer-file-name)) (setq name (file-name-nondirectory (substring buffer-file-name 0 -4))) (setq f (concat pwd name "_f.f90")) -(setq fh (concat pwd name "_fh.f90")) +(setq fh_func (concat pwd name "_fh_func.f90")) +(setq fh_type (concat pwd name "_fh_type.f90")) (setq c (concat pwd name ".c")) -(setq h (concat name ".h")) -(setq h_private (concat name "_private.h")) +(setq h_func (concat name "_func.h")) +(setq h_type (concat name "_type.h")) +(setq h_private_type (concat name "_private_type.h")) +(setq h_private_func (concat name "_private_func.h")) (setq c_test (concat pwd "test_" name ".c")) (setq f_test (concat pwd "test_" name "_f.f90")) (org-babel-lob-ingest "../tools/lib.org") diff --git a/tools/lib.org b/tools/lib.org index c59eaae..e6e5216 100644 --- a/tools/lib.org +++ b/tools/lib.org @@ -10,7 +10,7 @@ #+RESULTS: get_value * Table of function arguments - + #+NAME: test | qmckl_context | context | in | Global state | | char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed | @@ -24,32 +24,35 @@ | double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ | | int64_t | ldc | in | Leading dimension of array ~C~ | - + ** Fortran-C type conversions #+NAME:f_of_c #+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none" f_of_c_d = { '' : '' - , 'qmckl_context' : 'integer (c_int64_t)' - , 'int32_t' : 'integer (c_int32_t)' - , 'int64_t' : 'integer (c_int64_t)' - , 'float' : 'real (c_float )' - , 'double' : 'real (c_double )' - , 'char' : 'character' + , 'qmckl_context' : 'integer (qmckl_context)' + , 'qmckl_exit_code' : 'integer (qmckl_exit_code)' + , 'int32_t' : 'integer (c_int32_t)' + , 'int64_t' : 'integer (c_int64_t)' + , 'float' : 'real (c_float )' + , 'double' : 'real (c_double )' + , 'char' : 'character' } #+END_SRC - + #+NAME:c_of_f #+BEGIN_SRC python :var table=test :var rettyp="integer" :var fname=[] :results value :noweb yes :wrap "src f90 :tangle (eval f) :comments org :exports none" ctypeid_d = { '' : '' - , 'integer' : 'integer(c_int32_t)' - , 'integer*8' : 'integer(c_int64_t)' - , 'real' : 'real(c_float)' - , 'real*8' : 'real(c_double)' - , 'character' : 'character(c_char)' + , 'qmckl_context' : 'integer (qmckl_context)' + , 'qmckl_exit_code' : 'integer (qmckl_exit_code)' + , 'integer' : 'integer(c_int32_t)' + , 'integer*8' : 'integer(c_int64_t)' + , 'real' : 'real(c_float)' + , 'real*8' : 'real(c_double)' + , 'character' : 'character(c_char)' } #+END_SRC - + ** Parse the table #+NAME: parse_table @@ -79,7 +82,7 @@ def parse_table(table): else: d["name"] = d["name"].split('[')[0].strip() d["dims"] = [ x.replace(']','').strip() for x in dims[1:] ] - + result.append(d) return result @@ -88,7 +91,7 @@ def parse_table(table): ** Generates a C header #+NAME: generate_c_header - #+BEGIN_SRC python :var table=[] :var rettyp=[] :var fname=[] :results drawer :noweb yes :wrap "src c :tangle (eval h) :comments org" + #+BEGIN_SRC python :var table=[] :var rettyp=[] :var fname=[] :results drawer :noweb yes :wrap "src c :tangle (eval h_func) :comments org" <> results = [] @@ -105,7 +108,7 @@ for d in parse_table(table): const = "const" else: const = " " - + results += [ f" {const} {c_type} {name}" ] results=',\n'.join(results) @@ -116,11 +119,11 @@ return template #+END_SRC #+RESULTS: generate_c_header - #+begin_src c :tangle (eval h) :comments org + #+begin_src c :tangle (eval h_func) :comments org [] [] ( - ); + ); #+end_src - + ** Generates a C interface to the Fortran function #+NAME: generate_c_interface @@ -137,8 +140,9 @@ rettyp_c = ctypeid_d[rettyp.lower()] results = [ f"{rettyp_c} function {fname} &" , f" ({args}) &" , " bind(C) result(info)" -, "" +, "" , " use, intrinsic :: iso_c_binding" +, " import" , " implicit none" , "" ] @@ -166,7 +170,7 @@ for d in parse_table(table): results += [ "" , f" {rettyp_c}, external :: {fname}_f" , f" info = {fname}_f &" -, f" ({args})" +, f" ({args})" , "" , f"end function {fname}" ] @@ -193,6 +197,7 @@ results = [ f"interface" , f" ({args}) &" , " bind(C)" , " use, intrinsic :: iso_c_binding" +, " import" , " implicit none" , "" ]