mirror of
https://github.com/TREX-CoE/qmckl.git
synced 2024-12-22 20:36:01 +01:00
370 lines
10 KiB
Org Mode
370 lines
10 KiB
Org Mode
#+TITLE: Numerical precision
|
|
#+SETUPFILE: ../tools/theme.setup
|
|
#+INCLUDE: ../tools/lib.org
|
|
|
|
* Headers :noexport:
|
|
|
|
#+begin_src c :tangle (eval c_test) :noweb yes
|
|
#include "qmckl.h"
|
|
#include "assert.h"
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
int main() {
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval h_private_type)
|
|
#ifndef QMCKL_NUMPREC_HPT
|
|
#define QMCKL_NUMPREC_HPT
|
|
|
|
#ifdef HAVE_STDINT_H
|
|
#include <stdint.h>
|
|
#elif HAVE_INTTYPES_H
|
|
#include <inttypes.h>
|
|
#endif
|
|
#+end_src
|
|
|
|
#+begin_src c :tangle (eval c)
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#ifdef HAVE_STDINT_H
|
|
#include <stdint.h>
|
|
#elif HAVE_INTTYPES_H
|
|
#include <inttypes.h>
|
|
#endif
|
|
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
#ifdef HAVE_FPE
|
|
#define _GNU_SOURCE
|
|
#include <fenv.h>
|
|
#include <signal.h>
|
|
#include <stdio.h>
|
|
#include <execinfo.h>
|
|
|
|
|
|
#define MAX_BACKTRACE_SIZE 100
|
|
|
|
void floatingPointExceptionHandler(int signal) {
|
|
void* backtraceArray[MAX_BACKTRACE_SIZE];
|
|
int backtraceSize = backtrace(backtraceArray, MAX_BACKTRACE_SIZE);
|
|
char** backtraceSymbols = backtrace_symbols(backtraceArray, backtraceSize);
|
|
|
|
// Print the backtrace
|
|
for (int i = 0; i < backtraceSize; ++i) {
|
|
printf("[%d] %s\n", i, backtraceSymbols[i]);
|
|
}
|
|
|
|
// Clean up the memory used by backtrace_symbols
|
|
free(backtraceSymbols);
|
|
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
|
|
static void __attribute__ ((constructor))
|
|
trapfpe ()
|
|
{
|
|
/* Enable some exceptions. At startup all exceptions are masked. */
|
|
|
|
feenableexcept (FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
|
|
signal(SIGFPE, floatingPointExceptionHandler);
|
|
}
|
|
#endif
|
|
|
|
#include "qmckl.h"
|
|
#include "qmckl_context_private_type.h"
|
|
#+end_src
|
|
|
|
* Control of the numerical precision
|
|
|
|
Controlling numerical precision enables optimizations. Here, the
|
|
default parameters determining the target numerical precision and
|
|
range are defined. Following the IEEE Standard for Floating-Point
|
|
Arithmetic (IEEE 754), /precision/ refers to the number of significand bits and /range/
|
|
refers to the number of exponent bits.
|
|
|
|
#+NAME: table-precision
|
|
| ~QMCKL_DEFAULT_PRECISION~ | 53 |
|
|
| ~QMCKL_DEFAULT_RANGE~ | 11 |
|
|
|
|
# We need to force Emacs not to indent the Python code:
|
|
# -*- org-src-preserve-indentation: t
|
|
|
|
#+begin_src python :var table=table-precision :results drawer :exports results
|
|
""" This script generates the C and Fortran constants from the org-mode table.
|
|
"""
|
|
|
|
result = [ "#+begin_src c :comments org :tangle (eval h_type)" ]
|
|
for (text, code) in table:
|
|
text=text.replace("~","")
|
|
result += [ f"#define {text:30s} {code:d}" ]
|
|
result += [ "#+end_src" ]
|
|
|
|
result += [ "" ]
|
|
|
|
result += [ "#+begin_src f90 :comments org :tangle (eval fh_func) :exports none" ]
|
|
for (text, code) in table:
|
|
text=text.replace("~","")
|
|
result += [ f" integer, parameter :: {text:30s} = {code:d}" ]
|
|
result += [ "#+end_src" ]
|
|
|
|
return '\n'.join(result)
|
|
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
:results:
|
|
#+begin_src c :comments org :tangle (eval h_type)
|
|
#define QMCKL_DEFAULT_PRECISION 53
|
|
#define QMCKL_DEFAULT_RANGE 11
|
|
#+end_src
|
|
|
|
#+begin_src f90 :comments org :tangle (eval fh_func) :exports none
|
|
integer, parameter :: QMCKL_DEFAULT_PRECISION = 53
|
|
integer, parameter :: QMCKL_DEFAULT_RANGE = 11
|
|
#+end_src
|
|
:end:
|
|
|
|
#+begin_src c :comments org :tangle (eval h_private_type)
|
|
typedef struct qmckl_numprec_struct {
|
|
uint32_t precision;
|
|
uint32_t range;
|
|
} qmckl_numprec_struct;
|
|
#+end_src
|
|
|
|
The following functions set and get the required precision and
|
|
range. ~precision~ is an integer between 2 and 53, and ~range~ is an
|
|
integer between 2 and 11.
|
|
|
|
The setter functions functions return a new context as a 64-bit
|
|
integer. The getter functions return the value, as a 32-bit
|
|
integer. The update functions return ~QMCKL_SUCCESS~ or ~QMCKL_FAILURE~.
|
|
|
|
* Precision ~qmckl_context_set_numprec_precision~ modifies the parameter for the
|
|
numerical precision in the context.
|
|
|
|
# Header
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const int precision);
|
|
#+end_src
|
|
|
|
# Source
|
|
#+begin_src c :tangle (eval c)
|
|
qmckl_exit_code qmckl_set_numprec_precision(const qmckl_context context, const int precision) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT)
|
|
return QMCKL_INVALID_CONTEXT;
|
|
|
|
if (precision < 2) {
|
|
return qmckl_failwith(context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_update_numprec_precision",
|
|
"precision < 2");
|
|
}
|
|
|
|
if (precision > 53) {
|
|
return qmckl_failwith(context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_update_numprec_precision",
|
|
"precision > 53");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
|
|
|
|
/* This should be always true because the context is valid */
|
|
assert (ctx != NULL);
|
|
|
|
qmckl_lock(context);
|
|
{
|
|
ctx->numprec.precision = (uint32_t) precision;
|
|
}
|
|
qmckl_unlock(context);
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
# Fortran interface
|
|
|
|
#+begin_src f90 :tangle (eval fh_func)
|
|
interface
|
|
integer (qmckl_exit_code) function qmckl_set_numprec_precision(context, precision) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
integer (qmckl_context), intent(in), value :: context
|
|
integer (c_int32_t), intent(in), value :: precision
|
|
end function qmckl_set_numprec_precision
|
|
end interface
|
|
#+end_src ~qmckl_get_numprec_precision~ returns the value of the numerical precision in the context.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
int32_t qmckl_get_numprec_precision(const qmckl_context context);
|
|
#+end_src
|
|
|
|
# Source
|
|
#+begin_src c :tangle (eval c)
|
|
int qmckl_get_numprec_precision(const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith(context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_numprec_precision",
|
|
"");
|
|
}
|
|
|
|
const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
|
|
return ctx->numprec.precision;
|
|
}
|
|
#+end_src
|
|
|
|
# Fortran interface
|
|
#+begin_src f90 :tangle (eval fh_func)
|
|
interface
|
|
integer (qmckl_exit_code) function qmckl_get_numprec_precision(context) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
integer (qmckl_context), intent(in), value :: context
|
|
end function qmckl_get_numprec_precision
|
|
end interface
|
|
#+end_src
|
|
|
|
* Range ~qmckl_set_numprec_range~ modifies the parameter for the numerical
|
|
range in a given context.
|
|
|
|
# Header
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int range);
|
|
#+end_src
|
|
|
|
# Source
|
|
#+begin_src c :tangle (eval c)
|
|
qmckl_exit_code qmckl_set_numprec_range(const qmckl_context context, const int range) {
|
|
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT)
|
|
return QMCKL_INVALID_CONTEXT;
|
|
|
|
if (range < 2) {
|
|
return qmckl_failwith(context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_numprec_range",
|
|
"range < 2");
|
|
}
|
|
|
|
if (range > 11) {
|
|
return qmckl_failwith(context,
|
|
QMCKL_INVALID_ARG_2,
|
|
"qmckl_set_numprec_range",
|
|
"range > 11");
|
|
}
|
|
|
|
qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
|
|
|
|
/* This should be always true because the context is valid */
|
|
assert (ctx != NULL);
|
|
|
|
qmckl_lock(context);
|
|
{
|
|
ctx->numprec.range = (uint32_t) range;
|
|
}
|
|
qmckl_unlock(context);
|
|
|
|
return QMCKL_SUCCESS;
|
|
}
|
|
#+end_src
|
|
|
|
# Fortran interface
|
|
#+begin_src f90 :tangle (eval fh_func)
|
|
interface
|
|
integer (qmckl_exit_code) function qmckl_set_numprec_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_set_numprec_range
|
|
end interface
|
|
#+end_src ~qmckl_get_numprec_range~ returns the value of the numerical range in the context.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
int32_t qmckl_get_numprec_range(const qmckl_context context);
|
|
#+end_src
|
|
|
|
# Source
|
|
#+begin_src c :tangle (eval c)
|
|
int qmckl_get_numprec_range(const qmckl_context context) {
|
|
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
|
|
return qmckl_failwith(context,
|
|
QMCKL_INVALID_CONTEXT,
|
|
"qmckl_get_numprec_range",
|
|
"");
|
|
}
|
|
|
|
const qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
|
|
return ctx->numprec.range;
|
|
}
|
|
#+end_src
|
|
|
|
# Fortran interface
|
|
#+begin_src f90 :tangle (eval fh_func) :exports none
|
|
interface
|
|
integer (qmckl_exit_code) function qmckl_get_numprec_range(context) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
integer (qmckl_context), intent(in), value :: context
|
|
end function qmckl_get_numprec_range
|
|
end interface
|
|
#+end_src
|
|
|
|
* Helper functions ~qmckl_get_numprec_epsilon~ returns $\epsilon = 2^{1-n}$ where ~n~ is the precision.
|
|
We need to remove the sign bit from the precision.
|
|
|
|
#+begin_src c :comments org :tangle (eval h_func) :exports none
|
|
double qmckl_get_numprec_epsilon(const qmckl_context context);
|
|
#+end_src
|
|
|
|
# Source
|
|
#+begin_src c :tangle (eval c)
|
|
double qmckl_get_numprec_epsilon(const qmckl_context context) {
|
|
const int precision = qmckl_get_numprec_precision(context);
|
|
return 1. / (double) (1L << (precision-2));
|
|
}
|
|
#+end_src
|
|
|
|
# Fortran interface
|
|
#+begin_src f90 :tangle (eval fh_func) :exports none
|
|
interface
|
|
real (c_double) function qmckl_get_numprec_epsilon(context) bind(C)
|
|
use, intrinsic :: iso_c_binding
|
|
import
|
|
integer (qmckl_context), intent(in), value :: context
|
|
end function qmckl_get_numprec_epsilon
|
|
end interface
|
|
#+end_src
|
|
|
|
* End of files :noexport:
|
|
|
|
#+begin_src c :comments link :tangle (eval h_private_type)
|
|
#endif
|
|
#+end_src
|
|
|
|
*** Test
|
|
#+begin_src c :comments link :tangle (eval c_test)
|
|
return 0;
|
|
}
|
|
#+end_src
|
|
|