mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-06-18 11:15:38 +02:00
Easier context
This commit is contained in:
parent
09c8f9700f
commit
fac7e9d74f
|
|
@ -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
|
||||
|
|
|
|||
809
src/qmckl_ao.org
809
src/qmckl_ao.org
|
|
@ -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_<<filename()>>() {
|
||||
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
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -30,7 +30,7 @@ MunitResult test_<<filename()>>() {
|
|||
: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_<<filename()>>() {
|
|||
#+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
|
||||
|
|
|
|||
|
|
@ -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 <stdint.h>
|
||||
#include "qmckl_error.h"
|
||||
#include <string.h>
|
||||
#include <assert.h>
|
||||
#include <pthread.h>
|
||||
#include <errno.h>
|
||||
|
||||
#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_<<filename()>>() {
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :comments org :tangle (eval h)
|
||||
#include <errno.h>
|
||||
#include <string.h>
|
||||
#+end_src
|
||||
|
||||
*
|
||||
|
|
@ -35,10 +47,15 @@ MunitResult test_<<filename()>>() {
|
|||
|
||||
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[<<MAX_STRING_LENGTH()>>]);
|
||||
#+end_src
|
||||
|
|
@ -183,25 +200,157 @@ void qmckl_string_of_error_f(const qmckl_exit_code error, char result[<<MAX_STRI
|
|||
#+end_src
|
||||
|
||||
# Fortran interface
|
||||
#+begin_src f90 :tangle (eval fh) :noexport :noweb yes
|
||||
#+begin_src f90 :tangle (eval fh_func) :noexport :noweb yes
|
||||
interface
|
||||
subroutine qmckl_string_of_error (error, string) bind(C, name='qmckl_string_of_error_f')
|
||||
use, intrinsic :: iso_c_binding
|
||||
integer (c_int32_t), intent(in), value :: error
|
||||
import
|
||||
integer (qmckl_exit_code), intent(in), value :: error
|
||||
character, intent(out) :: string(<<MAX_STRING_LENGTH()>>)
|
||||
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
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -17,10 +17,13 @@ optimized libraries to fine-tune the memory allocation.
|
|||
#include <stdlib.h>
|
||||
#include <assert.h>
|
||||
|
||||
#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_<<filename()>>() {
|
|||
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;
|
||||
|
|
|
|||
328
src/qmckl_numprec.org
Normal file
328
src/qmckl_numprec.org
Normal file
|
|
@ -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_<<filename()>>() {
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval h_private_type)
|
||||
#ifndef QMCKL_NUMPREC_HPT
|
||||
#define QMCKL_NUMPREC_HPT
|
||||
|
||||
#include <stdint.h>
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle (eval c)
|
||||
#include <stdint.h>
|
||||
#include <assert.h>
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#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
|
||||
|
||||
|
|
@ -2,6 +2,4 @@ qmckl.org
|
|||
qmckl_error.org
|
||||
qmckl_context.org
|
||||
qmckl_memory.org
|
||||
qmckl_distance.org
|
||||
qmckl_ao.org
|
||||
test_qmckl.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}
|
||||
|
|
|
|||
|
|
@ -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}
|
||||
|
|
|
|||
|
|
@ -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")
|
||||
|
|
|
|||
|
|
@ -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"
|
||||
<<parse_table>>
|
||||
|
||||
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"
|
||||
, ""
|
||||
]
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user