From f99afd4a66dc22a1d19a26e1445ba276de7d56f3 Mon Sep 17 00:00:00 2001 From: Eugen Betke Date: Tue, 10 Oct 2023 08:12:53 +0000 Subject: [PATCH 1/5] Compile time initializatin of IBM an IEEE tables --- src/grib_ibmfloat.cc | 102 ++-------------------------------- src/grib_ibmfloat.h | 74 +++++++++++++++++++++++++ src/grib_ieeefloat.cc | 125 +++--------------------------------------- src/grib_ieeefloat.h | 61 +++++++++++++++++++++ src/grib_scaling.h | 4 +- 5 files changed, 147 insertions(+), 219 deletions(-) create mode 100644 src/grib_ibmfloat.h diff --git a/src/grib_ibmfloat.cc b/src/grib_ibmfloat.cc index a4b54b358..47f0eda35 100644 --- a/src/grib_ibmfloat.cc +++ b/src/grib_ibmfloat.cc @@ -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; diff --git a/src/grib_ibmfloat.h b/src/grib_ibmfloat.h new file mode 100644 index 000000000..a15736331 --- /dev/null +++ b/src/grib_ibmfloat.h @@ -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_api_internal.h" +#include "grib_scaling.h" + +#include +#include +#include + + +/** +.. _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 uint8_t TABLESIZE = 128; + static constexpr uint32_t mantissa_min = 0x100000; + static constexpr uint32_t mantissa_max = 0xffffff; + +public: + static constexpr std::array e = []() { + std::array multiplier{}; + for (uint8_t i = 1; i < TABLESIZE; ++i) { + multiplier[i] = codes_power(i - 70, 16); + } + return multiplier; + }(); + static constexpr std::array v = []() { + std::array values{}; + for (uint8_t 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; +}; diff --git a/src/grib_ieeefloat.cc b/src/grib_ieeefloat.cc index e638ab4e9..c473f81e8 100644 --- a/src/grib_ieeefloat.cc +++ b/src/grib_ieeefloat.cc @@ -13,116 +13,9 @@ ***************************************************************************/ #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(); static void binary_search(const double xx[], const unsigned long n, double x, unsigned long* j) { @@ -144,13 +37,11 @@ static void binary_search(const double xx[], const unsigned long n, double x, un 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]; } @@ -163,8 +54,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 +73,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,7 +110,7 @@ double grib_ieeefloat_error(double x) { unsigned long e = 0; - init_table_if_needed(); + //init_table_if_needed(); if (x < 0) x = -x; @@ -237,7 +126,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 +145,7 @@ double grib_long_to_ieee(unsigned long x) Assert(0); } #endif - init_table_if_needed(); + //init_table_if_needed(); if (c == 0 && m == 0) return 0; @@ -288,7 +177,7 @@ unsigned long grib_ieee_nearest_smaller_to_long(double x) if (x == 0) return 0; - init_table_if_needed(); + //init_table_if_needed(); l = grib_ieee_to_long(x); y = grib_long_to_ieee(l); @@ -334,7 +223,7 @@ int grib_nearest_smaller_ieee_float(double a, double* ret) { unsigned long l = 0; - init_table_if_needed(); + //init_table_if_needed(); if (a > ieee_table.vmax) { grib_context* c = grib_context_get_default(); diff --git a/src/grib_ieeefloat.h b/src/grib_ieeefloat.h index 74944f2ac..b81ec7b21 100644 --- a/src/grib_ieeefloat.h +++ b/src/grib_ieeefloat.h @@ -11,5 +11,66 @@ #pragma once #include "grib_api_internal.h" +#include "grib_scaling.h" + +#include +#include +#include template 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) | ++-------+----------------+----------------------+ +| 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 +struct IeeeTable { +private: + static_assert(std::is_floating_point::value, "ValueType must be a floating point type"); + static constexpr uint8_t TABLESIZE = 255; + static constexpr uint32_t mantissa_min = 0x800000; + static constexpr uint32_t mantissa_max = 0xffffff; + +public: + static constexpr std::array e = []() { + std::array multiplier{}; + for (uint8_t i = 1; i < TABLESIZE; ++i) { + multiplier[i] = codes_power(i - 150, 2); + } + return multiplier; + }(); + static constexpr std::array v = []() { + std::array values{}; + for (uint8_t 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; +}; diff --git a/src/grib_scaling.h b/src/grib_scaling.h index 53e3b1d4d..c48b3aa2c 100644 --- a/src/grib_scaling.h +++ b/src/grib_scaling.h @@ -1,10 +1,8 @@ #pragma once -template T codes_power(long s, long n); - /* Return n to the power of s */ template -T codes_power(long s, long n) +constexpr T codes_power(long s, long n) { T divisor = 1.0; if (s == 0) From ced0d8187061b8a3776c0169b42dabf03b36929a Mon Sep 17 00:00:00 2001 From: Eugen Betke Date: Tue, 10 Oct 2023 10:59:22 +0000 Subject: [PATCH 2/5] Clean up unused grib_ieee_table_e() and grib_ieee_table_v() --- src/eccodes_prototypes.h | 2 -- src/grib_ieeefloat.cc | 10 ---------- 2 files changed, 12 deletions(-) diff --git a/src/eccodes_prototypes.h b/src/eccodes_prototypes.h index 278fdcf7e..3b405a05a 100644 --- a/src/eccodes_prototypes.h +++ b/src/eccodes_prototypes.h @@ -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); diff --git a/src/grib_ieeefloat.cc b/src/grib_ieeefloat.cc index c473f81e8..a8bfd557a 100644 --- a/src/grib_ieeefloat.cc +++ b/src/grib_ieeefloat.cc @@ -35,16 +35,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) -{ - return ieee_table.e[e]; -} - -double grib_ieee_table_v(unsigned long e) -{ - return ieee_table.v[e]; -} - unsigned long grib_ieee_to_long(double x) { unsigned long s = 0; From 119de10a9da811ca0ddd8cf06f6f4fac99e568bd Mon Sep 17 00:00:00 2001 From: Eugen Betke Date: Tue, 10 Oct 2023 14:31:51 +0000 Subject: [PATCH 3/5] Use Visual Studio 2022 --- .appveyor.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.appveyor.yml b/.appveyor.yml index 4b0a3eaaa..9320dead6 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -9,7 +9,7 @@ branches: - develop - master -image: Visual Studio 2015 +image: Visual Studio 2022 environment: CONDA: c:\Miniconda37-x64 @@ -63,7 +63,7 @@ platform: before_build: - cmd: call "C:\Program Files\Microsoft SDKs\Windows\v7.1\Bin\SetEnv.cmd" /x64 - - cmd: call "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat" x86_amd64 + - cmd: call "C:\Program Files (x86)\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvarsall.bat" x86_amd64 build_script: - cmd: cd %ECCODES_SRC% From a620a9ec4be3929695761d12a7c5d81181a19737 Mon Sep 17 00:00:00 2001 From: Eugen Betke Date: Tue, 10 Oct 2023 14:43:05 +0000 Subject: [PATCH 4/5] Remove path to the SetEnv.bat --- .appveyor.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.appveyor.yml b/.appveyor.yml index 9320dead6..ddb857542 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -62,7 +62,6 @@ platform: - x64 before_build: - - cmd: call "C:\Program Files\Microsoft SDKs\Windows\v7.1\Bin\SetEnv.cmd" /x64 - cmd: call "C:\Program Files (x86)\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvarsall.bat" x86_amd64 build_script: From bcb2606740b9cf22bb8d50a015a54745de4c0d64 Mon Sep 17 00:00:00 2001 From: Eugen Betke Date: Tue, 10 Oct 2023 14:49:46 +0000 Subject: [PATCH 5/5] Change path to vcvarsall.bat --- .appveyor.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.appveyor.yml b/.appveyor.yml index ddb857542..096e0e4c8 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -62,7 +62,7 @@ platform: - x64 before_build: - - cmd: call "C:\Program Files (x86)\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvarsall.bat" x86_amd64 + - cmd: call "C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvarsall.bat" x86_amd64 build_script: - cmd: cd %ECCODES_SRC%