qmckl/src/qmckl_context.org

#+TITLE: Context
#+SETUPFILE: ../docs/theme.setup

  This file is written in C because it is more natural to express the
  context in C than in Fortran.

  The context variable is a handle for the state of the library,
  and is stored in a data structure, which can't be seen outside of
  the library.  To simplify compatibility with other languages, the
  pointer to the internal data structure is converted into a 64-bit
  signed integer, defined in the ~qmckl_context~ type.  A value of ~0~
  for the context is equivalent to a ~NULL~ pointer.

   #+begin_src c :comments org :tangle (eval h)
typedef int64_t qmckl_context ;
   #+end_src

* Headers                                                          :noexport:

  #+NAME: filename
  #+begin_src elisp tangle: no 
(file-name-nondirectory (substring buffer-file-name 0 -4))
  #+end_src


  #+begin_src c :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "munit.h"
MunitResult test_<<filename()>>() {
  #+end_src

  #+begin_src c :tangle (eval h_private)
#ifndef __QMCKL_CONTEXT__
#define __QMCKL_CONTEXT__

#include <stdint.h>
#include "qmckl_error.h"

  #+end_src

  #+begin_src c :tangle (eval c) 
#include <stdint.h>
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "qmckl_error.h"
#include "qmckl_context_private.h"
#include "qmckl_context.h"
#include "qmckl_memory.h"
#include "qmckl_precision.h"

  #+end_src


* Context handling

  The tag is used internally to check if the memory domain pointed
  by a pointer is a valid <<<context>>>.

  #+begin_src c :comments org :tangle (eval h_private) :noweb yes
<<qmckl_error_struct>>
<<qmckl_ao_basis_struct>>

typedef struct qmckl_context_struct {

  struct qmckl_context_struct * prev;

 /* Molecular system */
  //  qmckl_nucleus_struct     * nucleus;
  //  qmckl_electron_struct    * electron;
  qmckl_ao_basis_struct      * ao_basis;
  //  qmckl_mo_struct          * mo;
  //  qmckl_determinant_struct * det;

 /* Numerical precision */
  uint32_t tag;
  int32_t precision;
  int32_t range;

 /* Error handling */
  qmckl_error_struct   * error;

} qmckl_context_struct;

#define VALID_TAG   0xBEEFFACE
#define INVALID_TAG 0xDEADBEEF
  #+end_src

** ~qmckl_context_check~

   Checks if the  domain pointed by the pointer is  a valid context.
   Returns  the input ~qmckl_context~ if  the context  is valid,  0
   otherwise.

   #+begin_src c :comments org :tangle (eval h)
qmckl_context qmckl_context_check(const qmckl_context context) ;
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context qmckl_context_check(const qmckl_context context) {

  if (context == (qmckl_context) 0) return (qmckl_context) 0;

  const qmckl_context_struct * ctx = (qmckl_context_struct*) context;

  if (ctx->tag != VALID_TAG) return (qmckl_context) 0;

  return context;
}
    #+end_src

** ~qmckl_context_create~

   To create a new context, use ~qmckl_context_create()~.
   - On success, returns a pointer to a context using the ~qmckl_context~ type
   - Returns ~0~ upon failure to allocate the internal data structure

     #+begin_src c :comments org :tangle (eval h)
qmckl_context qmckl_context_create();
     #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context qmckl_context_create() {

  qmckl_context_struct* context =
    (qmckl_context_struct*) qmckl_malloc ((qmckl_context) 0, sizeof(qmckl_context_struct));
  if (context == NULL) {
    return (qmckl_context) 0;
  }

  context->prev      = NULL;
  context->ao_basis  = NULL;
  context->precision = QMCKL_DEFAULT_PRECISION;
  context->range     = QMCKL_DEFAULT_RANGE;
  context->tag       = VALID_TAG;
  context->error     = NULL;

  return (qmckl_context) context;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int64_t) function qmckl_context_create() bind(C)
       use, intrinsic :: iso_c_binding
     end function qmckl_context_create
  end interface
    #+end_src

*** Test                                                           :noexport:
    #+begin_src c :comments link :tangle (eval c_test)
qmckl_context context = qmckl_context_create();
munit_assert_int64( context, !=, (qmckl_context) 0);
munit_assert_int64( qmckl_context_check(context), ==,  context);
    #+end_src

** ~qmckl_context_copy~

   This function makes a shallow copy of the current context.
   - Copying the 0-valued context returns 0
   - On success, returns a pointer to the new context using the ~qmckl_context~ type
   - Returns 0 upon failure to allocate the internal data structure
     for the new context

     #+begin_src c :comments org :tangle (eval h)
qmckl_context qmckl_context_copy(const qmckl_context context);
     #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context qmckl_context_copy(const qmckl_context context) {

  const qmckl_context checked_context = qmckl_context_check(context);

  if (checked_context == (qmckl_context) 0) {
    return (qmckl_context) 0;
  }

  qmckl_context_struct* old_context = (qmckl_context_struct*) checked_context;

  qmckl_context_struct* new_context =
    (qmckl_context_struct*) qmckl_malloc (context, sizeof(qmckl_context_struct));
  if (new_context == NULL) {
    return (qmckl_context) 0;
  }

  new_context->prev      = old_context;
  new_context->ao_basis  = old_context->ao_basis;
  new_context->precision = old_context->precision;
  new_context->range     = old_context->range;
  new_context->tag       = VALID_TAG;
  new_context->error     = old_context->error;

  return (qmckl_context) new_context;
}

    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int64_t) function qmckl_context_copy(context) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
     end function qmckl_context_copy
  end interface
    #+end_src

*** Test                                                           :noexport:
    #+begin_src c :comments link :tangle (eval c_test)
qmckl_context new_context = qmckl_context_copy(context);
munit_assert_int64(new_context, !=, (qmckl_context) 0);
munit_assert_int64(new_context, !=, context);
munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
    #+end_src

** ~qmckl_context_previous~

   Returns the previous context
   - On success, returns the ancestor of the current context
   - Returns 0 for the initial context
   - Returns 0 for the 0-valued context

     #+begin_src c :comments org :tangle (eval h)
qmckl_context qmckl_context_previous(const qmckl_context context);
     #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context qmckl_context_previous(const qmckl_context context) {

  const qmckl_context checked_context = qmckl_context_check(context);
  if (checked_context == (qmckl_context) 0) {
    return (qmckl_context) 0;
  }

  const qmckl_context_struct* ctx = (qmckl_context_struct*) checked_context;
  return qmckl_context_check((qmckl_context) ctx->prev);
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int64_t) function qmckl_context_previous(context) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
     end function qmckl_context_previous
  end interface
    #+end_src

*** Test                                                           :noexport:
    #+begin_src c :comments link :tangle (eval c_test)
munit_assert_int64(qmckl_context_previous(new_context), !=, (qmckl_context) 0);
munit_assert_int64(qmckl_context_previous(new_context), ==, context);
munit_assert_int64(qmckl_context_previous(context), ==, (qmckl_context) 0);
munit_assert_int64(qmckl_context_previous((qmckl_context) 0), ==, (qmckl_context) 0);
    #+end_src

** ~qmckl_context_destroy~

   Destroys the current context, leaving the ancestors untouched.
   - Succeeds if the current context is properly destroyed
   - Fails otherwise
   - Fails if the 0-valued context is given in argument
   - Fails if the the pointer is not a valid context

   #+begin_src c :comments org :tangle (eval h)
qmckl_exit_code qmckl_context_destroy(qmckl_context context);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_exit_code qmckl_context_destroy(const qmckl_context context) {

  const qmckl_context checked_context = qmckl_context_check(context);
  if (checked_context == (qmckl_context) 0) return QMCKL_FAILURE;

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  ctx->tag = INVALID_TAG;
  return qmckl_free(context,ctx);
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int32_t) function qmckl_context_destroy(context) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
     end function qmckl_context_destroy
  end interface
    #+end_src

*** Test                                                           :noexport:
    #+begin_src c :tangle (eval c_test)
munit_assert_int64(qmckl_context_check(new_context), ==, new_context);
munit_assert_int64(new_context, !=, (qmckl_context) 0);
munit_assert_int32(qmckl_context_destroy(new_context), ==, QMCKL_SUCCESS);
munit_assert_int64(qmckl_context_check(new_context), !=, new_context);
munit_assert_int64(qmckl_context_check(new_context), ==, (qmckl_context) 0);
munit_assert_int64(qmckl_context_destroy((qmckl_context) 0), ==, QMCKL_FAILURE);
    #+end_src

* Error handling
** Data structure
  
   #+NAME: qmckl_error_struct
   #+begin_src c :comments org 
#define  QMCKL_MAX_FUN_LEN  256
#define  QMCKL_MAX_MSG_LEN  1024

typedef struct qmckl_error_struct {

  qmckl_exit_code exit_code;
  char function[QMCKL_MAX_FUN_LEN];
  char message [QMCKL_MAX_MSG_LEN];

} qmckl_error_struct;
   #+end_src

** ~qmckl_context_update_error~

   #+begin_src c :comments org :tangle (eval h)
qmckl_exit_code
qmckl_context_update_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_exit_code
qmckl_context_update_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message)
{
  assert (context != 0);
  assert (function != NULL);
  assert (message != NULL);
  assert (exit_code > 0);
  assert (exit_code < QMCKL_INVALID_EXIT_CODE);

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  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;
  strcpy(error->function, function);
  strcpy(error->message, message);

  ctx->error = error;

  return QMCKL_SUCCESS;
}
    #+end_src

*** TODO Test

** ~qmckl_context_set_error~

   #+begin_src c :comments org :tangle (eval h)
qmckl_context
qmckl_context_set_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context
qmckl_context_set_error(qmckl_context context, const qmckl_exit_code exit_code, const char* function, const char* message)
{
  assert (context != 0);
  assert (function != NULL);
  assert (message != NULL);
  assert (exit_code > 0);
  assert (exit_code < QMCKL_INVALID_EXIT_CODE);

  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return context;

  if (qmckl_context_update_error(new_context, exit_code,
                                 function, message) != QMCKL_SUCCESS) {
    return context;
  }

  return new_context;
}
    #+end_src

*** TODO Test

*** Test                                                           :noexport:

** ~qmckl_failwith~

   To make a function fail, the <<<~qmckl_failwith~>>> function should be
   called, such that information about the failure is stored in
   the context. The desired exit code is given as an argument, as
   well as the name of the function and an error message. The return
   code of the function is the desired return code.

   #+begin_src c :comments org :tangle (eval h)
qmckl_exit_code qmckl_failwith(qmckl_context context,
                               const qmckl_exit_code exit_code,
                               const char* function,
                               const char* message) ;
   #+end_src

   #+begin_src c :comments org :tangle (eval c)
qmckl_exit_code qmckl_failwith(qmckl_context context,
                               const qmckl_exit_code exit_code,
                               const char* function,
                               const char* message) {

  if (context == 0) return QMCKL_NULL_CONTEXT;
  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);
  
  context = qmckl_context_set_error(context, exit_code, function, message);
  return exit_code;
}

   #+end_src

   For example, this function can be used as
   #+begin_src c :tangle no
if (x < 0) {
  return qmckl_failwith(context,
                        QMCKL_INVALID_ARG_2,
                        "qmckl_function", 
                        "Expected x >= 0");
 }
   #+end_src
    
   #   To decode the error messages, ~qmckl_strerror~ converts an
   #   error code into a string. 
* Basis set

   For H_2 with the following basis set,

   #+BEGIN_EXAMPLE
HYDROGEN
S   5
1         3.387000E+01           6.068000E-03
2         5.095000E+00           4.530800E-02
3         1.159000E+00           2.028220E-01
4         3.258000E-01           5.039030E-01
5         1.027000E-01           3.834210E-01
S   1
1         3.258000E-01           1.000000E+00
S   1
1         1.027000E-01           1.000000E+00
P   1
1         1.407000E+00           1.000000E+00
P   1
1         3.880000E-01           1.000000E+00
D   1
1         1.057000E+00           1.0000000
   #+END_EXAMPLE

   we have:

   #+BEGIN_EXAMPLE
type = 'G'
shell_num = 12
prim_num = 20
SHELL_CENTER = [1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
SHELL_ANG_MOM = ['S', 'S', 'S', 'P', 'P', 'D', 'S', 'S', 'S', 'P', 'P', 'D']
SHELL_PRIM_NUM = [5, 1, 1, 1, 1, 1, 5, 1, 1, 1, 1, 1]
prim_index = [1, 6, 7, 8, 9, 10, 11, 16, 17, 18, 19, 20]
EXPONENT = [ 33.87, 5.095, 1.159, 0.3258, 0.1027, 0.3258, 0.1027,
             1.407, 0.388, 1.057, 33.87, 5.095, 1.159, 0.3258, 0.1027,
             0.3258, 0.1027, 1.407, 0.388, 1.057]
COEFFICIENT = [ 0.006068, 0.045308, 0.202822, 0.503903, 0.383421,
                1.0, 1.0, 1.0, 1.0, 1.0, 0.006068, 0.045308, 0.202822,
                0.503903, 0.383421, 1.0, 1.0, 1.0, 1.0, 1.0]
   #+END_EXAMPLE

** Data structure

   #+NAME: qmckl_ao_basis_struct
   #+begin_src c :comments org :tangle no
typedef struct qmckl_ao_basis_struct {

  int64_t   shell_num;
  int64_t   prim_num;
  int64_t * shell_center;
  int32_t * shell_ang_mom;
  double  * shell_factor;
  double  * exponent    ;
  double  * coefficient ;
  int64_t * shell_prim_num;
  char      type;

} qmckl_ao_basis_struct;
   #+end_src

** ~qmckl_context_update_ao_basis~

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

   | ~type~                        | Gaussian or Slater                                                   |
   | ~shell_num~                   | Number of shells                                                     |
   | ~prim_num~                    | Total number of primitives                                           |
   | ~SHELL_CENTER(shell_num)~     | Id of the nucleus on which the shell is centered                     |
   | ~SHELL_ANG_MOM(shell_num)~    | Id of the nucleus on which the shell is centered                     |
   | ~SHELL_FACTOR(shell_num)~     | Normalization factor for the shell                                   |
   | ~SHELL_PRIM_NUM(shell_num)~   | Number of primitives in the shell                                    |
   | ~SHELL_PRIM_INDEX(shell_num)~ | Address of the first primitive of the shelll in the ~EXPONENT~ array |
   | ~EXPONENT(prim_num)~          | Array of exponents                                                   |
   | ~COEFFICIENT(prim_num)~       | Array of coefficients                                                |

   #+begin_src c :comments org :tangle (eval h)
qmckl_exit_code
qmckl_context_update_ao_basis(qmckl_context   context     , const char      type,
                              const int64_t   shell_num   , const int64_t   prim_num,
                              const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                              const double  * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                              const int64_t * SHELL_PRIM_INDEX,
                              const double  * EXPONENT    , const double  * COEFFICIENT);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_exit_code
qmckl_context_update_ao_basis(qmckl_context   context     , const char      type,
                              const int64_t   shell_num   , const int64_t   prim_num,
                              const int64_t * SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                              const double  * SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                              const int64_t * SHELL_PRIM_INDEX,
                              const double  * EXPONENT    , const double  * COEFFICIENT)
{

  int64_t i;

 /* Check input */

  if (type != 'G' && type != 'S') return QMCKL_FAILURE;
  if (shell_num <= 0) return QMCKL_FAILURE;
  if (prim_num  <= 0) return QMCKL_FAILURE;
  if (prim_num  < shell_num) return QMCKL_FAILURE;

  for (i=0 ; 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*) malloc (sizeof(qmckl_ao_basis_struct));
  if (basis == NULL) return QMCKL_FAILURE;


 /* Memory allocations */

  basis->shell_center  = (int64_t*) malloc (shell_num * sizeof(int64_t));
  if (basis->shell_center == NULL) {
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->shell_ang_mom = (int32_t*) malloc (shell_num * sizeof(int32_t));
  if (basis->shell_ang_mom == NULL) {
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->shell_prim_num= (int64_t*) malloc (shell_num * sizeof(int64_t));
  if (basis->shell_prim_num == NULL) {
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->shell_factor  = (double *) malloc (shell_num * sizeof(double ));
  if (basis->shell_factor == NULL) {
    qmckl_free(context, basis->shell_prim_num);
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->exponent      = (double *) malloc (prim_num  * sizeof(double ));
  if (basis->exponent == NULL) {
    qmckl_free(context, basis->shell_factor);
    qmckl_free(context, basis->shell_prim_num);
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    return QMCKL_FAILURE;
  }

  basis->coefficient   = (double *) malloc (prim_num  * sizeof(double ));
  if (basis->coefficient == NULL) {
    qmckl_free(context, basis->exponent);
    qmckl_free(context, basis->shell_factor);
    qmckl_free(context, basis->shell_prim_num);
    qmckl_free(context, basis->shell_ang_mom);
    qmckl_free(context, basis->shell_center);
    qmckl_free(context, basis);
    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;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int32_t) function qmckl_context_update_ao_basis(context, &
          typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
          SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
       character(c_char)  , intent(in), value :: typ
       integer (c_int64_t), intent(in), value :: shell_num
       integer (c_int64_t), intent(in), value :: prim_num
       integer (c_int64_t), intent(in)        :: SHELL_CENTER(shell_num)
       integer (c_int32_t), intent(in)        :: SHELL_ANG_MOM(shell_num)
       double precision   , intent(in)        :: SHELL_FACTOR(shell_num)
       integer (c_int64_t), intent(in)        :: SHELL_PRIM_NUM(shell_num)
       integer (c_int64_t), intent(in)        :: SHELL_PRIM_INDEX(shell_num)
       double precision   , intent(in)        :: EXPONENT(prim_num)
       double precision   , intent(in)        :: COEFFICIENT(prim_num)
     end function qmckl_context_update_ao_basis
  end interface
    #+end_src

*** TODO Test

** ~qmckl_context_set_ao_basis~

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

   | ~type~                        | Gaussian or Slater                                                   |
   | ~shell_num~                   | Number of shells                                                     |
   | ~prim_num~                    | Total number of primitives                                           |
   | ~SHELL_CENTER(shell_num)~     | Id of the nucleus on which the shell is centered                     |
   | ~SHELL_ANG_MOM(shell_num)~    | Id of the nucleus on which the shell is centered                     |
   | ~SHELL_FACTOR(shell_num)~     | Normalization factor for the shell                                   |
   | ~SHELL_PRIM_NUM(shell_num)~   | Number of primitives in the shell                                    |
   | ~SHELL_PRIM_INDEX(shell_num)~ | Address of the first primitive of the shelll in the ~EXPONENT~ array |
   | ~EXPONENT(prim_num)~          | Array of exponents                                                   |
   | ~COEFFICIENT(prim_num)~       | Array of coefficients                                                |

   #+begin_src c :comments org :tangle (eval h)
qmckl_context
qmckl_context_set_ao_basis(const qmckl_context context     , const char      type,
                           const int64_t       shell_num   , const int64_t   prim_num,
                           const int64_t *     SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                           const double  *     SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                           const int64_t *     SHELL_PRIM_INDEX,
                           const double  *     EXPONENT    , const double  * COEFFICIENT);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context
qmckl_context_set_ao_basis(const qmckl_context context     , const char      type,
                           const int64_t       shell_num   , const int64_t   prim_num,
                           const int64_t *     SHELL_CENTER, const int32_t * SHELL_ANG_MOM,
                           const double  *     SHELL_FACTOR, const int64_t * SHELL_PRIM_NUM,
                           const int64_t *     SHELL_PRIM_INDEX,
                           const double  *     EXPONENT    , const double  * COEFFICIENT)
{

  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return 0;

  if (qmckl_context_update_ao_basis(new_context, type, shell_num, prim_num,
                                    SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR,
                                    SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT,
                                    COEFFICIENT
                                    ) == QMCKL_FAILURE)
    return 0;

  return new_context;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int64_t) function qmckl_context_set_ao_basis(context, &
          typ, shell_num, prim_num, SHELL_CENTER, SHELL_ANG_MOM, SHELL_FACTOR, &
          SHELL_PRIM_NUM, SHELL_PRIM_INDEX, EXPONENT, COEFFICIENT) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
       character(c_char)  , intent(in), value :: typ
       integer (c_int64_t), intent(in), value :: shell_num
       integer (c_int64_t), intent(in), value :: prim_num
       integer (c_int64_t), intent(in)        :: SHELL_CENTER(shell_num)
       integer (c_int32_t), intent(in)        :: SHELL_ANG_MOM(shell_num)
       double precision   , intent(in)        :: SHELL_FACTOR(shell_num)
       integer (c_int64_t), intent(in)        :: SHELL_PRIM_NUM(shell_num)
       integer (c_int64_t), intent(in)        :: SHELL_PRIM_INDEX(shell_num)
       double precision   , intent(in)        :: EXPONENT(prim_num)
       double precision   , intent(in)        :: COEFFICIENT(prim_num)
     end function qmckl_context_set_ao_basis
  end interface
    #+end_src

*** TODO Test

* Precision

  The  following  functions  set   and  get  the  expected  required
  precision and  range. ~precision~ should  be an integer  between 2
  and 53, and ~range~ should be 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~.

** ~qmckl_context_update_precision~
   Modifies the parameter for the numerical precision in a given context.
   #+begin_src c :comments org :tangle (eval h)
qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_exit_code qmckl_context_update_precision(const qmckl_context context, const int precision) {

  if (precision <  2) return QMCKL_FAILURE;
  if (precision > 53) return QMCKL_FAILURE;

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  ctx->precision = precision;
  return QMCKL_SUCCESS;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int32_t) function qmckl_context_update_precision(context, precision) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
       integer (c_int32_t), intent(in), value :: precision
     end function qmckl_context_update_precision
  end interface
    #+end_src

*** TODO Tests                                                     :noexport:
** ~qmckl_context_update_range~
   Modifies the parameter for the numerical range in a given context.
   #+begin_src c :comments org :tangle (eval h)
qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range);
    #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_exit_code qmckl_context_update_range(const qmckl_context context, const int range) {

  if (range <  2) return QMCKL_FAILURE;
  if (range > 11) return QMCKL_FAILURE;

  qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  if (ctx == NULL) return QMCKL_FAILURE;

  ctx->range = range;
  return QMCKL_SUCCESS;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int32_t) function qmckl_context_update_range(context, range) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
       integer (c_int32_t), intent(in), value :: range
     end function qmckl_context_update_range
  end interface
    #+end_src

*** TODO Tests                                                     :noexport:
** ~qmckl_context_set_precision~
   Returns a copy of the context with a different precision parameter.
   #+begin_src c :comments org :tangle (eval h)
qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context qmckl_context_set_precision(const qmckl_context context, const int precision) {
  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return 0;

  if (qmckl_context_update_precision(new_context, precision) == QMCKL_FAILURE) return 0;

  return new_context;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int64_t) function qmckl_context_set_precision(context, precision) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
       integer (c_int32_t), intent(in), value :: precision
     end function qmckl_context_set_precision
  end interface
    #+end_src

*** TODO Tests                                                     :noexport:
** ~qmckl_context_set_range~
   Returns a copy of the context with a different precision parameter.
   #+begin_src c :comments org :tangle (eval h)
qmckl_context qmckl_context_set_range(const qmckl_context context, const int range);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
qmckl_context qmckl_context_set_range(const qmckl_context context, const int range) {
  qmckl_context new_context = qmckl_context_copy(context);
  if (new_context == 0) return 0;

  if (qmckl_context_update_range(new_context, range) == QMCKL_FAILURE) return 0;

  return new_context;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int64_t) function qmckl_context_set_range(context, range) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
       integer (c_int32_t), intent(in), value :: range
     end function qmckl_context_set_range
  end interface
    #+end_src

*** TODO Tests                                                     :noexport:

** ~qmckl_context_get_precision~
   Returns the value of the numerical precision in the context
   #+begin_src c :comments org :tangle (eval h)
int32_t qmckl_context_get_precision(const qmckl_context context);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
int qmckl_context_get_precision(const qmckl_context context) {
  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  return ctx->precision;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int32_t) function qmckl_context_get_precision(context) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
     end function qmckl_context_get_precision
  end interface
    #+end_src

*** TODO Tests                                                     :noexport:
** ~qmckl_context_get_range~
   Returns the value of the numerical range in the context
   #+begin_src c :comments org :tangle (eval h)
int32_t qmckl_context_get_range(const qmckl_context context);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
int qmckl_context_get_range(const qmckl_context context) {
  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  return ctx->range;
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     integer (c_int32_t) function qmckl_context_get_range(context) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
     end function qmckl_context_get_range
  end interface
    #+end_src

*** TODO Tests                                                     :noexport:

** ~qmckl_context_get_epsilon~
   Returns $\epsilon = 2^{1-n}$ where ~n~ is the precision
   #+begin_src c :comments org :tangle (eval h)
double qmckl_context_get_epsilon(const qmckl_context context);
   #+end_src

*** Source
    #+begin_src c :tangle (eval c)
double qmckl_context_get_epsilon(const qmckl_context context) {
  const qmckl_context_struct* ctx = (qmckl_context_struct*) context;
  return pow(2.0,(double) 1-ctx->precision);
}
    #+end_src

*** Fortran interface
    #+begin_src f90 :tangle (eval fh)
  interface
     real (c_double) function qmckl_context_get_epsilon(context) bind(C)
       use, intrinsic :: iso_c_binding
       integer (c_int64_t), intent(in), value :: context
     end function qmckl_context_get_epsilon
  end interface
    #+end_src

*** TODO Tests                                                     :noexport:

* End of files                                                     :noexport:

  #+begin_src c :comments link :tangle (eval h_private)

#endif
  #+end_src
  
*** Test
     #+begin_src c :comments link :tangle (eval c_test)
return MUNIT_OK;
}
     #+end_src



# -*- mode: org -*-
# vim: syntax=c