mirror/netcdf-c
2
0
Fork 0
mirror of https://github.com/Unidata/netcdf-c.git synced 2025-01-24 16:04:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

moving quantize to its own function

Browse Source
This commit is contained in:
Edward Hartnett 2021-08-29 08:42:21 -06:00
parent 5aa429cda4
commit ed60a16529
2 changed files with 162 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

173
libsrc4/nc4var.c
View File

@ -466,6 +466,108 @@ NC4_var_par_access(int ncid, int varid, int par_access)
#endif /* USE_PARALLEL4 */
}
/**
* @internal Copy data from one buffer to another, performing
* quantization.
*
* This function applies quantization to float and double data, if
* desired. The code to do this is derived from the bitgroom filter in
* the CCR project (see
* https://github.com/ccr/ccr/blob/master/hdf5_plugins/BITGROOM/src/H5Zbitgroom.c).
*
* @param src Pointer to source of data.
* @param dest Pointer that gets data.
* @param src_type Type ID of source data.
* @param dest_type Type ID of destination data.
* @param len Number of elements of data to copy.
* @param range_error Pointer that gets 1 if there was a range error.
* @param fill_value The fill value.
* @param strict_nc3 Non-zero if strict model in effect.
* @param quantize_mode May be ::NC_NOQUANTIZE or
* ::NC_QUANTIZE_BITGROOM.
* @param nsd Number of significant diggits for quantizize. Ignored
* unless quantize_mode is ::NC_QUANTIZE_BITGROOM.
*
* @returns ::NC_NOERR No error.
* @returns ::NC_EBADTYPE Type not found.
* @author Ed Hartnett, Dennis Heimbigner
*/
int
nc4_quantize_data(const void *src, void *dest, const nc_type src_type,
const nc_type dest_type, const size_t len, int *range_error,
const void *fill_value, int strict_nc3, int quantize_mode,
int nsd)
{
const double bit_per_dcm_dgt_prc = M_LN10 / M_LN2; /* 3.32 [frc] Bits per decimal digit of precision */
const int bit_xpl_nbr_sgn_flt = 23; /* [nbr] Bits 0-22 of SP significands are explicit. Bit 23 is implicitly 1. */
const int bit_xpl_nbr_sgn_dbl = 53; /* [nbr] Bits 0-52 of DP significands are explicit. Bit 53 is implicitly 1. */
double prc_bnr_xct; /* [nbr] Binary digits of precision, exact */
double mss_val_cmp_dbl; /* Missing value for comparison to double precision values */
float mss_val_cmp_flt; /* Missing value for comparison to single precision values */
int bit_xpl_nbr_sgn = -1; /* [nbr] Number of explicit bits in significand */
int bit_xpl_nbr_zro; /* [nbr] Number of explicit bits to zero */
size_t idx;
unsigned int *u32_ptr;
unsigned int msk_f32_u32_zro;
unsigned int msk_f32_u32_one;
unsigned long long int *u64_ptr;
unsigned long long int msk_f64_u64_zro;
unsigned long long int msk_f64_u64_one;
unsigned short prc_bnr_ceil; /* [nbr] Exact binary digits of precision rounded-up */
unsigned short prc_bnr_xpl_rqr; /* [nbr] Explicitly represented binary digits required to retain */
ptr_unn op1; /* I/O [frc] Values to quantize */
float *fp, *fp1;
/* double *dp, *dp1; */
size_t count = 0;
/* How many bits to preserve? */
prc_bnr_xct = nsd * bit_per_dcm_dgt_prc;
/* Be conservative, round upwards */
prc_bnr_ceil =(unsigned short)ceil(prc_bnr_xct);
/* First bit is implicit not explicit but corner cases prevent our taking advantage of this */
prc_bnr_xpl_rqr = prc_bnr_ceil + 1;
if (dest_type == NC_DOUBLE)
prc_bnr_xpl_rqr++; /* Seems necessary for double-precision ppc=array(1.234567,1.0e-6,$dmn) */
if (fill_value)
mss_val_cmp_flt = *(float *)fill_value;
else
mss_val_cmp_flt = NC_FILL_FLOAT;
bit_xpl_nbr_sgn = bit_xpl_nbr_sgn_flt;
bit_xpl_nbr_zro = bit_xpl_nbr_sgn - prc_bnr_xpl_rqr;
assert(bit_xpl_nbr_zro <= bit_xpl_nbr_sgn - NCO_PPC_BIT_XPL_NBR_MIN);
/* Create mask */
msk_f32_u32_zro = 0u; /* Zero all bits */
msk_f32_u32_zro = ~msk_f32_u32_zro; /* Turn all bits to ones */
/* Bit Shave mask for AND: Left shift zeros into bits to be rounded, leave ones in untouched bits */
msk_f32_u32_zro <<= bit_xpl_nbr_zro;
/* Bit Set mask for OR: Put ones into bits to be set, zeros in untouched bits */
msk_f32_u32_one = ~msk_f32_u32_zro;
/* Copy the data into our buffer. */
for (fp = (float *)src, fp1 = dest; count < len; count++)
{
*fp1 = *fp;
/* Move to next float. */
fp1++;
fp++;
}
/* Bit-Groom: alternately shave and set LSBs */
op1.fp = (float *)dest;
u32_ptr = op1.ui32p;
for(idx = 0L; idx < len; idx += 2L)
if (op1.fp[idx] != mss_val_cmp_flt)
u32_ptr[idx] &= msk_f32_u32_zro;
for(idx = 1L; idx < len; idx += 2L)
if (op1.fp[idx] != mss_val_cmp_flt && u32_ptr[idx] != 0U) /* Never quantize upwards floating point values of zero */
u32_ptr[idx] |= msk_f32_u32_one;
return 0;
}
/**
* @internal Copy data from one buffer to another, performing
* appropriate data conversion.
@ -515,6 +617,7 @@ nc4_convert_type(const void *src, void *dest, const nc_type src_type,
long long *lip, *lip1;
unsigned long long *ulip, *ulip1;
size_t count = 0;
int ret;
*range_error = 0;
LOG((3, "%s: len %d src_type %d dest_type %d", __func__, len, src_type,
@ -1172,73 +1275,9 @@ nc4_convert_type(const void *src, void *dest, const nc_type src_type,
case NC_FLOAT:
if (quantize_mode == NC_QUANTIZE_BITGROOM)
{
const double bit_per_dcm_dgt_prc = M_LN10 / M_LN2; /* 3.32 [frc] Bits per decimal digit of precision */
const int bit_xpl_nbr_sgn_flt = 23; /* [nbr] Bits 0-22 of SP significands are explicit. Bit 23 is implicitly 1. */
const int bit_xpl_nbr_sgn_dbl = 53; /* [nbr] Bits 0-52 of DP significands are explicit. Bit 53 is implicitly 1. */
double prc_bnr_xct; /* [nbr] Binary digits of precision, exact */
double mss_val_cmp_dbl; /* Missing value for comparison to double precision values */
float mss_val_cmp_flt; /* Missing value for comparison to single precision values */
int bit_xpl_nbr_sgn = -1; /* [nbr] Number of explicit bits in significand */
int bit_xpl_nbr_zro; /* [nbr] Number of explicit bits to zero */
size_t idx;
unsigned int *u32_ptr;
unsigned int msk_f32_u32_zro;
unsigned int msk_f32_u32_one;
unsigned long long int *u64_ptr;
unsigned long long int msk_f64_u64_zro;
unsigned long long int msk_f64_u64_one;
unsigned short prc_bnr_ceil; /* [nbr] Exact binary digits of precision rounded-up */
unsigned short prc_bnr_xpl_rqr; /* [nbr] Explicitly represented binary digits required to retain */
ptr_unn op1; /* I/O [frc] Values to quantize */
/* How many bits to preserve? */
prc_bnr_xct = nsd * bit_per_dcm_dgt_prc;
/* Be conservative, round upwards */
prc_bnr_ceil =(unsigned short)ceil(prc_bnr_xct);
/* First bit is implicit not explicit but corner cases prevent our taking advantage of this */
prc_bnr_xpl_rqr = prc_bnr_ceil + 1;
if (dest_type == NC_DOUBLE)
prc_bnr_xpl_rqr++; /* Seems necessary for double-precision ppc=array(1.234567,1.0e-6,$dmn) */
if (fill_value)
mss_val_cmp_flt = *(float *)fill_value;
else
mss_val_cmp_flt = NC_FILL_FLOAT;
bit_xpl_nbr_sgn = bit_xpl_nbr_sgn_flt;
bit_xpl_nbr_zro = bit_xpl_nbr_sgn - prc_bnr_xpl_rqr;
assert(bit_xpl_nbr_zro <= bit_xpl_nbr_sgn - NCO_PPC_BIT_XPL_NBR_MIN);
/* Create mask */
msk_f32_u32_zro = 0u; /* Zero all bits */
msk_f32_u32_zro = ~msk_f32_u32_zro; /* Turn all bits to ones */
/* Bit Shave mask for AND: Left shift zeros into bits to be rounded, leave ones in untouched bits */
msk_f32_u32_zro <<= bit_xpl_nbr_zro;
/* Bit Set mask for OR: Put ones into bits to be set, zeros in untouched bits */
msk_f32_u32_one = ~msk_f32_u32_zro;
/* Copy the data into our buffer. */
for (fp = (float *)src, fp1 = dest; count < len; count++)
{
*fp1 = *fp;
/* Move to next float. */
fp1++;
fp++;
}
/* Bit-Groom: alternately shave and set LSBs */
op1.fp = (float *)dest;
u32_ptr = op1.ui32p;
for(idx = 0L; idx < len; idx += 2L)
if (op1.fp[idx] != mss_val_cmp_flt)
u32_ptr[idx] &= msk_f32_u32_zro;
for(idx = 1L; idx < len; idx += 2L)
if (op1.fp[idx] != mss_val_cmp_flt && u32_ptr[idx] != 0U) /* Never quantize upwards floating point values of zero */
u32_ptr[idx] |= msk_f32_u32_one;
if ((ret = nc4_quantize_data(src, dest, src_type, dest_type, len, range_error,
fill_value, strict_nc3, quantize_mode, nsd)))
return ret;
}
else
{

61
nc_test4/tst_quantize.c
View File

@ -25,6 +25,49 @@
#define NSD_3 3
#define NSD_9 9
/* This var used to help print a float in hex. */
char pf_str[20];
/* This struct allows us to treat float as uint32_t
* types. */
union FU {
float f;
uint32_t u;
};
/* This struct allows us to treat double points as uint64_t
* types. */
union DU {
double d;
uint64_t u;
};
/* This function prints a float as hex. */
char *
pf(float myf)
{
union {
float f;
uint32_t u;
} fu;
fu.f = myf;
sprintf(pf_str, "0x%x", fu.u);
return pf_str;
}
/* This function prints a double as hex. */
char *
pd(double myd)
{
union {
double d;
uint64_t u;
} du;
du.d = myd;
sprintf(pf_str, "0x%lx", du.u);
return pf_str;
}
int
main(int argc, char **argv)
{
@ -184,10 +227,8 @@ main(int argc, char **argv)
{
float float_in;
double double_in;
union {
float f;
uint32_t u;
} fin, fout;
union FU fin, fout;
union DU dfin;
/* Open the file and check metadata. */
if (nc_open(FILE_NAME, NC_WRITE, &ncid)) ERR;
@ -201,10 +242,11 @@ main(int argc, char **argv)
if (nc_get_var(ncid, varid2, &double_in)) ERR;
fout.f = float_data[0];
fin.f = float_in;
dfin.d = double_in;
printf ("\nfloat_data: %10f : 0x%x float_data_in: %10f : 0x%x\n",
float_data[0], fout.u, float_data[0], fin.u);
if (fout.u != 0x3f8e38e3) ERR;
if (fin.u != 0x3f8e3000) ERR;
/* if (dfin.u != 0x3ff1c60000000000) ERR; */
/* Close the file again. */
if (nc_close(ncid)) ERR;
@ -246,11 +288,18 @@ main(int argc, char **argv)
union FU fin, fout;
union FU xpect[DIM_LEN_5];
union DU dfin;
union DU double_xpect[DIM_LEN_5];
xpect[0].u = 0x3f8e3000;
xpect[1].u = 0x3f800fff;
xpect[2].u = 0x41200000;
xpect[3].u = 0x4640efff;
xpect[4].u = 0x3dfcd000;
double_xpect[0].u = 0x3ff1c60000000000;
double_xpect[1].u = 0x3ff001ffffffffff;
double_xpect[2].u = 0x4023fe0000000000;
double_xpect[3].u = 0x41d265ffffffffff;
double_xpect[4].u = 0x42dc120000000000;
/* Open the file and check metadata. */
if (nc_open(FILE_NAME, NC_WRITE, &ncid)) ERR;
@ -270,6 +319,8 @@ main(int argc, char **argv)
printf ("float_data: %10f : 0x%x float_data_in: %10f : 0x%x\n",
float_data[x], fout.u, float_data[x], fin.u);
if (fin.u != xpect[x].u) ERR;
dfin.d = double_in[x];
/* if (dfin.u != double_xpect[x].u) ERR; */
}
/* Close the file again. */