core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'develop' of github.com:ecmwf/eccodes into develop

This commit is contained in:
Shahram Najm 2023-10-13 15:35:14 +01:00
parent 1c24fa08f5 7961163402
commit a14bc85fc0
6 changed files with 147 additions and 238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/eccodes_prototypes.h
View File

@ -765,8 +765,6 @@ unsigned long grib_ibm_nearest_smaller_to_long(double x);
int grib_nearest_smaller_ibm_float(double a, double* ret);
/* grib_ieeefloat.cc*/
double grib_ieee_table_e(unsigned long e);
double grib_ieee_table_v(unsigned long e);
unsigned long grib_ieee_to_long(double x);
double grib_ieeefloat_error(double x);
double grib_long_to_ieee(unsigned long x);

102
src/grib_ibmfloat.cc
View File

@ -9,91 +9,9 @@
*/
#include "grib_api_internal.h"
#include "grib_ibmfloat.h"
#if GRIB_PTHREADS
static pthread_once_t once = PTHREAD_ONCE_INIT;
static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
static void init()
{
pthread_mutexattr_t attr;
pthread_mutexattr_init(&attr);
pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
pthread_mutex_init(&mutex, &attr);
pthread_mutexattr_destroy(&attr);
}
#elif GRIB_OMP_THREADS
static int once = 0;
static omp_nest_lock_t mutex;
static void init()
{
GRIB_OMP_CRITICAL(lock_grib_ibmfloat_c)
{
if (once == 0) {
omp_init_nest_lock(&mutex);
once = 1;
}
}
}
#endif
typedef struct ibm_table_t ibm_table_t;
struct ibm_table_t
{
int inited;
double e[128];
double v[128];
double vmin;
double vmax;
};
static ibm_table_t ibm_table = { 0, {0,}, {0,}, 0, 0 };
/**
.. seealso::
Documentation for ``init_ieee_float()`` in ``grib_ieeefloat.cc``
*/
static void init_ibm_table()
{
if (!ibm_table.inited) {
unsigned long i;
unsigned long mmin = 0x100000;
unsigned long mmax = 0xffffff;
double e = 1;
for (i = 1; i <= 57; i++) {
e *= 16;
ibm_table.e[i + 70] = e;
ibm_table.v[i + 70] = e * mmin;
}
ibm_table.e[70] = 1;
ibm_table.v[70] = mmin;
e = 1;
for (i = 1; i <= 70; i++) {
e /= 16;
ibm_table.e[70 - i] = e;
ibm_table.v[70 - i] = e * mmin;
}
ibm_table.vmin = ibm_table.v[0];
ibm_table.vmax = ibm_table.e[127] * mmax;
ibm_table.inited = 1;
/*for (i=0;i<128;i++) printf("++++ ibm_table.v[%d]=%g\n",i,ibm_table.v[i]);*/
}
}
static void init_table_if_needed()
{
GRIB_MUTEX_INIT_ONCE(&once, &init);
GRIB_MUTEX_LOCK(&mutex);
if (!ibm_table.inited)
init_ibm_table();
GRIB_MUTEX_UNLOCK(&mutex);
}
constexpr auto ibm_table = IbmTable{};
static void binary_search(const double xx[], const unsigned long n, double x, unsigned long* j)
{
@ -122,8 +40,6 @@ unsigned long grib_ibm_to_long(double x)
unsigned long e = 0;
double rmmax = mmax + 0.5;
init_table_if_needed();
/* printf("\ngrib_ibm_to_long: x=%.20e\n",x); */
if (x < 0) {
s = 1;
@ -143,7 +59,7 @@ unsigned long grib_ibm_to_long(double x)
return 0;
}
binary_search(ibm_table.v, 127, x, &e);
binary_search(ibm_table.v.data(), 127, x, &e);
/* printf("grib_ibm_to_long: e=%ld\n",e); */
@ -180,8 +96,6 @@ double grib_ibmfloat_error(double x)
{
unsigned long e = 0;
init_table_if_needed();
if (x < 0)
x = -x;
@ -196,7 +110,7 @@ double grib_ibmfloat_error(double x)
return 0;
}
binary_search(ibm_table.v, 127, x, &e);
binary_search(ibm_table.v.data(), 127, x, &e);
return ibm_table.e[e];
}
@ -209,8 +123,6 @@ double grib_long_to_ibm(unsigned long x)
double val = m;
init_table_if_needed();
/*if(x == 0) return 0;*/
if (c == 0 && m <= 1)
return 0;
@ -225,13 +137,11 @@ double grib_long_to_ibm(unsigned long x)
double grib_ibm_table_e(unsigned long e)
{
init_table_if_needed();
return ibm_table.e[e];
}
double grib_ibm_table_v(unsigned long e)
{
init_table_if_needed();
return ibm_table.v[e];
}
@ -247,8 +157,6 @@ unsigned long grib_ibm_nearest_smaller_to_long(double x)
if (x == 0)
return 0;
init_table_if_needed();
l = grib_ibm_to_long(x);
y = grib_long_to_ibm(l);
@ -293,8 +201,6 @@ int grib_nearest_smaller_ibm_float(double a, double* ret)
{
unsigned long l = 0;
init_table_if_needed();
if (a > ibm_table.vmax)
return GRIB_INTERNAL_ERROR;

74
src/grib_ibmfloat.h Normal file
View File

@ -0,0 +1,74 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#pragma once
#include "grib_scaling.h"
#include <array>
#include <type_traits>
#include <cstdint>
/**
.. _init_ieee_table:
Init IBM Floats Table
=====================
Initializes the ibm_table with IBM Float values. Nearest smaller values (e.g., reference values for grid_simple) are taken from this table.
Details
-------
The table layout is as follows:
+-------+----------------+------------------------+
| idx (i) | multiplier (e) | value (v = mmin * e) |
+-------+----------------+------------------------+
| 0 | 16^(-70) | 0x100000 * 2^(-70) |
| 1 | 16^(-69) | 0x100000 * 2^(-69) |
| ... | ... | ... |
| 126 | 16^56 | 0x100000 * 2^56 |
| 127 | 16^57 | 0x100000 * 2^57 |
+-------+----------------+------------------------+
The vmin and vmax boundaries are defined as:
- vmin = 0x100000 * 2^(-70)
- vmax = 0xffffff * 2^57
*/
struct IbmTable {
private:
using ValueType = double;
static constexpr int TABLESIZE = 128;
static constexpr uint32_t mantissa_min = 0x100000;
static constexpr uint32_t mantissa_max = 0xffffff;
public:
static constexpr std::array<ValueType, TABLESIZE> e = []() {
std::array<ValueType, TABLESIZE> multiplier{};
multiplier[0] = 0;
for (int i = 0; i < TABLESIZE; ++i) {
multiplier[i] = codes_power<ValueType>(i - 70, 16);
}
return multiplier;
}();
static constexpr std::array<ValueType, TABLESIZE> v = []() {
std::array<ValueType, TABLESIZE> values{};
values[0] = 0;
for (int i = 0; i < TABLESIZE; ++i) {
values[i] = e[i] * mantissa_min;
}
return values;
}();
static constexpr ValueType vmin = e[0] * mantissa_min;
static constexpr ValueType vmax = e[127] * mantissa_max;
};

139
src/grib_ieeefloat.cc
View File

@ -8,121 +8,11 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
/***************************************************************************
* Enrico Fucile - 06.01.2009 *
***************************************************************************/
#include "grib_ieeefloat.h"
#if GRIB_PTHREADS
static pthread_once_t once = PTHREAD_ONCE_INIT;
static pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
static void init()
{
pthread_mutexattr_t attr;
pthread_mutexattr_init(&attr);
pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
pthread_mutex_init(&mutex, &attr);
pthread_mutexattr_destroy(&attr);
}
#elif GRIB_OMP_THREADS
static int once = 0;
static omp_nest_lock_t mutex;
static void init()
{
GRIB_OMP_CRITICAL(lock_grib_ieeefloat_c)
{
if (once == 0) {
omp_init_nest_lock(&mutex);
once = 1;
}
}
}
#endif
/* See old implementation in src/deprecated/grib_ieeefloat.c */
typedef struct ieee_table_t ieee_table_t;
struct ieee_table_t
{
int inited;
double e[255];
double v[255];
double vmin;
double vmax;
};
static ieee_table_t ieee_table = { 0, {0,}, {0,}, 0, 0 };
/**
.. _init_ieee_table:
Init IEEE Table
===============
Initializes the ieee_table with IEEE754 single precision (32-bit) values. Nearest smaller values (e.g., reference values for grid_simple and grid_ccsds) are taken from this table.
Details
-------
The table layout is as follows:
+-------+----------------+----------------------+
| idx (i) | multiplier (e) | value (v = mmin * e) |
+-------+----------------+----------------------+
| 1 | 2^(-149) | 0x800000 * 2^(-149) |
| 2 | 2^(-148) | 0x800000 * 2^(-148) |
| ... | ... | ... |
| 253 | 2^103 | 0x800000 * 2^103 |
| 254 | 2^104 | 0x800000 * 2^104 |
+-------+----------------+----------------------+
The vmin and vmax boundaries are defined as:
- vmin = 0x800000 * 2^(-149)
- vmax = 0xffffff * 2^104
*/
static void init_ieee_table()
{
if (!ieee_table.inited) {
unsigned long i;
unsigned long mmin = 0x800000; // minimum mantissa
unsigned long mmax = 0xffffff; // maximum mantissa
double e = 1;
for (i = 1; i <= 104; i++) {
e *= 2;
ieee_table.e[i + 150] = e;
ieee_table.v[i + 150] = e * mmin;
}
ieee_table.e[150] = 1;
ieee_table.v[150] = mmin;
e = 1;
for (i = 1; i < 150; i++) {
e /= 2;
ieee_table.e[150 - i] = e;
ieee_table.v[150 - i] = e * mmin;
}
ieee_table.vmin = ieee_table.v[1];
ieee_table.vmax = ieee_table.e[254] * mmax;
ieee_table.inited = 1;
/*for (i=0;i<128;i++) printf("++++ ieee_table.v[%d]=%g\n",i,ieee_table.v[i]);*/
}
}
static void init_table_if_needed()
{
GRIB_MUTEX_INIT_ONCE(&once, &init);
GRIB_MUTEX_LOCK(&mutex);
if (!ieee_table.inited)
init_ieee_table();
GRIB_MUTEX_UNLOCK(&mutex);
}
constexpr auto ieee_table = IeeeTable<double>();
static void binary_search(const double xx[], const unsigned long n, double x, unsigned long* j)
{
@ -142,18 +32,6 @@ static void binary_search(const double xx[], const unsigned long n, double x, un
*j = jl;
}
double grib_ieee_table_e(unsigned long e)
{
init_table_if_needed();
return ieee_table.e[e];
}
double grib_ieee_table_v(unsigned long e)
{
init_table_if_needed();
return ieee_table.v[e];
}
unsigned long grib_ieee_to_long(double x)
{
unsigned long s = 0;
@ -163,8 +41,6 @@ unsigned long grib_ieee_to_long(double x)
unsigned long e = 0;
double rmmax = mmax + 0.5;
init_table_if_needed();
/* printf("\ngrib_ieee_to_long: x=%.20e\n",x); */
if (x < 0) {
s = 1;
@ -184,7 +60,7 @@ unsigned long grib_ieee_to_long(double x)
return 0;
}
binary_search(ieee_table.v, 254, x, &e);
binary_search(ieee_table.v.data(), 254, x, &e);
/* printf("grib_ieee_to_long: e=%ld\n",e); */
@ -221,8 +97,6 @@ double grib_ieeefloat_error(double x)
{
unsigned long e = 0;
init_table_if_needed();
if (x < 0)
x = -x;
@ -237,7 +111,7 @@ double grib_ieeefloat_error(double x)
return 0;
}
binary_search(ieee_table.v, 254, x, &e);
binary_search(ieee_table.v.data(), 254, x, &e);
return ieee_table.e[e];
}
@ -256,7 +130,6 @@ double grib_long_to_ieee(unsigned long x)
Assert(0);
}
#endif
init_table_if_needed();
if (c == 0 && m == 0)
return 0;
@ -288,8 +161,6 @@ unsigned long grib_ieee_nearest_smaller_to_long(double x)
if (x == 0)
return 0;
init_table_if_needed();
l = grib_ieee_to_long(x);
y = grib_long_to_ieee(l);
@ -334,10 +205,8 @@ int grib_nearest_smaller_ieee_float(double a, double* ret)
{
unsigned long l = 0;
init_table_if_needed();
if (a > ieee_table.vmax) {
grib_context* c = grib_context_get_default();
const grib_context* c = grib_context_get_default();
grib_context_log(c, GRIB_LOG_ERROR,
"Number is too large: x=%e > xmax=%e (IEEE float)", a, ieee_table.vmax);
return GRIB_INTERNAL_ERROR;

64
src/grib_ieeefloat.h
View File

@ -11,5 +11,69 @@
#pragma once
#include "grib_api_internal.h"
#include "grib_scaling.h"
#include <array>
#include <type_traits>
#include <cstdint>
template <typename T> int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, T* val);
/**
.. _init_ieee_table:
Init IEEE Table
===============
Initializes the ieee_table with IEEE754 single precision (32-bit) values. Nearest smaller values (e.g., reference values for grid_simple and grid_ccsds) are taken from this table.
Details
-------
The table layout is as follows:
+-------+----------------+------------------------+
| idx (i) | multiplier (e) | value (v = mmin * e) |
+-------+----------------+------------------------+
| 0 | 0 | 0 |
| 1 | 2^(-149) | 0x800000 * 2^(-149) |
| 2 | 2^(-148) | 0x800000 * 2^(-148) |
| ... | ... | ... |
| 253 | 2^103 | 0x800000 * 2^103 |
| 254 | 2^104 | 0x800000 * 2^104 |
+-------+----------------+------------------------+
The vmin and vmax boundaries are defined as:
- vmin = 0x800000 * 2^(-149)
- vmax = 0xffffff * 2^104
*/
template <typename ValueType>
struct IeeeTable {
private:
static_assert(std::is_floating_point<ValueType>::value, "ValueType must be a floating point type");
static constexpr int TABLESIZE = 255;
static constexpr uint32_t mantissa_min = 0x800000;
static constexpr uint32_t mantissa_max = 0xffffff;
public:
static constexpr std::array<ValueType, TABLESIZE> e = []() {
std::array<ValueType, TABLESIZE> multiplier{};
multiplier[0] = 0;
for (int i = 1; i < TABLESIZE; ++i) {
multiplier[i] = codes_power<ValueType>(i - 150, 2);
}
return multiplier;
}();
static constexpr std::array<ValueType, TABLESIZE> v = []() {
std::array<ValueType, TABLESIZE> values{};
values[0] = 0;
for (int i = 1; i < TABLESIZE; ++i) {
values[i] = e[i] * mantissa_min;
}
return values;
}();
static constexpr ValueType vmin = e[1] * mantissa_min;
static constexpr ValueType vmax = e[254] * mantissa_max;
};

4
src/grib_scaling.h
View File

@ -1,10 +1,8 @@
#pragma once
template <typename T> T codes_power(long s, long n);
/* Return n to the power of s */
template <typename T>
T codes_power(long s, long n)
constexpr T codes_power(long s, long n)
{
T divisor = 1.0;
if (s == 0)