diff --git a/README.html b/README.html index d399ca4..0a3c2d5 100644 --- a/README.html +++ b/README.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + QMCkl source code documentation @@ -315,8 +315,8 @@ for the JavaScript code in this tag.
  • Error handling
  • Context
  • Memory management
  • -
  • Inter-particle distances
  • Atomic Orbitals
  • +
  • Inter-particle distances
  • Testing
  • @@ -348,7 +348,7 @@ and bug reports should be submitted at

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/index.html b/index.html index d399ca4..0a3c2d5 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + QMCkl source code documentation @@ -315,8 +315,8 @@ for the JavaScript code in this tag.
  • Error handling
  • Context
  • Memory management
  • -
  • Inter-particle distances
  • Atomic Orbitals
  • +
  • Inter-particle distances
  • Testing
  • @@ -348,7 +348,7 @@ and bug reports should be submitted at

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/qmckl.html b/qmckl.html index 7672341..4aa2f47 100644 --- a/qmckl.html +++ b/qmckl.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Introduction @@ -333,30 +333,30 @@ for the JavaScript code in this tag.

    Table of Contents

    -
    -

    1 Using QMCkl

    +
    +

    1 Using QMCkl

    The qmckl.h header file has to be included in C codes when @@ -385,12 +385,12 @@ Both files are located in the include/ directory.

    -
    -

    2 Developing in QMCkl

    +
    +

    2 Developing in QMCkl

    -
    -

    2.1 Literate programming

    +
    +

    2.1 Literate programming

    In a traditional source code, most of the lines of source files of a program @@ -435,9 +435,8 @@ interactively, in the same spirit as Jupyter notebooks.

    - -
    -

    2.2 Source code editing

    +
    +

    2.2 Source code editing

    For a tutorial on literate programming with org-mode, follow this link. @@ -468,9 +467,8 @@ org-mode.

    - -
    -

    2.3 Choice of the programming language

    +
    +

    2.3 Choice of the programming language

    Most of the codes of the TREX CoE are written in Fortran with some scripts in @@ -508,7 +506,7 @@ The Fortran source files should provide a C interface using the 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.

    @@ -518,25 +516,27 @@ For more guidelines on using Fortran to generate a C interface, see

    -
    -

    2.4 Coding rules

    +
    +

    2.4 Coding rules

    -The authors should follow the recommendations of the -SEI+CERT+C+Coding+Standard. +The authors should follow the recommendations of the C99 +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
    • -
    +

    +Compliance can be checked with cppcheck as: +

    + +
    +
    cppcheck --addon=cert --enable=all *.c &> cppcheck.out
    +
    +
    - -
    -

    2.5 Design of the library

    +
    +

    2.5 Design of the library

    The proposed API should allow the library to: deal with memory transfers @@ -547,8 +547,8 @@ functions (see below).

    -
    -

    2.6 Naming conventions

    +
    +

    2.6 Naming conventions

    To avoid namespace collisions, we use qmckl_ as a prefix for all exported @@ -573,8 +573,8 @@ form is allowed.

    -
    -

    2.7 Application programming interface

    +
    +

    2.7 Application programming interface

    In the C language, the number of bits used by the integer types can change @@ -595,8 +595,7 @@ than C, we restrict the data types in the interfaces to the following: and terminated by a '\0' character (C convention).

  • Complex numbers can be represented by an array of 2 floats.
  • Boolean variables are stored as integers, 1 for true and 0 for false
  • -
  • Floating point variables should be by default
  • -
  • double unless explicitly mentioned
  • +
  • Floating point variables should be by default double unless explicitly mentioned
  • integers used for counting should always be int64_t
  • @@ -607,15 +606,15 @@ bindings in other languages in other repositories.
    -
    -

    2.8 Global state

    +
    +

    2.8 Global state

    Global variables should be avoided in the library, because it is possible that one single program needs to use multiple instances of the library. To solve this problem we propose to use a pointer to a context variable, built by the library with the -qmckl_context_create function. The =context= contains the global +qmckl_context_create function. The =context= contains the global state of the library, and is used as the first argument of many QMCkl functions.

    @@ -624,36 +623,119 @@ QMCkl functions. 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_.

    -
    -

    2.9 Low-level functions

    +
    +

    2.9 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. +

    + + + + +++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    FileScopeDescription
    *_type.hPublicType definitions
    *_func.hPublicFunction definitions
    *_private_type.hPrivateType definitions
    *_private_func.hPrivateFunction definitions
    *fh_type.f90PublicFortran type definitions
    *fh_func.f90PublicFortran function definitions
    +
    +
    + +
    +

    2.10 Low-level functions

    +
    +

    Low-level functions are very simple functions which are leaves of the function call tree (they don't call any other QMCkl function).

    These functions are pure, and unaware of the QMCkl -context. They are not allowed to allocate/deallocate memory, and +context. They are not allowed to allocate/deallocate memory, and if they need temporary memory it should be provided in input.

    -
    -

    2.10 High-level functions

    -
    +
    +

    2.11 High-level functions

    +

    High-level functions are at the top of the function call tree. They are able to choose which lower-level function to call @@ -665,28 +747,28 @@ temporary storage, to simplify the use of accelerators.

    The high-level functions should be pure, unless the introduction of non-purity is justified. All the side effects should be made in -the context variable. +the context variable.

    -
    -

    2.11 Numerical precision

    -
    +
    +

    2.12 Numerical precision

    +

    The number of bits of precision required for a function should be given as an input of low-level computational functions. This input will be used to define the values of the different thresholds that might be used to avoid computing unnecessary noise. High-level -functions will use the precision specified in the context +functions will use the precision specified in the context variable.

    -
    -

    2.12 Algorithms

    -
    +
    +

    2.13 Algorithms

    +

    Reducing the scaling of an algorithm usually implies also reducing its arithmetic complexity (number of flops per byte). Therefore, @@ -697,23 +779,11 @@ implemented adapted to different problem sizes.

    - -
    -

    2.13 Rules for the API

    -
    -
      -
    • stdint should be used for integers (int32_t, int64_t)
    • -
    • integers used for counting should always be int64_t
    • -
    • floats should be by default double, unless explicitly mentioned
    • -
    • pointers are converted to int64_t to increase portability
    • -
    -
    -

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/qmckl_ao.html b/qmckl_ao.html index 04fa775..b290f25 100644 --- a/qmckl_ao.html +++ b/qmckl_ao.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Atomic Orbitals @@ -333,77 +333,56 @@ for the JavaScript code in this tag.

    Table of Contents

    -

    -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. -

    - -
    -

    1 Polynomial part

    +
    +

    1 Context

    -
    -
    -

    1.1 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 \] +The following arrays are stored in the context:

    @@ -418,47 +397,1122 @@ the \(n\) points: - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + + + + + + + + + + + + + + + + + + + + + + + + +
    contextinputGlobal statetype Gaussian ('G') or Slater ('S')
    ninputNumber of valuesshell_num Number of shells
    X(n)inputArray containing the input valuesprim_num Total number of primitives
    LMAX(n)inputArray containing the maximum power for each valueshell_center[shell_num]Id of the nucleus on which each shell is centered
    P(LDP,n)outputArray containing all the powers of Xshell_ang_mom[shell_num]Angular momentum of each shell
    LDPinputLeading dimension of array Pshell_prim_num[shell_num]Number of primitives in each shell
    shell_prim_index[shell_num]Address of the first primitive of each shell in the EXPONENT array
    shell_factor[shell_num]Normalization factor for each shell
    exponent[prim_num]Array of exponents
    coefficient[prim_num]Array of coefficients

    -Requirements: +For H2 with the following basis set,

    +
    +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
    +
    + +

    +we have: +

    + +
    +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_factor  = [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]
    +shell_prim_num = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
    +shell_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]
    +
    +
    + +
    +

    1.1 Data structure

    +
    +
    +
    typedef struct qmckl_ao_basis_struct {
    +  int32_t   provided;
    +  int32_t   uninitialized;
    +  int64_t   shell_num;
    +  int64_t   prim_num;
    +  int64_t * shell_center;
    +  char    * shell_ang_mom;
    +  int64_t * shell_prim_num;
    +  int64_t * shell_prim_index;
    +  double  * shell_factor;
    +  double  * exponent    ;
    +  double  * coefficient ;
    +  char      type;
    +} qmckl_ao_basis_struct;
    +
    +
    + +

    +The uninitialized integer contains one bit set to one for each +initialization function which has not bee called. It becomes equal +to zero after all initialization functions have been called. The +struct is then initialized and provided == 1. +

    +
    +
    + +
    +

    1.2 Access functions

    +
    +

    +Access to scalars copies the values at the passed address, and +for array values a pointer to the array is returned. +

    + +
    +
    char      qmckl_get_ao_basis_type             (const qmckl_context context);
    +int64_t   qmckl_get_ao_basis_shell_num        (const qmckl_context context);
    +int64_t   qmckl_get_ao_basis_prim_num         (const qmckl_context context);
    +int64_t*  qmckl_get_ao_basis_shell_center     (const qmckl_context context);
    +char*     qmckl_get_ao_basis_shell_ang_mom    (const qmckl_context context);
    +int64_t*  qmckl_get_ao_basis_shell_prim_num   (const qmckl_context context);
    +int64_t*  qmckl_get_ao_basis_shell_prim_index (const qmckl_context context);
    +double*   qmckl_get_ao_basis_shell_factor     (const qmckl_context context);
    +double*   qmckl_get_ao_basis_exponent         (const qmckl_context context);
    +double*   qmckl_get_ao_basis_coefficient      (const qmckl_context context);
    +
    +
    + +
    +
    int32_t   qmckl_ao_basis_provided           (const qmckl_context context);
    +
    +
    + +
    +
    if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +  return NULL;
    +}
    +
    +
    + +
    +
    char qmckl_get_ao_basis_type (const qmckl_context context) {
    +
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return (char) 0;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return (char) 0;
    +  }
    +
    +  assert (ctx->ao_basis.type != (char) 0);
    +  return ctx->ao_basis.type;
    +}
    +
    +
    +int64_t qmckl_get_ao_basis_shell_num (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return (int64_t) 0;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 1;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return (int64_t) 0;
    +  }
    +
    +  assert (ctx->ao_basis.shell_num > (int64_t) 0);
    +  return ctx->ao_basis.shell_num;
    +}
    +
    +
    +int64_t qmckl_get_ao_basis_prim_num (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return (int64_t) 0;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 2;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return (int64_t) 0;
    +  }
    +
    +  assert (ctx->ao_basis.prim_num > (int64_t) 0);
    +  return ctx->ao_basis.prim_num;
    +}
    +
    +
    +int64_t* qmckl_get_ao_basis_shell_center (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return NULL;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 3;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return NULL;
    +  }
    +
    +  assert (ctx->ao_basis.shell_center != NULL);
    +  return ctx->ao_basis.shell_center;
    +}
    +
    +
    +char* qmckl_get_ao_basis_shell_ang_mom (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return NULL;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 4;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return NULL;
    +  }
    +
    +  assert (ctx->ao_basis.shell_ang_mom != NULL);
    +  return ctx->ao_basis.shell_ang_mom;
    +}
    +
    +
    +int64_t* qmckl_get_ao_basis_shell_prim_num (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return NULL;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 5;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return NULL;
    +  }
    +
    +  assert (ctx->ao_basis.shell_prim_num != NULL);
    +  return ctx->ao_basis.shell_prim_num;
    +}
    +
    +
    +int64_t* qmckl_get_ao_basis_shell_prim_index (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return NULL;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 6;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return NULL;
    +  }
    +
    +  assert (ctx->ao_basis.shell_prim_index != NULL);
    +  return ctx->ao_basis.shell_prim_index;
    +}
    +
    +
    +double* qmckl_get_ao_basis_shell_factor (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return NULL;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 7;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return NULL;
    +  }
    +
    +  assert (ctx->ao_basis.shell_factor != NULL);
    +  return ctx->ao_basis.shell_factor;
    +}
    +
    +
    +double* qmckl_get_ao_basis_exponent (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return NULL;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +
    +  int32_t mask = 1 << 8;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return NULL;
    +  }
    +
    +  assert (ctx->ao_basis.exponent != NULL);
    +  return ctx->ao_basis.exponent;
    +}
    +
    +
    +double* qmckl_get_ao_basis_coefficient (const qmckl_context context) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return NULL;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  int32_t mask = 1 << 9;
    +
    +  if ( (ctx->ao_basis.uninitialized & mask) != 0) {
    +    return NULL;
    +  }
    +
    +  assert (ctx->ao_basis.coefficient != NULL);
    +  return ctx->ao_basis.coefficient;
    +}
    +
    +
    +int32_t qmckl_ao_basis_provided(const qmckl_context context) {
    +
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return 0;
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);
    +
    +  return ctx->ao_basis.provided;
    +}
    +
    +
    +
    +
    + +
    +

    1.3 Initialization functions

    +
    +

    +To set the basis set, all the following functions need to be +called. When +

    + +
    +
    qmckl_exit_code  qmckl_set_ao_basis_type             (qmckl_context context, const char      t);
    +qmckl_exit_code  qmckl_set_ao_basis_shell_num        (qmckl_context context, const int64_t   shell_num);
    +qmckl_exit_code  qmckl_set_ao_basis_prim_num         (qmckl_context context, const int64_t   prim_num);
    +qmckl_exit_code  qmckl_set_ao_basis_shell_prim_index (qmckl_context context, const int64_t * shell_prim_index);
    +qmckl_exit_code  qmckl_set_ao_basis_shell_center     (qmckl_context context, const int64_t * shell_center);
    +qmckl_exit_code  qmckl_set_ao_basis_shell_ang_mom    (qmckl_context context, const char    * shell_ang_mom);
    +qmckl_exit_code  qmckl_set_ao_basis_shell_prim_num   (qmckl_context context, const int64_t * shell_prim_num);
    +qmckl_exit_code  qmckl_set_ao_basis_shell_factor     (qmckl_context context, const double  * shell_factor);
    +qmckl_exit_code  qmckl_set_ao_basis_exponent         (qmckl_context context, const double  * exponent);
    +qmckl_exit_code  qmckl_set_ao_basis_coefficient      (qmckl_context context, const double  * coefficient);
    +
    +
    + +
    +
    if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +  return QMCKL_NULL_CONTEXT;
    + }
    +
    +qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +
    + +
    +
    ctx->ao_basis.uninitialized &= ~mask; 
    +
    +if (ctx->ao_basis.uninitialized == 0) {
    +  ctx->ao_basis.provided = 1;
    +}
    +
    +return QMCKL_SUCCESS;
    +
    +
    + + +
    +
    qmckl_exit_code qmckl_set_ao_basis_type(qmckl_context context, const char t) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  if (t != 'G' && t != 'S') {
    +    return qmckl_failwith( context,
    +                           QMCKL_INVALID_ARG_2,
    +                           "qmckl_set_ao_basis_type",
    +                           NULL);
    +  }
    +
    +  int32_t mask = 1;
    +  ctx->ao_basis.type = t;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +qmckl_exit_code qmckl_set_ao_basis_shell_num(qmckl_context context, const int64_t shell_num) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  if (shell_num <= 0) {
    +    return qmckl_failwith( context,
    +                           QMCKL_INVALID_ARG_2,
    +                           "qmckl_set_ao_basis_shell_num",
    +                           "shell_num <= 0");
    +  }
    +
    +  int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
    +
    +  if (0L < prim_num && prim_num < shell_num) {
    +    return qmckl_failwith( context,
    +                           QMCKL_INVALID_ARG_2,
    +                           "qmckl_set_ao_basis_shell_num",
    +                           "shell_num > prim_num");
    +  }
    +
    +  int32_t mask = 1 << 1;
    +  ctx->ao_basis.shell_num = shell_num;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +qmckl_exit_code  qmckl_set_ao_basis_prim_num(qmckl_context context, const int64_t prim_num) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  if (prim_num <= 0) {
    +    return qmckl_failwith( context,
    +                           QMCKL_INVALID_ARG_2,
    +                           "qmckl_set_ao_basis_shell_num",
    +                           "prim_num must be positive");
    +  }
    +
    +  int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
    +
    +  if (prim_num < shell_num) {
    +    return qmckl_failwith( context,
    +                           QMCKL_INVALID_ARG_2,
    +                           "qmckl_set_ao_basis_shell_num",
    +                           "prim_num < shell_num");
    +  }
    +
    +  int32_t mask = 1 << 2;
    +  ctx->ao_basis.prim_num = prim_num;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +qmckl_exit_code  qmckl_set_ao_basis_shell_center(qmckl_context context, const int64_t* shell_center) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  int32_t mask = 1 << 3;
    +
    +  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
    +  if (shell_num == 0L) {
    +    return qmckl_failwith( context,
    +                           QMCKL_FAILURE,
    +                           "qmckl_set_ao_basis_shell_center",
    +                           "shell_num is not set");
    +  }
    +
    +  if (ctx->ao_basis.shell_center != NULL) {
    +    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_center);
    +    if (rc != QMCKL_SUCCESS) {
    +      return qmckl_failwith( context, rc, 
    +                             "qmckl_set_ao_basis_shell_center",
    +                             NULL);
    +    }
    +  }
    +
    +  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    +  mem_info.size = shell_num * sizeof(int64_t);
    +  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
    +
    +  if (new_array == NULL) {
    +    return qmckl_failwith( context,
    +                           QMCKL_ALLOCATION_FAILED,
    +                           "qmckl_set_ao_basis_shell_center",
    +                           NULL);
    +  }
    +
    +  memcpy(new_array, shell_center, mem_info.size);
    +
    +  ctx->ao_basis.shell_center = new_array;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +qmckl_exit_code  qmckl_set_ao_basis_shell_ang_mom(qmckl_context context, const char* shell_ang_mom) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  int32_t mask = 1 << 4;
    +
    +  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
    +  if (shell_num == 0L) {
    +    return qmckl_failwith( context,
    +                           QMCKL_FAILURE,
    +                           "qmckl_set_ao_basis_shell_ang_mom",
    +                           "shell_num is not set");
    +  }
    +
    +  if (ctx->ao_basis.shell_ang_mom != NULL) {
    +    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_ang_mom);
    +    if (rc != QMCKL_SUCCESS) {
    +      return qmckl_failwith( context, rc, 
    +                             "qmckl_set_ao_basis_shell_ang_mom",
    +                             NULL);
    +    }
    +  }
    +
    +
    +  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    +  mem_info.size = shell_num * sizeof(char);
    +  char* new_array = (char*) qmckl_malloc(context, mem_info);
    +
    +  if (new_array == NULL) {
    +    return qmckl_failwith( context,
    +                           QMCKL_ALLOCATION_FAILED,
    +                           "qmckl_set_ao_basis_shell_ang_mom",
    +                           NULL);
    +  }
    +
    +  memcpy(new_array, shell_ang_mom, mem_info.size);
    +
    +  ctx->ao_basis.shell_ang_mom = new_array;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +qmckl_exit_code  qmckl_set_ao_basis_shell_prim_num(qmckl_context context, const int64_t* shell_prim_num) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  int32_t mask = 1 << 5;
    +
    +  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
    +  if (shell_num == 0L) {
    +    return qmckl_failwith( context,
    +                           QMCKL_FAILURE,
    +                           "qmckl_set_ao_basis_shell_prim_num",
    +                           "shell_num is not set");
    +  }
    +
    +  if (ctx->ao_basis.shell_prim_num != NULL) {
    +    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_num);
    +    if (rc != QMCKL_SUCCESS) {
    +      return qmckl_failwith( context, rc, 
    +                             "qmckl_set_ao_basis_shell_prim_num",
    +                             NULL);
    +    }
    +  }
    +
    +
    +  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    +  mem_info.size = shell_num * sizeof(int64_t);
    +  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
    +
    +  if (new_array == NULL) {
    +    return qmckl_failwith( context,
    +                           QMCKL_ALLOCATION_FAILED,
    +                           "qmckl_set_ao_basis_shell_prim_num",
    +                           NULL);
    +  }
    +
    +  memcpy(new_array, shell_prim_num, mem_info.size);
    +
    +  ctx->ao_basis.shell_prim_num = new_array;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +qmckl_exit_code  qmckl_set_ao_basis_shell_prim_index(qmckl_context context, const int64_t* shell_prim_index) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  int32_t mask = 1 << 6;
    +
    +  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
    +  if (shell_num == 0L) {
    +    return qmckl_failwith( context,
    +                           QMCKL_FAILURE,
    +                           "qmckl_set_ao_basis_shell_prim_index",
    +                           "shell_num is not set");
    +  }
    +
    +  if (ctx->ao_basis.shell_prim_index != NULL) {
    +    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_prim_index);
    +    if (rc != QMCKL_SUCCESS) {
    +      return qmckl_failwith( context, rc, 
    +                             "qmckl_set_ao_basis_shell_prim_index",
    +                             NULL);
    +    }
    +  }
    +
    +  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    +  mem_info.size = shell_num * sizeof(int64_t);
    +  int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info);
    +
    +  if (new_array == NULL) {
    +    return qmckl_failwith( context,
    +                           QMCKL_ALLOCATION_FAILED,
    +                           "qmckl_set_ao_basis_shell_prim_index",
    +                           NULL);
    +  }
    +
    +  memcpy(new_array, shell_prim_index, mem_info.size);
    +
    +  ctx->ao_basis.shell_prim_index = new_array;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +qmckl_exit_code  qmckl_set_ao_basis_shell_factor(qmckl_context context, const double* shell_factor) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  int32_t mask = 1 << 7;
    +
    +  const int64_t shell_num = qmckl_get_ao_basis_shell_num(context);
    +  if (shell_num == 0L) {
    +    return qmckl_failwith( context,
    +                           QMCKL_FAILURE,
    +                           "qmckl_set_ao_basis_shell_factor",
    +                           "shell_num is not set");
    +  }
    +
    +
    +  if (ctx->ao_basis.shell_factor != NULL) {
    +    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.shell_factor);
    +    if (rc != QMCKL_SUCCESS) {
    +      return qmckl_failwith( context, rc, 
    +                             "qmckl_set_ao_basis_shell_factor",
    +                             NULL);
    +    }
    +  }
    +
    +  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    +  mem_info.size = shell_num * sizeof(double);
    +  double* new_array = (double*) qmckl_malloc(context, mem_info);
    +
    +  if (new_array == NULL) {
    +    return qmckl_failwith( context,
    +                           QMCKL_ALLOCATION_FAILED,
    +                           "qmckl_set_ao_basis_shell_factor",
    +                           NULL);
    +  }
    +
    +  memcpy(new_array, shell_factor, mem_info.size);
    +
    +  ctx->ao_basis.shell_factor = new_array;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +qmckl_exit_code  qmckl_set_ao_basis_exponent(qmckl_context context, const double* exponent) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  int32_t mask = 1 << 8;
    +
    +  const int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
    +  if (prim_num == 0L) {
    +    return qmckl_failwith( context,
    +                           QMCKL_FAILURE,
    +                           "qmckl_set_ao_basis_exponent",
    +                           "prim_num is not set");
    +  }
    +
    +  if (ctx->ao_basis.exponent != NULL) {
    +    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.exponent);
    +    if (rc != QMCKL_SUCCESS) {
    +      return qmckl_failwith( context, rc, 
    +                             "qmckl_set_ao_basis_exponent",
    +                             NULL);
    +    }
    +  }
    +
    +  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    +  mem_info.size = prim_num * sizeof(double);
    +  double* new_array = (double*) qmckl_malloc(context, mem_info);
    +
    +  if (new_array == NULL) {
    +    return qmckl_failwith( context,
    +                           QMCKL_ALLOCATION_FAILED,
    +                           "qmckl_set_ao_basis_exponent",
    +                           NULL);
    +  }
    +
    +  memcpy(new_array, exponent, mem_info.size);
    +
    +  ctx->ao_basis.exponent = new_array;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +qmckl_exit_code  qmckl_set_ao_basis_coefficient(qmckl_context context, const double* coefficient) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    +    return QMCKL_NULL_CONTEXT;
    +   }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  int32_t mask = 1 << 9;
    +
    +  const int64_t prim_num = qmckl_get_ao_basis_prim_num(context);
    +  if (prim_num == 0L) {
    +    return qmckl_failwith( context,
    +                           QMCKL_FAILURE,
    +                           "qmckl_set_ao_basis_coefficient",
    +                           "prim_num is not set");
    +  }
    +
    +  if (ctx->ao_basis.coefficient != NULL) {
    +    qmckl_exit_code rc = qmckl_free(context, ctx->ao_basis.coefficient);
    +    if (rc != QMCKL_SUCCESS) {
    +      return qmckl_failwith( context, rc, 
    +                             "qmckl_set_ao_basis_coefficient",
    +                             NULL);
    +    }
    +  }
    +
    +  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    +  mem_info.size = prim_num * sizeof(double);
    +  double* new_array = (double*) qmckl_malloc(context, mem_info);
    +
    +  if (new_array == NULL) {
    +    return qmckl_failwith( context,
    +                           QMCKL_ALLOCATION_FAILED,
    +                           "qmckl_set_ao_basis_coefficient",
    +                           NULL);
    +  }
    +
    +  memcpy(new_array, coefficient, mem_info.size);
    +
    +  ctx->ao_basis.coefficient = new_array;
    +
    +  ctx->ao_basis.uninitialized &= ~mask; 
    +
    +  if (ctx->ao_basis.uninitialized == 0) {
    +    ctx->ao_basis.provided = 1;
    +  }
    +
    +  return QMCKL_SUCCESS;
    +}
    +
    +
    +
    +
    +
    + + +
    +

    1.4 Fortran interfaces

    +
    + + + +++ ++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    qmcklcontextcontextinGlobal state
    int64tninNumber of values
    doubleX[n]inArray containing the input values
    int32tLMAX[n]inArray containing the maximum power for each value
    doubleP[n][ldp]outArray containing all the powers of X
    int64tldpinLeading dimension of array P
    +
    +
    + +
    +

    1.5 Test

    +
    +
    +
    /* Reference input data */
    +
    +char    typ = 'G';
    +#define shell_num   ((int64_t) 12)
    +#define prim_num    ((int64_t) 20)
    +
    +int64_t shell_center [shell_num] =
    +  { 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2 };
    +
    +char shell_ang_mom [shell_num] =
    +  { 'S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D' };
    +
    +double shell_factor [shell_num] =
    +  { 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1. };
    +
    +int64_t shell_prim_num [shell_num] =
    +  {5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1};
    +
    +int64_t shell_prim_index [shell_num] =
    +  {1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20};
    +
    +double exponent [prim_num] =
    +  { 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 };
    +
    +double coefficient [prim_num] =
    +  { 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 };
    +/* --- */
    +
    +qmckl_exit_code rc;
    +
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_type (context, typ);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_shell_num (context, shell_num);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_prim_num (context, prim_num);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_shell_center (context, shell_center);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_shell_factor  (context, shell_factor);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_shell_center  (context, shell_prim_num);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_exponent      (context, exponent);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 0);
    +
    +rc = qmckl_set_ao_basis_coefficient   (context, coefficient);
    +munit_assert_int64(rc, ==, QMCKL_SUCCESS);
    +munit_assert_int(qmckl_ao_basis_provided(context), ==, 1);
    +
    +
    +
    +
    +
    +
    + +
    +

    2 Polynomial part

    +
    +
    +
    +

    2.1 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 \] +

    + + + + +++ ++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    qmcklcontextcontextinGlobal state
    int64tninNumber of values
    doubleX[n]inArray containing the input values
    int32tLMAX[n]inArray containing the maximum power for each value
    doubleP[n][ldp]outArray containing all the powers of X
    int64tldpinLeading dimension of array P
    +
    + +
    +

    2.1.1 Requirements:

    +
    • context is not QMCKL_NULL_CONTEXT
    • n > 0
    • @@ -467,18 +1521,28 @@ Requirements:
    • P is allocated with at least \(n \times \max_i \text{LMAX}_i \times 8\) bytes
    • LDP >= \(\max_i\) LMAX[i]
    - -
    -
    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);
    -
    +
    +
    +

    2.1.2 C Header

    +
    +
    +
    qmckl_exit_code qmckl_ao_power (
    +              const qmckl_context context,
    +      const int64_t n,
    +      const double* X,
    +      const int32_t* LMAX,
    +      double* const P,
    +      const int64_t ldp ); 
    +
    +
    +
    +
    + +
    +

    2.1.3 Source

    +
    integer function qmckl_ao_power_f(context, n, X, LMAX, P, ldp) result(info)
       use qmckl
    @@ -525,13 +1589,25 @@ Requirements:
     end function qmckl_ao_power_f
     
    +
    +
    +
    +

    2.1.4 C interface

    +
    +
    +

    2.1.5 Fortran interface

    +
    + +
    +

    2.1.6 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(qmckl_context), intent(in), value :: context
     
       integer*8                     :: n, LDP 
       integer, allocatable          :: LMAX(:) 
    @@ -539,7 +1615,7 @@ Requirements:
       integer*8                     :: i,j
       double precision              :: epsilon
     
    -  epsilon = qmckl_context_get_epsilon(context)
    +  epsilon = qmckl_get_numprec_epsilon(context)
     
       n = 100;
       LDP = 10;
    @@ -552,9 +1628,9 @@ Requirements:
       end do
     
       test_qmckl_ao_power = qmckl_ao_power(context, n, X, LMAX, P, LDP) 
    -  if (test_qmckl_ao_power /= 0) return
    +  if (test_qmckl_ao_power /= QMCKL_SUCCESS) return
     
    -  test_qmckl_ao_power = -1
    +  test_qmckl_ao_power = QMCKL_FAILURE
     
       do j=1,n
          do i=1,LMAX(j)
    @@ -566,17 +1642,18 @@ Requirements:
          end do
       end do
     
    -  test_qmckl_ao_power = 0
    +  test_qmckl_ao_power = QMCKL_SUCCESS
       deallocate(X,P,LMAX)
     end function test_qmckl_ao_power
     
    +
    -
    -

    1.2 Value, Gradient and Laplacian of a polynomial

    -
    +
    +

    2.2 Value, Gradient and Laplacian of a polynomial

    +

    A polynomial is centered on a nucleus \(\mathbf{R}_i\)

    @@ -620,7 +1697,7 @@ Laplacians at a given point in space, of all polynomials with an angular momentum up to lmax.

    - +
    @@ -628,69 +1705,80 @@ angular momentum up to lmax. ++ - - + + + - - + + + - - + + + - - + + + - - + + + - - + + + - - + + + - - + + + - - + + +
    contextinputqmcklcontextcontextin Global state
    X(3)inputdoubleX[3]in Array containing the coordinates of the points
    R(3)inputdoubleR[3]in Array containing the x,y,z coordinates of the center
    lmaxinputint32tlmaxin Maximum angular momentum
    noutputint64tninout Number of computed polynomials
    L(ldl,n)outputint32tL[n][ldl]out Contains a,b,c for all n results
    ldlinputint64tldlin Leading dimension of L
    VGL(ldv,n)outputdoubleVGL[n][ldv]out Value, gradients and Laplacian of the polynomials
    ldvinputint64tldvin Leading dimension of array VGL
    +
    -

    -Requirements: -

    - +
    +

    2.2.1 Requirements:

    +
    • context is not QMCKL_NULL_CONTEXT
    • n > 0
    • @@ -711,21 +1799,31 @@ 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"
    - -
    -
    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);
    -
    +
    +
    +

    2.2.2 C Header

    +
    +
    +
    qmckl_exit_code qmckl_ao_polynomial_vgl (
    +              const qmckl_context context,
    +      const double* X,
    +      const double* R,
    +      const int32_t lmax,
    +      int64_t* n,
    +      int32_t* const L,
    +      const int64_t ldl,
    +      double* const VGL,
    +      const int64_t ldv ); 
    +
    +
    +
    +
    + +
    +

    2.2.3 Source

    +
    integer function qmckl_ao_polynomial_vgl_f(context, X, R, lmax, n, L, ldl, VGL, ldv) result(info)
       use qmckl
    @@ -856,8 +1954,20 @@ For example, with a=0, b=2 and c=1 the string is "yyz"
     end function qmckl_ao_polynomial_vgl_f
     
    +
    +
    +
    +

    2.2.4 C interface

    +
    +
    +

    2.2.5 Fortran interface

    +
    + +
    +

    2.2.6 Test

    +
    integer(c_int32_t) function test_qmckl_ao_polynomial_vgl(context) bind(C)
       use qmckl
    @@ -873,7 +1983,7 @@ For example, with a=0, b=2 and c=1 the string is "yyz"
       double precision              :: w
       double precision              :: epsilon
     
    -  epsilon = qmckl_context_get_epsilon(context)
    +  epsilon = qmckl_get_numprec_epsilon(context)
     
       X = (/ 1.1 , 2.2 ,  3.3 /)
       R = (/ 0.1 , 1.2 , -2.3 /)
    @@ -959,10 +2069,15 @@ munit_assert_int(0, ==, test_qmckl_ao_polynomial_vgl(context));
     
    +
    -
    -

    2 Gaussian basis functions

    -
    +
    +

    3 Radial part

    +
    +
    +
    +

    3.1 Gaussian basis functions

    +

    qmckl_ao_gaussian_vgl computes the values, gradients and Laplacians at a given point of n Gaussian functions centered at @@ -1135,7 +2250,7 @@ Requirements : double precision, allocatable :: VGL(:,:), A(:) double precision :: epsilon - epsilon = qmckl_context_get_epsilon(context) + epsilon = qmckl_get_numprec_epsilon(context) X = (/ 1.1 , 2.2 , 3.3 /) R = (/ 0.1 , 1.2 , -2.3 /) @@ -1193,13 +2308,21 @@ Requirements :

    -
    -

    3 TODO Slater basis functions

    +
    +

    3.2 TODO Slater basis functions

    +
    + +
    +

    3.3 TODO Radial functions on a grid

    +
    +
    +
    +

    4 Combining radial and polynomial parts

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/qmckl_context.html b/qmckl_context.html index 9b8f21c..7af6afa 100644 --- a/qmckl_context.html +++ b/qmckl_context.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Context @@ -299,28 +299,6 @@ for the JavaScript code in this tag. } /*]]>*///--> - - + +
    +

    1 Context handling

    +

    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 +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 @@ -396,21 +338,16 @@ A value of QMCKL_NULL_CONTEXT for the context is equivalent to a

    -
    typedef int64_t qmckl_context ;
    +
    typedef int64_t qmckl_context ;
     #define QMCKL_NULL_CONTEXT (qmckl_context) 0
     
    -
    -

    1 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. +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.

    @@ -419,52 +356,15 @@ and ctx is a qmckl_context_struct* pointer.

    -
    -

    1.1 Data structure

    +
    +

    1.1 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. +When a new element is added to the context, the functions +qmcklcontextcreate, qmcklcontextdestroy and qmcklcontextcopy +should be updated inorder to make deep copies.

    -
    -
    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;
    -
    -  /* Validity checking */
    -  uint32_t                   tag;
    -
    -} qmckl_context_struct;
    -
    -
    -

    A tag is used internally to check if the memory domain pointed @@ -496,7 +396,7 @@ if the context is valid, QMCKL_NULL_CONTEXT otherwise. 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) { @@ -510,8 +410,8 @@ if the context is valid, QMCKL_NULL_CONTEXT otherwise.

    -
    -

    1.2 Creation

    +
    +

    1.2 Creation

    To create a new context, qmckl_context_create() should be used. @@ -524,57 +424,69 @@ To create a new context, qmckl_context_create() should be used.

    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;
       }
     
    -  /* Set all pointers to NULL */
    -  memset(ctx, 0, sizeof(qmckl_context_struct));
    +  /* Set all pointers and values to NULL */
    +  {
    +    memset(ctx, 0, sizeof(qmckl_context_struct));
    +  }
     
       /* Initialize lock */
    -  init_lock(&(ctx->mutex));
    +  {
    +    pthread_mutexattr_t attr;
    +    int rc;
    +
    +    rc = pthread_mutexattr_init(&attr); 
    +    assert (rc == 0);
    +
    +    (void) pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
    +
    +    rc = pthread_mutex_init ( &(ctx->mutex), &attr);
    +    assert (rc == 0);
    +
    +    (void) pthread_mutexattr_destroy(&attr);
    +  }  
     
       /* 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;
    -}
    -
    -
    -
    -
    -
    -

    1.3 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. -

    + ctx->numprec.precision = QMCKL_DEFAULT_PRECISION; + ctx->numprec.range = QMCKL_DEFAULT_RANGE; -
    -
    qmckl_context qmckl_context_previous(const qmckl_context context) {
    +  ctx->ao_basis.uninitialized = (1 << 10) - 1;
     
    -  const qmckl_context checked_context = qmckl_context_check(context);
    -  if (checked_context == (qmckl_context) 0) {
    -    return (qmckl_context) 0;
    +  /* Allocate qmckl_memory_struct */
    +  {
    +    const size_t size = 128L;
    +    qmckl_memory_info_struct * new_array = calloc(size, sizeof(qmckl_memory_info_struct));
    +    if (new_array == NULL) {
    +      free(ctx);
    +      return QMCKL_NULL_CONTEXT; 
    +    }
    +    memset( &(new_array[0]), 0, size * sizeof(qmckl_memory_info_struct) );
    +
    +    ctx->memory.element = new_array;
    +    ctx->memory.array_size = size;
    +    ctx->memory.n_allocated = (size_t) 0;
       }
     
    -  const qmckl_context_struct* ctx = (qmckl_context_struct*) checked_context;
    -  return qmckl_context_check((qmckl_context) ctx->prev);
    +  return (qmckl_context) ctx;
     }
     
    -
    -

    1.4 Locking

    -
    +
    +

    1.3 Locking

    +

    For thread safety, the context may be locked/unlocked. The lock is initialized with the PTHREAD_MUTEX_RECURSIVE attribute, and the @@ -583,47 +495,32 @@ number of times the thread has locked it is saved in the

    -
    void init_lock(pthread_mutex_t* mutex) {
    -  pthread_mutexattr_t attr;
    -  int rc;
    -
    -  rc = pthread_mutexattr_init(&attr); 
    -  assert (rc == 0);
    -
    -  (void) pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
    -
    -  rc = pthread_mutex_init ( mutex, &attr);
    -  assert (rc == 0);
    -
    -  (void)pthread_mutexattr_destroy(&attr);
    -}
    -
    -void qmckl_lock(qmckl_context context) {
    +
    void qmckl_lock(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);
     */
    @@ -633,231 +530,99 @@ number of times the thread has locked it is saved in the
     
    -
    -

    1.5 Copy

    -
    +
    +

    1.4 TODO Copy

    +

    -qmckl_context_copy makes a shallow copy of a context. It returns +qmckl_context_copy makes a deep copy of a context. It returns QMCKL_NULL_CONTEXT upon failure.

    qmckl_context qmckl_context_copy(const qmckl_context context) {
     
    -  qmckl_lock(context);
    -
       const qmckl_context checked_context = qmckl_context_check(context);
     
       if (checked_context == QMCKL_NULL_CONTEXT) {
    -    qmckl_unlock(context);
         return QMCKL_NULL_CONTEXT;
       }
     
    +  /*
    +  qmckl_lock(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) malloc (context, sizeof(qmckl_context_struct));
     
    -  if (new_ctx == NULL) {
    -    qmckl_unlock(context);
    +    if (new_ctx == NULL) {
    +      qmckl_unlock(context);
    +      return QMCKL_NULL_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));
    +
    +    qmckl_unlock( (qmckl_context) new_ctx );
    +
    +    return (qmckl_context) new_ctx;
    +  }
    +  qmckl_unlock(context);
    +*/
         return QMCKL_NULL_CONTEXT;
    -  }
    -
    -  /* Copy the old context on the new one */
    -  memcpy(new_ctx, old_ctx, sizeof(qmckl_context_struct));
    -
    -  new_ctx->prev = old_ctx;
    -
    -  qmckl_unlock( (qmckl_context) new_ctx );
    -  qmckl_unlock( (qmckl_context) old_ctx );
    -
    -  return (qmckl_context) new_ctx;
     }
     
     
    -
    -

    1.6 Destroy

    -
    +
    +

    1.5 Destroy

    +

    The context is destroyed with qmckl_context_destroy, leaving the ancestors untouched. It frees the context, and returns the previous context.

    -
    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_context_struct* const ctx = (qmckl_context_struct* const) context;
    +  assert (ctx != NULL);  /* Shouldn't be possible because the context is valid */
     
       qmckl_lock(context);
    -
    -  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
    -  assert (ctx != NULL);  /* Shouldn't be true 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;
    +  {
    +    /* Memory: Remove all allocated data */
    +    for (size_t pos = (size_t) 0 ; pos < ctx->memory.array_size ; ++pos) {
    +      if (ctx->memory.element[pos].pointer != NULL) {
    +        free(ctx->memory.element[pos].pointer);
    +        memset( &(ctx->memory.element[pos]), 0, sizeof(qmckl_memory_info_struct) );
    +        ctx->memory.n_allocated -= 1;
    +      }
         }
    +    assert (ctx->memory.n_allocated == (size_t) 0);
    +    free(ctx->memory.element);
    +    ctx->memory.element = NULL;
    +    ctx->memory.array_size = (size_t) 0;
       }
    -
    -  qmckl_exit_code rc;
    -  rc = qmckl_context_remove_memory(context,ctx);
    -  assert (rc == QMCKL_SUCCESS);
    +  qmckl_unlock(context);
     
       ctx->tag = INVALID_TAG;
     
       const int rc_destroy = pthread_mutex_destroy( &(ctx->mutex) );
       if (rc_destroy != 0) {
    +/* DEBUG */
          fprintf(stderr, "qmckl_context_destroy: %s (count = %d)\n", strerror(rc_destroy), ctx->lock_count);
          abort();
       }
     
    -  rc = qmckl_free(context,ctx);
    -  assert (rc == QMCKL_SUCCESS);
    -
    -  return prev_context;
    -}
    -
    -
    -
    -
    -
    -
    -

    2 Memory allocation handling

    -
    -
    -
    -

    2.1 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. -

    - -
    -
    typedef struct qmckl_memory_struct {
    -  struct qmckl_memory_struct * next    ;
    -  void                       * pointer ;
    -  size_t                       size    ;
    -} qmckl_memory_struct;
    -
    -
    -
    -
    - -
    -

    2.2 Append memory

    -
    -

    -The following function, called in 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. -

    - -
    -
    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;
    -
    -}
    -
    -
    -
    -
    - -
    -

    2.3 Remove memory

    -
    -

    -The following function, called in 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. -

    - -
    -
    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);
    +  free(ctx);
     
       return QMCKL_SUCCESS;
     }
    @@ -866,982 +631,10 @@ immediately with QMCKL_SUCCESS.
     
    - -
    -

    3 Error handling

    -
    -
    -
    -

    3.1 Data structure

    -
    -
    -
    #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;
    -
    -
    -
    -
    - -
    -

    3.2 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. -

    - -
    -
    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;
    -}
    -
    -
    - -

    -The qmckl_context_set_error function returns a new context with -the error domain updated. -

    - -
    -
    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;
    -}
    -
    -
    - - -

    -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. -

    - -
    -
    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;
    -}
    -
    -
    -
    - -

    -For example, this function can be used as -

    -
    -
    if (x < 0) {
    -  return qmckl_failwith(context,
    -                        QMCKL_INVALID_ARG_2,
    -                        "qmckl_function", 
    -                        "Expected x >= 0");
    - }
    -
    -
    -
    -
    -
    - -
    -

    4 Control of the numerical precision

    -
    -

    -Controlling numerical precision enables optimizations. Here, the -default parameters determining the target numerical precision and -range are defined. -

    - - - - --- -- - - - - - - - - - - - -
    QMCKL_DEFAULT_PRECISION53
    QMCKL_DEFAULT_RANGE11
    - -
    -
    """ 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)
    -
    -
    -
    - -
    -
    typedef struct qmckl_precision_struct {
    -  int  precision;
    -  int  range;
    -} qmckl_precision_struct;
    -
    -
    - -

    -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. -

    -
    - -
    -

    4.1 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. -

    - -
    -
    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;
    -}
    -
    -
    - -
    -
    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
    -
    -
    - -

    -qmckl_context_set_precision returns a copy of the context with a -different precision parameter. -

    - -
    -
    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;
    -}
    -
    -
    - -

    -qmckl_context_get_precision returns the value of the numerical precision in the context. -

    - -
    -
    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;
    -}
    -
    -
    - -
    -
    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
    -
    -
    -
    -
    - -
    -

    4.2 Range

    -
    -

    -qmckl_context_update_range modifies the parameter for the numerical range in a given context. -

    - -
    -
    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;
    -}
    -
    -
    - -
    -
    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
    -
    -
    - -

    -qmckl_context_set_range returns a copy of the context with a different precision parameter. -

    - -
    -
    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;
    -}
    -
    -
    - -

    -qmckl_context_get_range returns the value of the numerical range in the context. -

    - -
    -
    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;
    -}
    -
    -
    -
    -
    -
    -

    4.3 Helper functions

    -
    -

    -qmckl_context_get_epsilon returns \(\epsilon = 2^{1-n}\) where n is the precision. -

    - -
    -
    double qmckl_context_get_epsilon(const qmckl_context context) {
    -  const int precision = qmckl_context_get_precision(context);
    -  return 1. /  (double) (1L << (precision-1));
    -}
    -
    -
    -
    -
    -
    -
    -

    5 TODO Basis set

    -
    -

    -For H2 with the following basis set, -

    - -
    -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
    -
    - -

    -we have: -

    - -
    -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]
    -
    -
    - -
    -

    5.1 Data structure

    -
    -
    -
    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;
    -
    -
    -
    -
    - -
    -

    5.2 qmckl_context_update_ao_basis

    -
    -

    -Updates the data describing the AO basis set into the context. -

    - - - - --- -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    typeGaussian or Slater
    shell_numNumber of shells
    prim_numTotal 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
    - -
    -
    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);
    -
    -
    -
    - -
    -

    5.2.1 Source

    -
    -
    -
    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 ; i<shell_num ; i++) {
    -    if (SHELL_CENTER[i] <= 0) return QMCKL_FAILURE;
    -    if (SHELL_PRIM_NUM[i] <= 0) return QMCKL_FAILURE;
    -    if (SHELL_ANG_MOM[i] < 0) return QMCKL_FAILURE;
    -    if (SHELL_PRIM_INDEX[i] < 0) return QMCKL_FAILURE;
    -  }
    -
    -  for (i=0 ; i<prim_num ; i++) {
    -    if (EXPONENT[i] <= 0) return QMCKL_FAILURE;
    -  }
    -
    -  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
    -  if (ctx == NULL) return QMCKL_FAILURE;
    -
    -  qmckl_ao_basis_struct* basis =
    -    (qmckl_ao_basis_struct*) qmckl_malloc (context, sizeof(qmckl_ao_basis_struct));
    -  if (basis == NULL) return QMCKL_ALLOCATION_FAILED;
    -
    -
    -  /* Memory allocations */
    -
    -  assert (basis->shell_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 ; i<shell_num ; i++) {
    -    basis->shell_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 ; i<prim_num ; i++) {
    -    basis->exponent   [i] = EXPONENT[i];
    -    basis->coefficient[i] = COEFFICIENT[i];
    -  }
    -
    -  ctx->ao_basis = basis;
    -  return QMCKL_SUCCESS;
    -}
    -
    -
    -
    -
    - -
    -

    5.2.2 Fortran interface

    -
    -
    -
    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
    -
    -
    -
    -
    - -
    -

    5.2.3 TODO Test

    -
    -
    - -
    -

    5.3 qmckl_context_set_ao_basis

    -
    -

    -Sets the data describing the AO basis set into the context. -

    - - - - --- -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    typeGaussian or Slater
    shell_numNumber of shells
    prim_numTotal 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
    - -
    -
    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);
    -
    -
    -
    - -
    -

    5.3.1 Source

    -
    -
    -
    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;
    -}
    -
    -
    -
    -
    - -
    -

    5.3.2 Fortran interface

    -
    -
    -
    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
    -
    -
    -
    -
    - -
    -

    5.3.3 TODO Test

    -
    -
    -

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/qmckl_distance.html b/qmckl_distance.html index 6f4d952..0be56d9 100644 --- a/qmckl_distance.html +++ b/qmckl_distance.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Inter-particle distances @@ -333,14 +333,14 @@ for the JavaScript code in this tag.

    Table of Contents

    -

    -Functions for the computation of distances between particles. -

    -
    -

    1 Squared distance

    +
    +

    1 Squared distance

    -
    -

    1.1 qmckl_distance_sq

    +
    +

    1.1 qmckl_distance_sq

    qmckl_distance_sq computes the matrix of the squared distances @@ -370,7 +367,7 @@ between all pairs of points in two sets, one point within each set: \]

    - +
    @@ -420,7 +417,7 @@ between all pairs of points in two sets, one point within each set: - + @@ -434,7 +431,7 @@ between all pairs of points in two sets, one point within each set: - + @@ -463,8 +460,8 @@ between all pairs of points in two sets, one point within each set:
    doubleA[3][lda]A[][lda] in Array containing the \(m \times 3\) matrix \(A\)
    doubleB[3][ldb]B[][ldb] in Array containing the \(n \times 3\) matrix \(B\)
    -
    -

    1.1.1 Requirements

    +
    +

    1.1.1 Requirements

    • context is not QMCKL_NULL_CONTEXT
    • @@ -482,8 +479,8 @@ between all pairs of points in two sets, one point within each set:
    -
    -

    1.1.2 C header

    +
    +

    1.1.2 C header

    qmckl_exit_code qmckl_distance_sq (
    @@ -496,21 +493,21 @@ between all pairs of points in two sets, one point within each set:
           const int64_t lda,
           const double* B,
           const int64_t ldb,
    -            double* C,
    +      double* const C,
           const int64_t ldc ); 
     
    -
    -

    1.1.3 Source

    +
    +

    1.1.3 Source

    integer function qmckl_distance_sq_f(context, transa, transb, m, n, A, LDA, B, LDB, C, LDC) result(info)
       use qmckl
       implicit none
    -  integer*8  , intent(in)  :: context
    +  integer(qmckl_context)  , intent(in)  :: context
       character  , intent(in)  :: transa, transb
       integer*8  , intent(in)  :: m, n
       integer*8  , intent(in)  :: lda
    @@ -637,8 +634,8 @@ between all pairs of points in two sets, one point within each set:
     
    -
    -

    1.1.4 Performance

    +
    +

    1.1.4 Performance

    This function might be more efficient when A and B are @@ -651,7 +648,7 @@ transposed.

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/qmckl_error.html b/qmckl_error.html index 6347d80..af6e739 100644 --- a/qmckl_error.html +++ b/qmckl_error.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Error handling @@ -311,14 +311,18 @@ for the JavaScript code in this tag.

    Table of Contents

    -
    -

    -
    +
    +

    +

    The library should never make the calling programs abort, nor perform any input/output operations. This decision has to be taken @@ -329,10 +333,11 @@ by the developer of the code calling the library. All the functions return with an exit code, defined as

    -
    typedef int32_t qmckl_exit_code;
    +
    typedef int32_t qmckl_exit_code;
     
    +

    The exit code returns the completion status of the function to the calling program. When a function call completed successfully, @@ -345,7 +350,7 @@ error code is returned to the program. Here is the complete list of exit codes.

    - +
    @@ -468,10 +473,20 @@ Here is the complete list of exit codes.

    -The qmckl_strerror converts an exit code into a string. The +The qmckl_string_of_error converts an exit code into a string. The string is assumed to be large enough to contain the error message (typically 128 characters).

    +
    +
    + +
    +

    1 Decoding errors

    +
    +

    +To decode the error messages, qmckl_string_of_error converts an +error code into a string. +

     128
    @@ -550,7 +565,8 @@ The text strings are extracted from the previous table.
     
    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(128)
        end subroutine qmckl_string_of_error
     end interface
    @@ -558,10 +574,143 @@ The text strings are extracted from the previous table.
     
    + +
    +

    2 Data structure in context

    +
    +

    +The strings are declared with a maximum fixed size to avoid +dynamic memory allocation. +

    + +
    +
    #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;
    +
    +
    +
    +
    + +
    +

    3 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. +

    + +
    +
    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;
    +}
    +
    +
    +
    +
    + +
    +

    4 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. If the +message is NULL, then the default message obtained by +qmckl_string_of_error is used. The return code of the function is +the desired return code. +Upon failure, a QMCKL_NULL_CONTEXT is returned. +

    + +
    +
    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 (strlen(function) < QMCKL_MAX_FUN_LEN);
    +  if (message != NULL) {
    +    assert (strlen(message)  < QMCKL_MAX_MSG_LEN);
    +  }
    +
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT)
    +    return QMCKL_INVALID_CONTEXT;
    +
    +  if (message == NULL) {
    +    qmckl_exit_code rc = 
    +      qmckl_set_error(context, exit_code, function, qmckl_string_of_error(exit_code));
    +    assert (rc == QMCKL_SUCCESS);
    +  } else {
    +    qmckl_exit_code rc = 
    +      qmckl_set_error(context, exit_code, function, message);
    +    assert (rc == QMCKL_SUCCESS);
    +  }
    +
    +  return exit_code;
    +}
    +
    +
    +
    + +

    +For example, this function can be used as +

    +
    +
    if (x < 0) {
    +  return qmckl_failwith(context,
    +                        QMCKL_INVALID_ARG_2,
    +                        "qmckl_function", 
    +                        "Expected x >= 0");
    + }
    +
    +
    +
    +

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/qmckl_memory.html b/qmckl_memory.html index 8daecb6..c946278 100644 --- a/qmckl_memory.html +++ b/qmckl_memory.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Memory management @@ -311,19 +311,69 @@ for the JavaScript code in this tag.

    Table of Contents

    + +
    +

    1 Memory data structure for the context

    +

    -We override the allocation functions to enable the possibility of -optimized libraries to fine-tune the memory allocation. +Every time a new block of memory is allocated, the information +relative to the allocation is stored in a new qmckl_memory_info_struct. +A qmckl_memory_info_struct contains the pointer to the memory block, +its size in bytes, and extra implementation-specific information such as +alignment, pinning, if the memory should be allocated on CPU or GPU +etc.

    +
    +
    typedef struct qmckl_memory_info_struct {
    +  size_t size;
    +  void*  pointer;
    +} qmckl_memory_info_struct;
     
    -
    -

    -
    +static const qmckl_memory_info_struct qmckl_memory_info_struct_zero = + { + .size = (size_t) 0, + .pointer = NULL + }; +
    +
    + +

    +The memory element of the context is a data structure which +contains an array of qmckl_memory_info_struct, the size of the +array, and the number of allocated blocks. +

    + +
    +
    typedef struct qmckl_memory_struct {
    +  size_t                    n_allocated;
    +  size_t                    array_size;
    +  qmckl_memory_info_struct* element;
    +} qmckl_memory_struct;
    +
    +
    +
    +
    + +
    +

    2 Passing info to allocation routines

    +
    +

    +Passing information to the allocation routine should be done by +passing an instance of a qmckl_memory_info_struct. +

    +
    +
    + +
    +

    3 Allocation/deallocation functions

    +

    Memory allocation inside the library should be done with qmckl_malloc. It lets the library choose how the memory will be @@ -337,51 +387,87 @@ If the allocation failed, the NULL pointer is returned.

    void* qmckl_malloc(qmckl_context context,
    -                   const size_t size);
    +                   const qmckl_memory_info_struct info);
     
    -

    -In this implementation, we use calloc because it initializes the -memory block to zero, so structs will have NULL-initialized pointers. -

    -
    -
    void* qmckl_malloc(qmckl_context context, const size_t size) {
    +
    void* qmckl_malloc(qmckl_context context, const qmckl_memory_info_struct info) {
     
    -  void * pointer = calloc(size, (size_t) 1);
    +  assert (qmckl_context_check(context) != QMCKL_NULL_CONTEXT);
     
    -  if (qmckl_context_check(context) != QMCKL_NULL_CONTEXT) {
    -    qmckl_exit_code rc;
    -    rc = qmckl_context_append_memory(context, pointer, size);
    -    assert (rc == QMCKL_SUCCESS);
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  /* Allocate memory and zero it */
    +  void * pointer = malloc(info.size);
    +  if (pointer == NULL) {
    +    return NULL;
       }
    +  memset(pointer, 0, info.size);
    +
    +  qmckl_lock(context);
    +  {
    +    /* If qmckl_memory_struct is full, reallocate a larger one */
    +    if (ctx->memory.n_allocated == ctx->memory.array_size) {
    +      const size_t old_size = ctx->memory.array_size;
    +      qmckl_memory_info_struct * new_array = reallocarray(ctx->memory.element,
    +                                                          2L * old_size,
    +                                                          sizeof(qmckl_memory_info_struct));
    +      if (new_array == NULL) {
    +        qmckl_unlock(context);
    +        free(pointer);
    +        return NULL;
    +      }
    +
    +      memset( &(new_array[old_size]), 0, old_size * sizeof(qmckl_memory_info_struct) );
    +      ctx->memory.element = new_array;
    +      ctx->memory.array_size = 2L * old_size;
    +    }
    +
    +    /* Find first NULL entry */
    +    size_t pos = (size_t) 0;
    +    while ( pos < ctx->memory.array_size && ctx->memory.element[pos].size > (size_t) 0) {
    +      pos += (size_t) 1;
    +    }
    +    assert (ctx->memory.element[pos].size == (size_t) 0);
    +
    +    /* Copy info at the new location */
    +    ctx->memory.element[pos].size    = info.size;
    +    ctx->memory.element[pos].pointer = pointer;
    +    ctx->memory.n_allocated += (size_t) 1;
    +  }
    +  qmckl_unlock(context);
     
       return pointer;
     }
     
    -
    -
    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
    -   end function qmckl_malloc
    -end interface
    -
    -
    -
    qmckl_context context = qmckl_context_create();    
    +
    /* Create a context */
    +qmckl_context context = qmckl_context_create();    
     
    -int *a = (int*) qmckl_malloc(context, 3*sizeof(int));
    +qmckl_memory_info_struct info = qmckl_memory_info_struct_zero;
    +info.size = (size_t) 3;
    +
    +/* Allocate an array of ints */
    +int *a = (int*) qmckl_malloc(context, info);
    +
    +/* Check that array of ints is OK */
     munit_assert(a != NULL);
    -
     a[0] = 1;  munit_assert_int(a[0], ==, 1);
     a[1] = 2;  munit_assert_int(a[1], ==, 2);
     a[2] = 3;  munit_assert_int(a[2], ==, 3);
    +
    +/* Allocate another array of ints */
    +int *b = (int*) qmckl_malloc(context, info);
    +
    +/* Check that array of ints is OK */
    +munit_assert(b != NULL);
    +b[0] = 1;  munit_assert_int(b[0], ==, 1);
    +b[1] = 2;  munit_assert_int(b[1], ==, 2);
    +b[2] = 3;  munit_assert_int(b[2], ==, 3);
     
    @@ -393,38 +479,53 @@ allocation and needs to be updated.
    qmckl_exit_code qmckl_free(qmckl_context context,
    -                           void *ptr);
    +                           void * const ptr);
     
    -
    interface
    -   integer (c_int32_t) function qmckl_free (context, ptr) bind(C)
    -     use, intrinsic :: iso_c_binding
    -     integer (c_int64_t), intent(in), value :: context
    -     type (c_ptr), intent(in), value :: ptr
    -   end function qmckl_free
    -end interface
    -
    -
    +
    qmckl_exit_code qmckl_free(qmckl_context context, void * const ptr) {
     
    -
    -
    qmckl_exit_code qmckl_free(qmckl_context context, void *ptr) {
    -  if (qmckl_context_check(context) != QMCKL_NULL_CONTEXT) {
    -
    -    if (ptr == NULL) {
    +  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
           return qmckl_failwith(context,
    -                            QMCKL_INVALID_ARG_2,
    +                            QMCKL_INVALID_CONTEXT,
                                 "qmckl_free",
    -                            "NULL pointer");
    +                            NULL);
    +  }
    +
    +  if (ptr == NULL) {
    +    return qmckl_failwith(context,
    +                          QMCKL_INVALID_ARG_2,
    +                          "qmckl_free",
    +                          "NULL pointer");
    +  }
    +
    +  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
    +
    +  qmckl_lock(context);
    +  {
    +    /* Find pointer in array of saved pointers */
    +    size_t pos = (size_t) 0;
    +    while ( pos < ctx->memory.array_size && ctx->memory.element[pos].pointer != ptr) {
    +      pos += (size_t) 1;
         }
     
    -    qmckl_exit_code rc;
    -    rc = qmckl_context_remove_memory(context, ptr);
    +    if (pos >= ctx->memory.array_size) {
    +      /* Not found */
    +      qmckl_unlock(context);
    +      return qmckl_failwith(context,
    +                            QMCKL_FAILURE,
    +                            "qmckl_free",
    +                            "Pointer not found in context");
    +    }
     
    -    assert (rc == QMCKL_SUCCESS);
    +    free(ptr);
    +
    +    memset( &(ctx->memory.element[pos]), 0, sizeof(qmckl_memory_info_struct) );
    +    ctx->memory.n_allocated -= (size_t) 1;
       }
    -  free(ptr);
    +  qmckl_unlock(context);
    +
       return QMCKL_SUCCESS;
     }
     
    @@ -434,7 +535,7 @@ allocation and needs to be updated.

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/qmckl_numprec.html b/qmckl_numprec.html new file mode 100644 index 0000000..e3af9fb --- /dev/null +++ b/qmckl_numprec.html @@ -0,0 +1,588 @@ + + + + + + + +Numerical precision + + + + + + + + + + + + + +
    + UP + | + HOME +
    +

    Numerical precision

    + + +
    +

    1 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. +

    + + + + +++ ++ + + + + + + + + + + + +
    QMCKL_DEFAULT_PRECISION53
    QMCKL_DEFAULT_RANGE11
    + +
    +
    typedef struct qmckl_numprec_struct {
    +  uint32_t  precision;
    +  uint32_t  range;
    +} qmckl_numprec_struct;
    +
    +
    + +

    +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. +

    +
    +
    + +
    +

    2 Precision

    +
    +

    +qmckl_context_set_numprec_precision modifies the parameter for the +numerical precision in the context. +

    + +
    +
    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;
    +}
    +
    +
    + +
    +
    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
    +
    +
    + +

    +qmckl_get_numprec_precision returns the value of the numerical precision in the context. +

    + +
    +
    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;
    +}
    +
    +
    + +
    +
    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
    +
    +
    +
    +
    + +
    +

    3 Range

    +
    +

    +qmckl_set_numprec_range modifies the parameter for the numerical +range in a given context. +

    + +
    +
    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;
    +}
    +
    +
    + +
    +
    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
    +
    +
    + +

    +qmckl_get_numprec_range returns the value of the numerical range in the context. +

    + +
    +
    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;
    +}
    +
    +
    +
    +
    +
    +

    4 Helper functions

    +
    +

    +qmckl_get_numprec_epsilon returns \(\epsilon = 2^{1-n}\) where n is the precision. +We need to remove the sign bit from the precision. +

    + +
    +
    double qmckl_get_numprec_epsilon(const qmckl_context context) {
    +  const int precision = qmckl_get_numprec_precision(context);
    +  return 1. /  (double) (1L << (precision-2));
    +}
    +
    +
    +
    +
    +
    +
    +

    Author: TREX CoE

    +

    Created: 2021-04-20 Tue 22:01

    +

    Validate

    +
    + + diff --git a/test_qmckl.html b/test_qmckl.html index 0b79223..4e4e9e0 100644 --- a/test_qmckl.html +++ b/test_qmckl.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Testing @@ -311,7 +311,7 @@ for the JavaScript code in this tag.

    Author: TREX CoE

    -

    Created: 2021-03-28 Sun 23:40

    +

    Created: 2021-04-20 Tue 22:01

    Validate

    diff --git a/theme.setup b/theme.setup index da0d396..bbb2ba6 100644 --- a/theme.setup +++ b/theme.setup @@ -13,3 +13,4 @@ #+LANGUAGE: en +