NetCDF 4.9.2
Loading...
Searching...
No Matches
nc4var.c
Go to the documentation of this file.
1/* Copyright 2003-2018, University Corporation for Atmospheric
2 * Research. See COPYRIGHT file for copying and redistribution
3 * conditions.*/
13#include "config.h"
14#include "nc4internal.h"
15#include "nc4dispatch.h"
16#ifdef USE_HDF5
17#include "hdf5internal.h"
18#endif
19#include <math.h>
20
22#define DEFAULT_1D_UNLIM_SIZE (4096)
23
24/* Define log_e for 10 and 2. Prefer constants defined in math.h,
25 * however, GCC environments can have hard time defining M_LN10/M_LN2
26 * despite finding math.h */
27#ifndef M_LN10
28# define M_LN10 2.30258509299404568402
29#endif /* M_LN10 */
30#ifndef M_LN2
31# define M_LN2 0.69314718055994530942
32#endif /* M_LN2 */
33
38#define BIT_XPL_NBR_SGN_FLT (23)
39
44#define BIT_XPL_NBR_SGN_DBL (52)
45
47typedef union { /* ptr_unn */
48 float *fp;
49 double *dp;
50 unsigned int *ui32p;
51 unsigned long long *ui64p;
52 void *vp;
53} ptr_unn;
54
71int
72NC4_get_var_chunk_cache(int ncid, int varid, size_t *sizep,
73 size_t *nelemsp, float *preemptionp)
74{
75 NC *nc;
76 NC_GRP_INFO_T *grp;
77 NC_FILE_INFO_T *h5;
78 NC_VAR_INFO_T *var;
79 int retval;
80
81 /* Find info for this file and group, and set pointer to each. */
82 if ((retval = nc4_find_nc_grp_h5(ncid, &nc, &grp, &h5)))
83 return retval;
84 assert(nc && grp && h5);
85
86 /* Find the var. */
87 var = (NC_VAR_INFO_T*)ncindexith(grp->vars,varid);
88 if(!var)
89 return NC_ENOTVAR;
90 assert(var && var->hdr.id == varid);
91
92 /* Give the user what they want. */
93 if (sizep)
94 *sizep = var->chunkcache.size;
95 if (nelemsp)
96 *nelemsp = var->chunkcache.nelems;
97 if (preemptionp)
98 *preemptionp = var->chunkcache.preemption;
99
100 return NC_NOERR;
101}
102
119int
120nc_get_var_chunk_cache_ints(int ncid, int varid, int *sizep,
121 int *nelemsp, int *preemptionp)
122{
123 size_t real_size, real_nelems;
124 float real_preemption;
125 int ret;
126
127 if ((ret = NC4_get_var_chunk_cache(ncid, varid, &real_size,
128 &real_nelems, &real_preemption)))
129 return ret;
130
131 if (sizep)
132 *sizep = real_size / MEGABYTE;
133 if (nelemsp)
134 *nelemsp = (int)real_nelems;
135 if(preemptionp)
136 *preemptionp = (int)(real_preemption * 100);
137
138 return NC_NOERR;
139}
140
178int
179NC4_inq_var_all(int ncid, int varid, char *name, nc_type *xtypep,
180 int *ndimsp, int *dimidsp, int *nattsp,
181 int *shufflep, int *deflatep, int *deflate_levelp,
182 int *fletcher32p, int *storagep, size_t *chunksizesp,
183 int *no_fill, void *fill_valuep, int *endiannessp,
184 unsigned int *idp, size_t *nparamsp, unsigned int *params)
185{
186 NC_GRP_INFO_T *grp;
187 NC_FILE_INFO_T *h5;
188 NC_VAR_INFO_T *var;
189 int d;
190 int retval;
191
192 LOG((2, "%s: ncid 0x%x varid %d", __func__, ncid, varid));
193
194 /* Find info for this file and group, and set pointer to each. */
195 if ((retval = nc4_find_nc_grp_h5(ncid, NULL, &grp, &h5)))
196 return retval;
197 assert(grp && h5);
198
199 /* If the varid is -1, find the global atts and call it a day. */
200 if (varid == NC_GLOBAL && nattsp)
201 {
202 *nattsp = ncindexcount(grp->att);
203 return NC_NOERR;
204 }
205
206 /* Find the var. */
207 if (!(var = (NC_VAR_INFO_T *)ncindexith(grp->vars, varid)))
208 return NC_ENOTVAR;
209 assert(var && var->hdr.id == varid);
210
211 /* Copy the data to the user's data buffers. */
212 if (name)
213 strcpy(name, var->hdr.name);
214 if (xtypep)
215 *xtypep = var->type_info->hdr.id;
216 if (ndimsp)
217 *ndimsp = var->ndims;
218 if (dimidsp)
219 for (d = 0; d < var->ndims; d++)
220 dimidsp[d] = var->dimids[d];
221 if (nattsp)
222 *nattsp = ncindexcount(var->att);
223
224 /* Did the user want the chunksizes? */
225 if (var->storage == NC_CHUNKED && chunksizesp)
226 {
227 for (d = 0; d < var->ndims; d++)
228 {
229 chunksizesp[d] = var->chunksizes[d];
230 LOG((4, "chunksizesp[%d]=%d", d, chunksizesp[d]));
231 }
232 }
233
234 /* Did the user inquire about the storage? */
235 if (storagep)
236 *storagep = var->storage;
237
238 /* Filter stuff. */
239 if (shufflep) {
240 retval = nc_inq_var_filter_info(ncid,varid,H5Z_FILTER_SHUFFLE,0,NULL);
241 if(retval && retval != NC_ENOFILTER) return retval;
242 *shufflep = (retval == NC_NOERR?1:0);
243 }
244 if (fletcher32p) {
245 retval = nc_inq_var_filter_info(ncid,varid,H5Z_FILTER_FLETCHER32,0,NULL);
246 if(retval && retval != NC_ENOFILTER) return retval;
247 *fletcher32p = (retval == NC_NOERR?1:0);
248 }
249 if (deflatep)
250 return NC_EFILTER;
251
252 if (idp) {
253 return NC_EFILTER;
254 }
255
256 /* Fill value stuff. */
257 if (no_fill)
258 *no_fill = (int)var->no_fill;
259
260 /* Don't do a thing with fill_valuep if no_fill mode is set for
261 * this var, or if fill_valuep is NULL. */
262 if (!var->no_fill && fill_valuep)
263 {
264 /* Do we have a fill value for this var? */
265 if (var->fill_value)
266#ifdef SEPDATA
267 {
268 if (var->type_info->nc_type_class == NC_STRING)
269 {
270 assert(*(char **)var->fill_value);
271 /* This will allocate memory and copy the string. */
272 if (!(*(char **)fill_valuep = strdup(*(char **)var->fill_value)))
273 {
274 free(*(char **)fill_valuep);
275 return NC_ENOMEM;
276 }
277 }
278 else
279 {
280 assert(var->type_info->size);
281 memcpy(fill_valuep, var->fill_value, var->type_info->size);
282 }
283 }
284#else
285 {
286 int xtype = var->type_info->hdr.id;
287 if((retval = nc_copy_data(ncid,xtype,var->fill_value,1,fill_valuep))) return retval;
288 }
289#endif
290 else
291 {
292#ifdef SEPDATA
293 if (var->type_info->nc_type_class == NC_STRING)
294 {
295 if (!(*(char **)fill_valuep = calloc(1, sizeof(char *))))
296 return NC_ENOMEM;
297
298 if ((retval = nc4_get_default_fill_value(var->type_info->hdr.ud, (char **)fill_valuep)))
299 {
300 free(*(char **)fill_valuep);
301 return retval;
302 }
303 }
304 else
305 {
306 if ((retval = nc4_get_default_fill_value(var->type_info->hdr.id, fill_valuep)))
307 return retval;
308 }
309#else
310 if ((retval = nc4_get_default_fill_value(var->type_info, fill_valuep)))
311 return retval;
312#endif
313 }
314 }
315
316 /* Does the user want the endianness of this variable? */
317 if (endiannessp)
318 *endiannessp = var->endianness;
319
320 return NC_NOERR;
321}
322
339int
340nc_inq_var_chunking_ints(int ncid, int varid, int *storagep, int *chunksizesp)
341{
342 NC_VAR_INFO_T *var;
343 size_t *cs = NULL;
344 int i, retval;
345
346 /* Get pointer to the var. */
347 if ((retval = nc4_find_grp_h5_var(ncid, varid, NULL, NULL, &var)))
348 return retval;
349 assert(var);
350
351 /* Allocate space for the size_t copy of the chunksizes array. */
352 if (var->ndims)
353 if (!(cs = malloc(var->ndims * sizeof(size_t))))
354 return NC_ENOMEM;
355
356 /* Call the netcdf-4 version directly. */
357 retval = NC4_inq_var_all(ncid, varid, NULL, NULL, NULL, NULL, NULL,
358 NULL, NULL, NULL, NULL, storagep, cs, NULL,
359 NULL, NULL, NULL, NULL, NULL);
360
361 /* Copy from size_t array. */
362 if (!retval && chunksizesp && var->storage == NC_CHUNKED)
363 {
364 for (i = 0; i < var->ndims; i++)
365 {
366 chunksizesp[i] = (int)cs[i];
367 if (cs[i] > NC_MAX_INT)
368 retval = NC_ERANGE;
369 }
370 }
371
372 if (var->ndims)
373 free(cs);
374 return retval;
375}
376
389int
390NC4_inq_varid(int ncid, const char *name, int *varidp)
391{
392 NC *nc;
393 NC_GRP_INFO_T *grp;
394 NC_VAR_INFO_T *var;
395 char norm_name[NC_MAX_NAME + 1];
396 int retval;
397
398 if (!name)
399 return NC_EINVAL;
400 if (!varidp)
401 return NC_NOERR;
402
403 LOG((2, "%s: ncid 0x%x name %s", __func__, ncid, name));
404
405 /* Find info for this file and group, and set pointer to each. */
406 if ((retval = nc4_find_nc_grp_h5(ncid, &nc, &grp, NULL)))
407 return retval;
408
409 /* Normalize name. */
410 if ((retval = nc4_normalize_name(name, norm_name)))
411 return retval;
412
413 /* Find var of this name. */
414 var = (NC_VAR_INFO_T*)ncindexlookup(grp->vars,norm_name);
415 if(var)
416 {
417 *varidp = var->hdr.id;
418 return NC_NOERR;
419 }
420 return NC_ENOTVAR;
421}
422
441int
442NC4_var_par_access(int ncid, int varid, int par_access)
443{
444#ifndef USE_PARALLEL4
445 NC_UNUSED(ncid);
446 NC_UNUSED(varid);
447 NC_UNUSED(par_access);
448 return NC_ENOPAR;
449#else
450 NC *nc;
451 NC_GRP_INFO_T *grp;
452 NC_FILE_INFO_T *h5;
453 NC_VAR_INFO_T *var;
454 int retval;
455
456 LOG((1, "%s: ncid 0x%x varid %d par_access %d", __func__, ncid,
457 varid, par_access));
458
459 if (par_access != NC_INDEPENDENT && par_access != NC_COLLECTIVE)
460 return NC_EINVAL;
461
462 /* Find info for this file and group, and set pointer to each. */
463 if ((retval = nc4_find_nc_grp_h5(ncid, &nc, &grp, &h5)))
464 return retval;
465
466 /* This function only for files opened with nc_open_par or nc_create_par. */
467 if (!h5->parallel)
468 return NC_ENOPAR;
469
470 /* Find the var, and set its preference. */
471 var = (NC_VAR_INFO_T*)ncindexith(grp->vars,varid);
472 if (!var) return NC_ENOTVAR;
473 assert(var->hdr.id == varid);
474
475 /* If zlib, shuffle, or fletcher32 filters are in use, then access
476 * must be collective. Fail an attempt to set such a variable to
477 * independent access. */
478 if (nclistlength((NClist*)var->filters) > 0 &&
479 par_access == NC_INDEPENDENT)
480 return NC_EINVAL;
481
482 if (par_access)
483 var->parallel_access = NC_COLLECTIVE;
484 else
485 var->parallel_access = NC_INDEPENDENT;
486 return NC_NOERR;
487#endif /* USE_PARALLEL4 */
488}
489
522int
523nc4_convert_type(const void *src, void *dest, const nc_type src_type,
524 const nc_type dest_type, const size_t len, int *range_error,
525 const void *fill_value, int strict_nc3, int quantize_mode,
526 int nsd)
527{
528 /* These vars are used with quantize feature. */
529 const double bit_per_dgt = M_LN10 / M_LN2; /* 3.32 [frc] Bits per decimal digit of precision = log2(10) */
530 const double dgt_per_bit= M_LN2 / M_LN10; /* 0.301 [frc] Decimal digits per bit of precision = log10(2) */
531 double mnt; /* [frc] Mantissa, 0.5 <= mnt < 1.0 */
532 double mnt_fabs; /* [frc] fabs(mantissa) */
533 double mnt_log10_fabs; /* [frc] log10(fabs(mantissa))) */
534 double val; /* [frc] Copy of input value to avoid indirection */
535 double mss_val_cmp_dbl; /* Missing value for comparison to double precision values */
536 float mss_val_cmp_flt; /* Missing value for comparison to single precision values */
537 int bit_xpl_nbr_zro; /* [nbr] Number of explicit bits to zero */
538 int dgt_nbr; /* [nbr] Number of digits before decimal point */
539 int qnt_pwr; /* [nbr] Power of two in quantization mask: qnt_msk = 2^qnt_pwr */
540 int xpn_bs2; /* [nbr] Binary exponent xpn_bs2 in val = sign(val) * 2^xpn_bs2 * mnt, 0.5 < mnt <= 1.0 */
541 size_t idx;
542 unsigned int *u32_ptr;
543 unsigned int msk_f32_u32_zro;
544 unsigned int msk_f32_u32_one;
545 unsigned int msk_f32_u32_hshv;
546 unsigned long long int *u64_ptr;
547 unsigned long long int msk_f64_u64_zro;
548 unsigned long long int msk_f64_u64_one;
549 unsigned long long int msk_f64_u64_hshv;
550 unsigned short prc_bnr_xpl_rqr; /* [nbr] Explicitly represented binary digits required to retain */
551 ptr_unn op1; /* I/O [frc] Values to quantize */
552
553 char *cp, *cp1;
554 float *fp, *fp1;
555 double *dp, *dp1;
556 int *ip, *ip1;
557 short *sp, *sp1;
558 signed char *bp, *bp1;
559 unsigned char *ubp, *ubp1;
560 unsigned short *usp, *usp1;
561 unsigned int *uip, *uip1;
562 long long *lip, *lip1;
563 unsigned long long *ulip, *ulip1;
564 size_t count = 0;
565
566 *range_error = 0;
567 LOG((3, "%s: len %d src_type %d dest_type %d", __func__, len, src_type,
568 dest_type));
569
570 /* If quantize is in use, set up some values. Quantize can only be
571 * used when the destination type is NC_FLOAT or NC_DOUBLE. */
572 if (quantize_mode != NC_NOQUANTIZE)
573 {
574 assert(dest_type == NC_FLOAT || dest_type == NC_DOUBLE);
575
576 /* Parameters shared by all quantization codecs */
577 if (dest_type == NC_FLOAT)
578 {
579 /* Determine the fill value. */
580 if (fill_value)
581 mss_val_cmp_flt = *(float *)fill_value;
582 else
583 mss_val_cmp_flt = NC_FILL_FLOAT;
584
585 }
586 else
587 {
588
589 /* Determine the fill value. */
590 if (fill_value)
591 mss_val_cmp_dbl = *(double *)fill_value;
592 else
593 mss_val_cmp_dbl = NC_FILL_DOUBLE;
594
595 }
596
597 /* Set parameters used by BitGroom and BitRound here, outside value loop.
598 Equivalent parameters used by GranularBR are set inside value loop,
599 since keep bits and thus masks can change for every value. */
600 if (quantize_mode == NC_QUANTIZE_BITGROOM ||
601 quantize_mode == NC_QUANTIZE_BITROUND )
602 {
603
604 if (quantize_mode == NC_QUANTIZE_BITGROOM){
605
606 /* BitGroom interprets nsd as number of significant decimal digits
607 * Must convert that to number of significant bits to preserve
608 * How many bits to preserve? Being conservative, we round up the
609 * exact binary digits of precision. Add one because the first bit
610 * is implicit not explicit but corner cases prevent our taking
611 * advantage of this. */
612 prc_bnr_xpl_rqr = (unsigned short)ceil(nsd * bit_per_dgt) + 1;
613
614 }else if (quantize_mode == NC_QUANTIZE_BITROUND){
615
616 /* BitRound interprets nsd as number of significant binary digits (bits) */
617 prc_bnr_xpl_rqr = nsd;
618
619 }
620
621 if (dest_type == NC_FLOAT)
622 {
623
624 bit_xpl_nbr_zro = BIT_XPL_NBR_SGN_FLT - prc_bnr_xpl_rqr;
625
626 /* Create mask */
627 msk_f32_u32_zro = 0u; /* Zero all bits */
628 msk_f32_u32_zro = ~msk_f32_u32_zro; /* Turn all bits to ones */
629
630 /* BitShave mask for AND: Left shift zeros into bits to be
631 * rounded, leave ones in untouched bits. */
632 msk_f32_u32_zro <<= bit_xpl_nbr_zro;
633
634 /* BitSet mask for OR: Put ones into bits to be set, zeros in
635 * untouched bits. */
636 msk_f32_u32_one = ~msk_f32_u32_zro;
637
638 /* BitRound mask for ADD: Set one bit: the MSB of LSBs */
639 msk_f32_u32_hshv=msk_f32_u32_one & (msk_f32_u32_zro >> 1);
640
641 }
642 else
643 {
644
645 bit_xpl_nbr_zro = BIT_XPL_NBR_SGN_DBL - prc_bnr_xpl_rqr;
646 /* Create mask. */
647 msk_f64_u64_zro = 0ul; /* Zero all bits. */
648 msk_f64_u64_zro = ~msk_f64_u64_zro; /* Turn all bits to ones. */
649
650 /* BitShave mask for AND: Left shift zeros into bits to be
651 * rounded, leave ones in untouched bits. */
652 msk_f64_u64_zro <<= bit_xpl_nbr_zro;
653
654 /* BitSet mask for OR: Put ones into bits to be set, zeros in
655 * untouched bits. */
656 msk_f64_u64_one =~ msk_f64_u64_zro;
657
658 /* BitRound mask for ADD: Set one bit: the MSB of LSBs */
659 msk_f64_u64_hshv = msk_f64_u64_one & (msk_f64_u64_zro >> 1);
660
661 }
662
663 }
664
665 } /* endif quantize */
666
667 /* OK, this is ugly. If you can think of anything better, I'm open
668 to suggestions!
669
670 Note that we don't use a default fill value for type
671 NC_BYTE. This is because Lord Voldemort cast a nofilleramous spell
672 at Harry Potter, but it bounced off his scar and hit the netcdf-4
673 code.
674 */
675 switch (src_type)
676 {
677 case NC_CHAR:
678 switch (dest_type)
679 {
680 case NC_CHAR:
681 for (cp = (char *)src, cp1 = dest; count < len; count++)
682 *cp1++ = *cp++;
683 break;
684 default:
685 LOG((0, "%s: Unknown destination type.", __func__));
686 }
687 break;
688
689 case NC_BYTE:
690 switch (dest_type)
691 {
692 case NC_BYTE:
693 for (bp = (signed char *)src, bp1 = dest; count < len; count++)
694 *bp1++ = *bp++;
695 break;
696 case NC_UBYTE:
697 for (bp = (signed char *)src, ubp = dest; count < len; count++)
698 {
699 if (*bp < 0)
700 (*range_error)++;
701 *ubp++ = *bp++;
702 }
703 break;
704 case NC_SHORT:
705 for (bp = (signed char *)src, sp = dest; count < len; count++)
706 *sp++ = *bp++;
707 break;
708 case NC_USHORT:
709 for (bp = (signed char *)src, usp = dest; count < len; count++)
710 {
711 if (*bp < 0)
712 (*range_error)++;
713 *usp++ = *bp++;
714 }
715 break;
716 case NC_INT:
717 for (bp = (signed char *)src, ip = dest; count < len; count++)
718 *ip++ = *bp++;
719 break;
720 case NC_UINT:
721 for (bp = (signed char *)src, uip = dest; count < len; count++)
722 {
723 if (*bp < 0)
724 (*range_error)++;
725 *uip++ = *bp++;
726 }
727 break;
728 case NC_INT64:
729 for (bp = (signed char *)src, lip = dest; count < len; count++)
730 *lip++ = *bp++;
731 break;
732 case NC_UINT64:
733 for (bp = (signed char *)src, ulip = dest; count < len; count++)
734 {
735 if (*bp < 0)
736 (*range_error)++;
737 *ulip++ = *bp++;
738 }
739 break;
740 case NC_FLOAT:
741 for (bp = (signed char *)src, fp = dest; count < len; count++)
742 *fp++ = *bp++;
743 break;
744 case NC_DOUBLE:
745 for (bp = (signed char *)src, dp = dest; count < len; count++)
746 *dp++ = *bp++;
747 break;
748 default:
749 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
750 __func__, src_type, dest_type));
751 return NC_EBADTYPE;
752 }
753 break;
754
755 case NC_UBYTE:
756 switch (dest_type)
757 {
758 case NC_BYTE:
759 for (ubp = (unsigned char *)src, bp = dest; count < len; count++)
760 {
761 if (!strict_nc3 && *ubp > X_SCHAR_MAX)
762 (*range_error)++;
763 *bp++ = *ubp++;
764 }
765 break;
766 case NC_SHORT:
767 for (ubp = (unsigned char *)src, sp = dest; count < len; count++)
768 *sp++ = *ubp++;
769 break;
770 case NC_UBYTE:
771 for (ubp = (unsigned char *)src, ubp1 = dest; count < len; count++)
772 *ubp1++ = *ubp++;
773 break;
774 case NC_USHORT:
775 for (ubp = (unsigned char *)src, usp = dest; count < len; count++)
776 *usp++ = *ubp++;
777 break;
778 case NC_INT:
779 for (ubp = (unsigned char *)src, ip = dest; count < len; count++)
780 *ip++ = *ubp++;
781 break;
782 case NC_UINT:
783 for (ubp = (unsigned char *)src, uip = dest; count < len; count++)
784 *uip++ = *ubp++;
785 break;
786 case NC_INT64:
787 for (ubp = (unsigned char *)src, lip = dest; count < len; count++)
788 *lip++ = *ubp++;
789 break;
790 case NC_UINT64:
791 for (ubp = (unsigned char *)src, ulip = dest; count < len; count++)
792 *ulip++ = *ubp++;
793 break;
794 case NC_FLOAT:
795 for (ubp = (unsigned char *)src, fp = dest; count < len; count++)
796 *fp++ = *ubp++;
797 break;
798 case NC_DOUBLE:
799 for (ubp = (unsigned char *)src, dp = dest; count < len; count++)
800 *dp++ = *ubp++;
801 break;
802 default:
803 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
804 __func__, src_type, dest_type));
805 return NC_EBADTYPE;
806 }
807 break;
808
809 case NC_SHORT:
810 switch (dest_type)
811 {
812 case NC_UBYTE:
813 for (sp = (short *)src, ubp = dest; count < len; count++)
814 {
815 if (*sp > X_UCHAR_MAX || *sp < 0)
816 (*range_error)++;
817 *ubp++ = *sp++;
818 }
819 break;
820 case NC_BYTE:
821 for (sp = (short *)src, bp = dest; count < len; count++)
822 {
823 if (*sp > X_SCHAR_MAX || *sp < X_SCHAR_MIN)
824 (*range_error)++;
825 *bp++ = *sp++;
826 }
827 break;
828 case NC_SHORT:
829 for (sp = (short *)src, sp1 = dest; count < len; count++)
830 *sp1++ = *sp++;
831 break;
832 case NC_USHORT:
833 for (sp = (short *)src, usp = dest; count < len; count++)
834 {
835 if (*sp < 0)
836 (*range_error)++;
837 *usp++ = *sp++;
838 }
839 break;
840 case NC_INT:
841 for (sp = (short *)src, ip = dest; count < len; count++)
842 *ip++ = *sp++;
843 break;
844 case NC_UINT:
845 for (sp = (short *)src, uip = dest; count < len; count++)
846 {
847 if (*sp < 0)
848 (*range_error)++;
849 *uip++ = *sp++;
850 }
851 break;
852 case NC_INT64:
853 for (sp = (short *)src, lip = dest; count < len; count++)
854 *lip++ = *sp++;
855 break;
856 case NC_UINT64:
857 for (sp = (short *)src, ulip = dest; count < len; count++)
858 {
859 if (*sp < 0)
860 (*range_error)++;
861 *ulip++ = *sp++;
862 }
863 break;
864 case NC_FLOAT:
865 for (sp = (short *)src, fp = dest; count < len; count++)
866 *fp++ = *sp++;
867 break;
868 case NC_DOUBLE:
869 for (sp = (short *)src, dp = dest; count < len; count++)
870 *dp++ = *sp++;
871 break;
872 default:
873 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
874 __func__, src_type, dest_type));
875 return NC_EBADTYPE;
876 }
877 break;
878
879 case NC_USHORT:
880 switch (dest_type)
881 {
882 case NC_UBYTE:
883 for (usp = (unsigned short *)src, ubp = dest; count < len; count++)
884 {
885 if (*usp > X_UCHAR_MAX)
886 (*range_error)++;
887 *ubp++ = *usp++;
888 }
889 break;
890 case NC_BYTE:
891 for (usp = (unsigned short *)src, bp = dest; count < len; count++)
892 {
893 if (*usp > X_SCHAR_MAX)
894 (*range_error)++;
895 *bp++ = *usp++;
896 }
897 break;
898 case NC_SHORT:
899 for (usp = (unsigned short *)src, sp = dest; count < len; count++)
900 {
901 if (*usp > X_SHORT_MAX)
902 (*range_error)++;
903 *sp++ = *usp++;
904 }
905 break;
906 case NC_USHORT:
907 for (usp = (unsigned short *)src, usp1 = dest; count < len; count++)
908 *usp1++ = *usp++;
909 break;
910 case NC_INT:
911 for (usp = (unsigned short *)src, ip = dest; count < len; count++)
912 *ip++ = *usp++;
913 break;
914 case NC_UINT:
915 for (usp = (unsigned short *)src, uip = dest; count < len; count++)
916 *uip++ = *usp++;
917 break;
918 case NC_INT64:
919 for (usp = (unsigned short *)src, lip = dest; count < len; count++)
920 *lip++ = *usp++;
921 break;
922 case NC_UINT64:
923 for (usp = (unsigned short *)src, ulip = dest; count < len; count++)
924 *ulip++ = *usp++;
925 break;
926 case NC_FLOAT:
927 for (usp = (unsigned short *)src, fp = dest; count < len; count++)
928 *fp++ = *usp++;
929 break;
930 case NC_DOUBLE:
931 for (usp = (unsigned short *)src, dp = dest; count < len; count++)
932 *dp++ = *usp++;
933 break;
934 default:
935 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
936 __func__, src_type, dest_type));
937 return NC_EBADTYPE;
938 }
939 break;
940
941 case NC_INT:
942 switch (dest_type)
943 {
944 case NC_UBYTE:
945 for (ip = (int *)src, ubp = dest; count < len; count++)
946 {
947 if (*ip > X_UCHAR_MAX || *ip < 0)
948 (*range_error)++;
949 *ubp++ = *ip++;
950 }
951 break;
952 case NC_BYTE:
953 for (ip = (int *)src, bp = dest; count < len; count++)
954 {
955 if (*ip > X_SCHAR_MAX || *ip < X_SCHAR_MIN)
956 (*range_error)++;
957 *bp++ = *ip++;
958 }
959 break;
960 case NC_SHORT:
961 for (ip = (int *)src, sp = dest; count < len; count++)
962 {
963 if (*ip > X_SHORT_MAX || *ip < X_SHORT_MIN)
964 (*range_error)++;
965 *sp++ = *ip++;
966 }
967 break;
968 case NC_USHORT:
969 for (ip = (int *)src, usp = dest; count < len; count++)
970 {
971 if (*ip > X_USHORT_MAX || *ip < 0)
972 (*range_error)++;
973 *usp++ = *ip++;
974 }
975 break;
976 case NC_INT: /* src is int */
977 for (ip = (int *)src, ip1 = dest; count < len; count++)
978 {
979 if (*ip > X_INT_MAX || *ip < X_INT_MIN)
980 (*range_error)++;
981 *ip1++ = *ip++;
982 }
983 break;
984 case NC_UINT:
985 for (ip = (int *)src, uip = dest; count < len; count++)
986 {
987 if (*ip > X_UINT_MAX || *ip < 0)
988 (*range_error)++;
989 *uip++ = *ip++;
990 }
991 break;
992 case NC_INT64:
993 for (ip = (int *)src, lip = dest; count < len; count++)
994 *lip++ = *ip++;
995 break;
996 case NC_UINT64:
997 for (ip = (int *)src, ulip = dest; count < len; count++)
998 {
999 if (*ip < 0)
1000 (*range_error)++;
1001 *ulip++ = *ip++;
1002 }
1003 break;
1004 case NC_FLOAT:
1005 for (ip = (int *)src, fp = dest; count < len; count++)
1006 *fp++ = *ip++;
1007 break;
1008 case NC_DOUBLE:
1009 for (ip = (int *)src, dp = dest; count < len; count++)
1010 *dp++ = *ip++;
1011 break;
1012 default:
1013 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1014 __func__, src_type, dest_type));
1015 return NC_EBADTYPE;
1016 }
1017 break;
1018
1019 case NC_UINT:
1020 switch (dest_type)
1021 {
1022 case NC_UBYTE:
1023 for (uip = (unsigned int *)src, ubp = dest; count < len; count++)
1024 {
1025 if (*uip > X_UCHAR_MAX)
1026 (*range_error)++;
1027 *ubp++ = *uip++;
1028 }
1029 break;
1030 case NC_BYTE:
1031 for (uip = (unsigned int *)src, bp = dest; count < len; count++)
1032 {
1033 if (*uip > X_SCHAR_MAX)
1034 (*range_error)++;
1035 *bp++ = *uip++;
1036 }
1037 break;
1038 case NC_SHORT:
1039 for (uip = (unsigned int *)src, sp = dest; count < len; count++)
1040 {
1041 if (*uip > X_SHORT_MAX)
1042 (*range_error)++;
1043 *sp++ = *uip++;
1044 }
1045 break;
1046 case NC_USHORT:
1047 for (uip = (unsigned int *)src, usp = dest; count < len; count++)
1048 {
1049 if (*uip > X_USHORT_MAX)
1050 (*range_error)++;
1051 *usp++ = *uip++;
1052 }
1053 break;
1054 case NC_INT:
1055 for (uip = (unsigned int *)src, ip = dest; count < len; count++)
1056 {
1057 if (*uip > X_INT_MAX)
1058 (*range_error)++;
1059 *ip++ = *uip++;
1060 }
1061 break;
1062 case NC_UINT:
1063 for (uip = (unsigned int *)src, uip1 = dest; count < len; count++)
1064 {
1065 if (*uip > X_UINT_MAX)
1066 (*range_error)++;
1067 *uip1++ = *uip++;
1068 }
1069 break;
1070 case NC_INT64:
1071 for (uip = (unsigned int *)src, lip = dest; count < len; count++)
1072 *lip++ = *uip++;
1073 break;
1074 case NC_UINT64:
1075 for (uip = (unsigned int *)src, ulip = dest; count < len; count++)
1076 *ulip++ = *uip++;
1077 break;
1078 case NC_FLOAT:
1079 for (uip = (unsigned int *)src, fp = dest; count < len; count++)
1080 *fp++ = *uip++;
1081 break;
1082 case NC_DOUBLE:
1083 for (uip = (unsigned int *)src, dp = dest; count < len; count++)
1084 *dp++ = *uip++;
1085 break;
1086 default:
1087 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1088 __func__, src_type, dest_type));
1089 return NC_EBADTYPE;
1090 }
1091 break;
1092
1093 case NC_INT64:
1094 switch (dest_type)
1095 {
1096 case NC_UBYTE:
1097 for (lip = (long long *)src, ubp = dest; count < len; count++)
1098 {
1099 if (*lip > X_UCHAR_MAX || *lip < 0)
1100 (*range_error)++;
1101 *ubp++ = *lip++;
1102 }
1103 break;
1104 case NC_BYTE:
1105 for (lip = (long long *)src, bp = dest; count < len; count++)
1106 {
1107 if (*lip > X_SCHAR_MAX || *lip < X_SCHAR_MIN)
1108 (*range_error)++;
1109 *bp++ = *lip++;
1110 }
1111 break;
1112 case NC_SHORT:
1113 for (lip = (long long *)src, sp = dest; count < len; count++)
1114 {
1115 if (*lip > X_SHORT_MAX || *lip < X_SHORT_MIN)
1116 (*range_error)++;
1117 *sp++ = *lip++;
1118 }
1119 break;
1120 case NC_USHORT:
1121 for (lip = (long long *)src, usp = dest; count < len; count++)
1122 {
1123 if (*lip > X_USHORT_MAX || *lip < 0)
1124 (*range_error)++;
1125 *usp++ = *lip++;
1126 }
1127 break;
1128 case NC_UINT:
1129 for (lip = (long long *)src, uip = dest; count < len; count++)
1130 {
1131 if (*lip > X_UINT_MAX || *lip < 0)
1132 (*range_error)++;
1133 *uip++ = *lip++;
1134 }
1135 break;
1136 case NC_INT:
1137 for (lip = (long long *)src, ip = dest; count < len; count++)
1138 {
1139 if (*lip > X_INT_MAX || *lip < X_INT_MIN)
1140 (*range_error)++;
1141 *ip++ = *lip++;
1142 }
1143 break;
1144 case NC_INT64:
1145 for (lip = (long long *)src, lip1 = dest; count < len; count++)
1146 *lip1++ = *lip++;
1147 break;
1148 case NC_UINT64:
1149 for (lip = (long long *)src, ulip = dest; count < len; count++)
1150 {
1151 if (*lip < 0)
1152 (*range_error)++;
1153 *ulip++ = *lip++;
1154 }
1155 break;
1156 case NC_FLOAT:
1157 for (lip = (long long *)src, fp = dest; count < len; count++)
1158 *fp++ = *lip++;
1159 break;
1160 case NC_DOUBLE:
1161 for (lip = (long long *)src, dp = dest; count < len; count++)
1162 *dp++ = *lip++;
1163 break;
1164 default:
1165 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1166 __func__, src_type, dest_type));
1167 return NC_EBADTYPE;
1168 }
1169 break;
1170
1171 case NC_UINT64:
1172 switch (dest_type)
1173 {
1174 case NC_UBYTE:
1175 for (ulip = (unsigned long long *)src, ubp = dest; count < len; count++)
1176 {
1177 if (*ulip > X_UCHAR_MAX)
1178 (*range_error)++;
1179 *ubp++ = *ulip++;
1180 }
1181 break;
1182 case NC_BYTE:
1183 for (ulip = (unsigned long long *)src, bp = dest; count < len; count++)
1184 {
1185 if (*ulip > X_SCHAR_MAX)
1186 (*range_error)++;
1187 *bp++ = *ulip++;
1188 }
1189 break;
1190 case NC_SHORT:
1191 for (ulip = (unsigned long long *)src, sp = dest; count < len; count++)
1192 {
1193 if (*ulip > X_SHORT_MAX)
1194 (*range_error)++;
1195 *sp++ = *ulip++;
1196 }
1197 break;
1198 case NC_USHORT:
1199 for (ulip = (unsigned long long *)src, usp = dest; count < len; count++)
1200 {
1201 if (*ulip > X_USHORT_MAX)
1202 (*range_error)++;
1203 *usp++ = *ulip++;
1204 }
1205 break;
1206 case NC_UINT:
1207 for (ulip = (unsigned long long *)src, uip = dest; count < len; count++)
1208 {
1209 if (*ulip > X_UINT_MAX)
1210 (*range_error)++;
1211 *uip++ = *ulip++;
1212 }
1213 break;
1214 case NC_INT:
1215 for (ulip = (unsigned long long *)src, ip = dest; count < len; count++)
1216 {
1217 if (*ulip > X_INT_MAX)
1218 (*range_error)++;
1219 *ip++ = *ulip++;
1220 }
1221 break;
1222 case NC_INT64:
1223 for (ulip = (unsigned long long *)src, lip = dest; count < len; count++)
1224 {
1225 if (*ulip > X_INT64_MAX)
1226 (*range_error)++;
1227 *lip++ = *ulip++;
1228 }
1229 break;
1230 case NC_UINT64:
1231 for (ulip = (unsigned long long *)src, ulip1 = dest; count < len; count++)
1232 *ulip1++ = *ulip++;
1233 break;
1234 case NC_FLOAT:
1235 for (ulip = (unsigned long long *)src, fp = dest; count < len; count++)
1236 *fp++ = *ulip++;
1237 break;
1238 case NC_DOUBLE:
1239 for (ulip = (unsigned long long *)src, dp = dest; count < len; count++)
1240 *dp++ = *ulip++;
1241 break;
1242 default:
1243 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1244 __func__, src_type, dest_type));
1245 return NC_EBADTYPE;
1246 }
1247 break;
1248
1249 case NC_FLOAT:
1250 switch (dest_type)
1251 {
1252 case NC_UBYTE:
1253 for (fp = (float *)src, ubp = dest; count < len; count++)
1254 {
1255 if (*fp > X_UCHAR_MAX || *fp < 0)
1256 (*range_error)++;
1257 *ubp++ = *fp++;
1258 }
1259 break;
1260 case NC_BYTE:
1261 for (fp = (float *)src, bp = dest; count < len; count++)
1262 {
1263 if (*fp > (double)X_SCHAR_MAX || *fp < (double)X_SCHAR_MIN)
1264 (*range_error)++;
1265 *bp++ = *fp++;
1266 }
1267 break;
1268 case NC_SHORT:
1269 for (fp = (float *)src, sp = dest; count < len; count++)
1270 {
1271 if (*fp > (double)X_SHORT_MAX || *fp < (double)X_SHORT_MIN)
1272 (*range_error)++;
1273 *sp++ = *fp++;
1274 }
1275 break;
1276 case NC_USHORT:
1277 for (fp = (float *)src, usp = dest; count < len; count++)
1278 {
1279 if (*fp > X_USHORT_MAX || *fp < 0)
1280 (*range_error)++;
1281 *usp++ = *fp++;
1282 }
1283 break;
1284 case NC_UINT:
1285 for (fp = (float *)src, uip = dest; count < len; count++)
1286 {
1287 if (*fp > X_UINT_MAX || *fp < 0)
1288 (*range_error)++;
1289 *uip++ = *fp++;
1290 }
1291 break;
1292 case NC_INT:
1293 for (fp = (float *)src, ip = dest; count < len; count++)
1294 {
1295 if (*fp > (double)X_INT_MAX || *fp < (double)X_INT_MIN)
1296 (*range_error)++;
1297 *ip++ = *fp++;
1298 }
1299 break;
1300 case NC_INT64:
1301 for (fp = (float *)src, lip = dest; count < len; count++)
1302 {
1303 if (*fp > X_INT64_MAX || *fp <X_INT64_MIN)
1304 (*range_error)++;
1305 *lip++ = *fp++;
1306 }
1307 break;
1308 case NC_UINT64:
1309 for (fp = (float *)src, lip = dest; count < len; count++)
1310 {
1311 if (*fp > X_UINT64_MAX || *fp < 0)
1312 (*range_error)++;
1313 *lip++ = *fp++;
1314 }
1315 break;
1316 case NC_FLOAT:
1317 for (fp = (float *)src, fp1 = dest; count < len; count++)
1318 *fp1++ = *fp++;
1319 break;
1320 case NC_DOUBLE:
1321 for (fp = (float *)src, dp = dest; count < len; count++)
1322 *dp++ = *fp++;
1323 break;
1324 default:
1325 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1326 __func__, src_type, dest_type));
1327 return NC_EBADTYPE;
1328 }
1329 break;
1330
1331 case NC_DOUBLE:
1332 switch (dest_type)
1333 {
1334 case NC_UBYTE:
1335 for (dp = (double *)src, ubp = dest; count < len; count++)
1336 {
1337 if (*dp > X_UCHAR_MAX || *dp < 0)
1338 (*range_error)++;
1339 *ubp++ = *dp++;
1340 }
1341 break;
1342 case NC_BYTE:
1343 for (dp = (double *)src, bp = dest; count < len; count++)
1344 {
1345 if (*dp > X_SCHAR_MAX || *dp < X_SCHAR_MIN)
1346 (*range_error)++;
1347 *bp++ = *dp++;
1348 }
1349 break;
1350 case NC_SHORT:
1351 for (dp = (double *)src, sp = dest; count < len; count++)
1352 {
1353 if (*dp > X_SHORT_MAX || *dp < X_SHORT_MIN)
1354 (*range_error)++;
1355 *sp++ = *dp++;
1356 }
1357 break;
1358 case NC_USHORT:
1359 for (dp = (double *)src, usp = dest; count < len; count++)
1360 {
1361 if (*dp > X_USHORT_MAX || *dp < 0)
1362 (*range_error)++;
1363 *usp++ = *dp++;
1364 }
1365 break;
1366 case NC_UINT:
1367 for (dp = (double *)src, uip = dest; count < len; count++)
1368 {
1369 if (*dp > X_UINT_MAX || *dp < 0)
1370 (*range_error)++;
1371 *uip++ = *dp++;
1372 }
1373 break;
1374 case NC_INT:
1375 for (dp = (double *)src, ip = dest; count < len; count++)
1376 {
1377 if (*dp > X_INT_MAX || *dp < X_INT_MIN)
1378 (*range_error)++;
1379 *ip++ = *dp++;
1380 }
1381 break;
1382 case NC_INT64:
1383 for (dp = (double *)src, lip = dest; count < len; count++)
1384 {
1385 if (*dp > X_INT64_MAX || *dp < X_INT64_MIN)
1386 (*range_error)++;
1387 *lip++ = *dp++;
1388 }
1389 break;
1390 case NC_UINT64:
1391 for (dp = (double *)src, lip = dest; count < len; count++)
1392 {
1393 if (*dp > X_UINT64_MAX || *dp < 0)
1394 (*range_error)++;
1395 *lip++ = *dp++;
1396 }
1397 break;
1398 case NC_FLOAT:
1399 for (dp = (double *)src, fp = dest; count < len; count++)
1400 {
1401 if (isgreater(*dp, X_FLOAT_MAX) || isless(*dp, X_FLOAT_MIN))
1402 (*range_error)++;
1403 *fp++ = *dp++;
1404 }
1405 break;
1406 case NC_DOUBLE:
1407 for (dp = (double *)src, dp1 = dest; count < len; count++)
1408 *dp1++ = *dp++;
1409 break;
1410 default:
1411 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1412 __func__, src_type, dest_type));
1413 return NC_EBADTYPE;
1414 }
1415 break;
1416
1417 default:
1418 LOG((0, "%s: unexpected src type. src_type %d, dest_type %d",
1419 __func__, src_type, dest_type));
1420 return NC_EBADTYPE;
1421 }
1422
1423 /* If quantize is in use, determine masks, copy the data, do the
1424 * quantization. */
1425 if (quantize_mode == NC_QUANTIZE_BITGROOM)
1426 {
1427 if (dest_type == NC_FLOAT)
1428 {
1429 /* BitGroom: alternately shave and set LSBs */
1430 op1.fp = (float *)dest;
1431 u32_ptr = op1.ui32p;
1432 for (idx = 0L; idx < len; idx += 2L)
1433 if (op1.fp[idx] != mss_val_cmp_flt)
1434 u32_ptr[idx] &= msk_f32_u32_zro;
1435 for (idx = 1L; idx < len; idx += 2L)
1436 if (op1.fp[idx] != mss_val_cmp_flt && u32_ptr[idx] != 0U) /* Never quantize upwards floating point values of zero */
1437 u32_ptr[idx] |= msk_f32_u32_one;
1438 }
1439 else
1440 {
1441 /* BitGroom: alternately shave and set LSBs. */
1442 op1.dp = (double *)dest;
1443 u64_ptr = op1.ui64p;
1444 for (idx = 0L; idx < len; idx += 2L)
1445 if (op1.dp[idx] != mss_val_cmp_dbl)
1446 u64_ptr[idx] &= msk_f64_u64_zro;
1447 for (idx = 1L; idx < len; idx += 2L)
1448 if (op1.dp[idx] != mss_val_cmp_dbl && u64_ptr[idx] != 0ULL) /* Never quantize upwards floating point values of zero */
1449 u64_ptr[idx] |= msk_f64_u64_one;
1450 }
1451 } /* endif BitGroom */
1452
1453 if (quantize_mode == NC_QUANTIZE_BITROUND)
1454 {
1455 if (dest_type == NC_FLOAT)
1456 {
1457 /* BitRound: Quantize to user-specified NSB with IEEE-rounding */
1458 op1.fp = (float *)dest;
1459 u32_ptr = op1.ui32p;
1460 for (idx = 0L; idx < len; idx++){
1461 if (op1.fp[idx] != mss_val_cmp_flt){
1462 u32_ptr[idx] += msk_f32_u32_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1463 u32_ptr[idx] &= msk_f32_u32_zro; /* Shave it */
1464 }
1465 }
1466 }
1467 else
1468 {
1469 /* BitRound: Quantize to user-specified NSB with IEEE-rounding */
1470 op1.dp = (double *)dest;
1471 u64_ptr = op1.ui64p;
1472 for (idx = 0L; idx < len; idx++){
1473 if (op1.dp[idx] != mss_val_cmp_dbl){
1474 u64_ptr[idx] += msk_f64_u64_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1475 u64_ptr[idx] &= msk_f64_u64_zro; /* Shave it */
1476 }
1477 }
1478 }
1479 } /* endif BitRound */
1480
1481 if (quantize_mode == NC_QUANTIZE_GRANULARBR)
1482 {
1483 if (dest_type == NC_FLOAT)
1484 {
1485 /* Granular BitRound */
1486 op1.fp = (float *)dest;
1487 u32_ptr = op1.ui32p;
1488 for (idx = 0L; idx < len; idx++)
1489 {
1490
1491 if((val = op1.fp[idx]) != mss_val_cmp_flt && u32_ptr[idx] != 0U)
1492 {
1493 mnt = frexp(val, &xpn_bs2); /* DGG19 p. 4102 (8) */
1494 mnt_fabs = fabs(mnt);
1495 mnt_log10_fabs = log10(mnt_fabs);
1496 /* 20211003 Continuous determination of dgt_nbr improves CR by ~10% */
1497 dgt_nbr = (int)floor(xpn_bs2 * dgt_per_bit + mnt_log10_fabs) + 1; /* DGG19 p. 4102 (8.67) */
1498 qnt_pwr = (int)floor(bit_per_dgt * (dgt_nbr - nsd)); /* DGG19 p. 4101 (7) */
1499 prc_bnr_xpl_rqr = mnt_fabs == 0.0 ? 0 : abs((int)floor(xpn_bs2 - bit_per_dgt*mnt_log10_fabs) - qnt_pwr); /* Protect against mnt = -0.0 */
1500 prc_bnr_xpl_rqr--; /* 20211003 Reduce formula result by 1 bit: Passes all tests, improves CR by ~10% */
1501
1502 bit_xpl_nbr_zro = BIT_XPL_NBR_SGN_FLT - prc_bnr_xpl_rqr;
1503 msk_f32_u32_zro = 0u; /* Zero all bits */
1504 msk_f32_u32_zro = ~msk_f32_u32_zro; /* Turn all bits to ones */
1505 /* Bit Shave mask for AND: Left shift zeros into bits to be rounded, leave ones in untouched bits */
1506 msk_f32_u32_zro <<= bit_xpl_nbr_zro;
1507 /* Bit Set mask for OR: Put ones into bits to be set, zeros in untouched bits */
1508 msk_f32_u32_one = ~msk_f32_u32_zro;
1509 msk_f32_u32_hshv = msk_f32_u32_one & (msk_f32_u32_zro >> 1); /* Set one bit: the MSB of LSBs */
1510 u32_ptr[idx] += msk_f32_u32_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1511 u32_ptr[idx] &= msk_f32_u32_zro; /* Shave it */
1512
1513 } /* !mss_val_cmp_flt */
1514
1515 }
1516 }
1517 else
1518 {
1519 /* Granular BitRound */
1520 op1.dp = (double *)dest;
1521 u64_ptr = op1.ui64p;
1522 for (idx = 0L; idx < len; idx++)
1523 {
1524
1525 if((val = op1.dp[idx]) != mss_val_cmp_dbl && u64_ptr[idx] != 0ULL)
1526 {
1527 mnt = frexp(val, &xpn_bs2); /* DGG19 p. 4102 (8) */
1528 mnt_fabs = fabs(mnt);
1529 mnt_log10_fabs = log10(mnt_fabs);
1530 /* 20211003 Continuous determination of dgt_nbr improves CR by ~10% */
1531 dgt_nbr = (int)floor(xpn_bs2 * dgt_per_bit + mnt_log10_fabs) + 1; /* DGG19 p. 4102 (8.67) */
1532 qnt_pwr = (int)floor(bit_per_dgt * (dgt_nbr - nsd)); /* DGG19 p. 4101 (7) */
1533 prc_bnr_xpl_rqr = mnt_fabs == 0.0 ? 0 : abs((int)floor(xpn_bs2 - bit_per_dgt*mnt_log10_fabs) - qnt_pwr); /* Protect against mnt = -0.0 */
1534 prc_bnr_xpl_rqr--; /* 20211003 Reduce formula result by 1 bit: Passes all tests, improves CR by ~10% */
1535
1536 bit_xpl_nbr_zro = BIT_XPL_NBR_SGN_DBL - prc_bnr_xpl_rqr;
1537 msk_f64_u64_zro = 0ull; /* Zero all bits */
1538 msk_f64_u64_zro = ~msk_f64_u64_zro; /* Turn all bits to ones */
1539 /* Bit Shave mask for AND: Left shift zeros into bits to be rounded, leave ones in untouched bits */
1540 msk_f64_u64_zro <<= bit_xpl_nbr_zro;
1541 /* Bit Set mask for OR: Put ones into bits to be set, zeros in untouched bits */
1542 msk_f64_u64_one = ~msk_f64_u64_zro;
1543 msk_f64_u64_hshv = msk_f64_u64_one & (msk_f64_u64_zro >> 1); /* Set one bit: the MSB of LSBs */
1544 u64_ptr[idx] += msk_f64_u64_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1545 u64_ptr[idx] &= msk_f64_u64_zro; /* Shave it */
1546
1547 } /* !mss_val_cmp_dbl */
1548
1549 }
1550 }
1551 } /* endif GranularBR */
1552
1553 return NC_NOERR;
1554}
1555
1567int
1568nc4_get_fill_value(NC_FILE_INFO_T *h5, NC_VAR_INFO_T *var, void **fillp)
1569{
1570 size_t size;
1571 int retval;
1572
1573 /* Find out how much space we need for this type's fill value. */
1574 if (var->type_info->nc_type_class == NC_VLEN)
1575 size = sizeof(nc_vlen_t);
1576 else if (var->type_info->nc_type_class == NC_STRING)
1577 size = sizeof(char *);
1578 else
1579 {
1580 if ((retval = nc4_get_typelen_mem(h5, var->type_info->hdr.id, &size)))
1581 return retval;
1582 }
1583 assert(size);
1584
1585 /* Allocate the space. */
1586 if (!((*fillp) = calloc(1, size)))
1587 return NC_ENOMEM;
1588
1589 /* If the user has set a fill_value for this var, use, otherwise
1590 * find the default fill value. */
1591 if (var->fill_value)
1592 {
1593 LOG((4, "Found a fill value for var %s", var->hdr.name));
1594 if (var->type_info->nc_type_class == NC_VLEN)
1595 {
1596 nc_vlen_t *in_vlen = (nc_vlen_t *)(var->fill_value), *fv_vlen = (nc_vlen_t *)(*fillp);
1597 size_t basetypesize = 0;
1598
1599 if((retval=nc4_get_typelen_mem(h5, var->type_info->u.v.base_nc_typeid, &basetypesize)))
1600 return retval;
1601
1602 fv_vlen->len = in_vlen->len;
1603 if (!(fv_vlen->p = malloc(basetypesize * in_vlen->len)))
1604 {
1605 free(*fillp);
1606 *fillp = NULL;
1607 return NC_ENOMEM;
1608 }
1609 memcpy(fv_vlen->p, in_vlen->p, in_vlen->len * basetypesize);
1610 }
1611 else if (var->type_info->nc_type_class == NC_STRING)
1612 {
1613 if (*(char **)var->fill_value)
1614 if (!(**(char ***)fillp = strdup(*(char **)var->fill_value)))
1615 {
1616 free(*fillp);
1617 *fillp = NULL;
1618 return NC_ENOMEM;
1619 }
1620 }
1621 else
1622 memcpy((*fillp), var->fill_value, size);
1623 }
1624 else
1625 {
1626 if (nc4_get_default_fill_value(var->type_info, *fillp))
1627 {
1628 /* Note: release memory, but don't return error on failure */
1629 free(*fillp);
1630 *fillp = NULL;
1631 }
1632 }
1633
1634 return NC_NOERR;
1635}
1636
1649int
1650nc4_get_typelen_mem(NC_FILE_INFO_T *h5, nc_type xtype, size_t *len)
1651{
1652 NC_TYPE_INFO_T *type;
1653 int retval;
1654
1655 LOG((4, "%s xtype: %d", __func__, xtype));
1656 assert(len);
1657
1658 /* If this is an atomic type, the answer is easy. */
1659 switch (xtype)
1660 {
1661 case NC_BYTE:
1662 case NC_CHAR:
1663 case NC_UBYTE:
1664 *len = sizeof(char);
1665 return NC_NOERR;
1666 case NC_SHORT:
1667 case NC_USHORT:
1668 *len = sizeof(short);
1669 return NC_NOERR;
1670 case NC_INT:
1671 case NC_UINT:
1672 *len = sizeof(int);
1673 return NC_NOERR;
1674 case NC_FLOAT:
1675 *len = sizeof(float);
1676 return NC_NOERR;
1677 case NC_DOUBLE:
1678 *len = sizeof(double);
1679 return NC_NOERR;
1680 case NC_INT64:
1681 case NC_UINT64:
1682 *len = sizeof(long long);
1683 return NC_NOERR;
1684 case NC_STRING:
1685 *len = sizeof(char *);
1686 return NC_NOERR;
1687 }
1688
1689 /* See if var is compound type. */
1690 if ((retval = nc4_find_type(h5, xtype, &type)))
1691 return retval;
1692
1693 if (!type)
1694 return NC_EBADTYPE;
1695
1696 *len = type->size;
1697
1698 LOG((5, "type->size: %d", type->size));
1699
1700 return NC_NOERR;
1701}
1702
1703
1717int
1718nc4_check_chunksizes(NC_GRP_INFO_T *grp, NC_VAR_INFO_T *var, const size_t *chunksizes)
1719{
1720 double dprod;
1721 size_t type_len;
1722 int d;
1723 int retval;
1724
1725 if ((retval = nc4_get_typelen_mem(grp->nc4_info, var->type_info->hdr.id, &type_len)))
1726 return retval;
1727 if (var->type_info->nc_type_class == NC_VLEN)
1728 dprod = (double)sizeof(nc_vlen_t);
1729 else
1730 dprod = (double)type_len;
1731 for (d = 0; d < var->ndims; d++)
1732 dprod *= (double)chunksizes[d];
1733
1734 if (dprod > (double) NC_MAX_UINT)
1735 return NC_EBADCHUNK;
1736
1737 return NC_NOERR;
1738}
1739
1751int
1752nc4_find_default_chunksizes2(NC_GRP_INFO_T *grp, NC_VAR_INFO_T *var)
1753{
1754 int d;
1755 size_t type_size;
1756 float num_values = 1, num_unlim = 0;
1757 int retval;
1758 size_t suggested_size;
1759#ifdef LOGGING
1760 double total_chunk_size;
1761#endif
1762
1763 if (var->type_info->nc_type_class == NC_STRING)
1764 type_size = sizeof(char *);
1765 else
1766 type_size = var->type_info->size;
1767
1768#ifdef LOGGING
1769 /* Later this will become the total number of bytes in the default
1770 * chunk. */
1771 total_chunk_size = (double) type_size;
1772#endif
1773
1774 if(var->chunksizes == NULL) {
1775 if((var->chunksizes = calloc(1,sizeof(size_t)*var->ndims)) == NULL)
1776 return NC_ENOMEM;
1777 }
1778
1779 /* How many values in the variable (or one record, if there are
1780 * unlimited dimensions). */
1781 for (d = 0; d < var->ndims; d++)
1782 {
1783 assert(var->dim[d]);
1784 if (! var->dim[d]->unlimited)
1785 num_values *= (float)var->dim[d]->len;
1786 else {
1787 num_unlim++;
1788 var->chunksizes[d] = 1; /* overwritten below, if all dims are unlimited */
1789 }
1790 }
1791 /* Special case to avoid 1D vars with unlim dim taking huge amount
1792 of space (DEFAULT_CHUNK_SIZE bytes). Instead we limit to about
1793 4KB */
1794 if (var->ndims == 1 && num_unlim == 1) {
1795 if (DEFAULT_CHUNK_SIZE / type_size <= 0)
1796 suggested_size = 1;
1797 else if (DEFAULT_CHUNK_SIZE / type_size > DEFAULT_1D_UNLIM_SIZE)
1798 suggested_size = DEFAULT_1D_UNLIM_SIZE;
1799 else
1800 suggested_size = DEFAULT_CHUNK_SIZE / type_size;
1801 var->chunksizes[0] = suggested_size / type_size;
1802 LOG((4, "%s: name %s dim %d DEFAULT_CHUNK_SIZE %d num_values %f type_size %d "
1803 "chunksize %ld", __func__, var->hdr.name, d, DEFAULT_CHUNK_SIZE, num_values, type_size, var->chunksizes[0]));
1804 }
1805 if (var->ndims > 1 && var->ndims == num_unlim) { /* all dims unlimited */
1806 suggested_size = pow((double)DEFAULT_CHUNK_SIZE/type_size, 1.0/(double)(var->ndims));
1807 for (d = 0; d < var->ndims; d++)
1808 {
1809 var->chunksizes[d] = suggested_size ? suggested_size : 1;
1810 LOG((4, "%s: name %s dim %d DEFAULT_CHUNK_SIZE %d num_values %f type_size %d "
1811 "chunksize %ld", __func__, var->hdr.name, d, DEFAULT_CHUNK_SIZE, num_values, type_size, var->chunksizes[d]));
1812 }
1813 }
1814
1815 /* Pick a chunk length for each dimension, if one has not already
1816 * been picked above. */
1817 for (d = 0; d < var->ndims; d++)
1818 if (!var->chunksizes[d])
1819 {
1820 suggested_size = (pow((double)DEFAULT_CHUNK_SIZE/(num_values * type_size),
1821 1.0/(double)(var->ndims - num_unlim)) * var->dim[d]->len - .5);
1822 if (suggested_size > var->dim[d]->len)
1823 suggested_size = var->dim[d]->len;
1824 var->chunksizes[d] = suggested_size ? suggested_size : 1;
1825 LOG((4, "%s: name %s dim %d DEFAULT_CHUNK_SIZE %d num_values %f type_size %d "
1826 "chunksize %ld", __func__, var->hdr.name, d, DEFAULT_CHUNK_SIZE, num_values, type_size, var->chunksizes[d]));
1827 }
1828
1829#ifdef LOGGING
1830 /* Find total chunk size. */
1831 for (d = 0; d < var->ndims; d++)
1832 total_chunk_size *= (double) var->chunksizes[d];
1833 LOG((4, "total_chunk_size %f", total_chunk_size));
1834#endif
1835
1836 /* But did this result in a chunk that is too big? */
1837 retval = nc4_check_chunksizes(grp, var, var->chunksizes);
1838 if (retval)
1839 {
1840 /* Other error? */
1841 if (retval != NC_EBADCHUNK)
1842 return retval;
1843
1844 /* Chunk is too big! Reduce each dimension by half and try again. */
1845 for ( ; retval == NC_EBADCHUNK; retval = nc4_check_chunksizes(grp, var, var->chunksizes))
1846 for (d = 0; d < var->ndims; d++)
1847 var->chunksizes[d] = var->chunksizes[d]/2 ? var->chunksizes[d]/2 : 1;
1848 }
1849
1850 /* Do we have any big data overhangs? They can be dangerous to
1851 * babies, the elderly, or confused campers who have had too much
1852 * beer. */
1853 for (d = 0; d < var->ndims; d++)
1854 {
1855 size_t num_chunks;
1856 size_t overhang;
1857 assert(var->chunksizes[d] > 0);
1858 num_chunks = (var->dim[d]->len + var->chunksizes[d] - 1) / var->chunksizes[d];
1859 if(num_chunks > 0) {
1860 overhang = (num_chunks * var->chunksizes[d]) - var->dim[d]->len;
1861 var->chunksizes[d] -= overhang / num_chunks;
1862 }
1863 }
1864
1865
1866 return NC_NOERR;
1867}
1868
1880int
1881nc4_get_default_fill_value(NC_TYPE_INFO_T* tinfo, void *fill_value)
1882{
1883 if(tinfo->hdr.id > NC_NAT && tinfo->hdr.id <= NC_MAX_ATOMIC_TYPE)
1884 return nc4_get_default_atomic_fill_value(tinfo->hdr.id,fill_value);
1885#ifdef USE_NETCDF4
1886 switch(tinfo->nc_type_class) {
1887 case NC_ENUM:
1888 return nc4_get_default_atomic_fill_value(tinfo->u.e.base_nc_typeid,fill_value);
1889 case NC_OPAQUE:
1890 case NC_VLEN:
1891 case NC_COMPOUND:
1892 if(fill_value)
1893 memset(fill_value,0,tinfo->size);
1894 break;
1895 default: return NC_EBADTYPE;
1896 }
1897#endif
1898 return NC_NOERR;
1899}
1900
1912int
1913nc4_get_default_atomic_fill_value(nc_type xtype, void *fill_value)
1914{
1915 switch (xtype)
1916 {
1917 case NC_CHAR:
1918 *(char *)fill_value = NC_FILL_CHAR;
1919 break;
1920
1921 case NC_STRING:
1922 *(char **)fill_value = strdup(NC_FILL_STRING);
1923 break;
1924
1925 case NC_BYTE:
1926 *(signed char *)fill_value = NC_FILL_BYTE;
1927 break;
1928
1929 case NC_SHORT:
1930 *(short *)fill_value = NC_FILL_SHORT;
1931 break;
1932
1933 case NC_INT:
1934 *(int *)fill_value = NC_FILL_INT;
1935 break;
1936
1937 case NC_UBYTE:
1938 *(unsigned char *)fill_value = NC_FILL_UBYTE;
1939 break;
1940
1941 case NC_USHORT:
1942 *(unsigned short *)fill_value = NC_FILL_USHORT;
1943 break;
1944
1945 case NC_UINT:
1946 *(unsigned int *)fill_value = NC_FILL_UINT;
1947 break;
1948
1949 case NC_INT64:
1950 *(long long *)fill_value = NC_FILL_INT64;
1951 break;
1952
1953 case NC_UINT64:
1954 *(unsigned long long *)fill_value = NC_FILL_UINT64;
1955 break;
1956
1957 case NC_FLOAT:
1958 *(float *)fill_value = NC_FILL_FLOAT;
1959 break;
1960
1961 case NC_DOUBLE:
1962 *(double *)fill_value = NC_FILL_DOUBLE;
1963 break;
1964
1965 default:
1966 return NC_EINVAL;
1967 }
1968 return NC_NOERR;
1969}
1970
EXTERNL int nc_inq_var_filter_info(int ncid, int varid, unsigned int id, size_t *nparams, unsigned int *params)
Find the the param info about filter (if any) associated with a variable and with specified id.
Definition dfilter.c:97
#define M_LN10
log_e 10
Definition nc4var.c:28
#define BIT_XPL_NBR_SGN_DBL
Used in quantize code.
Definition nc4var.c:44
#define M_LN2
log_e 2
Definition nc4var.c:31
#define BIT_XPL_NBR_SGN_FLT
Used in quantize code.
Definition nc4var.c:38
#define NC_EBADTYPE
Not a netcdf data type.
Definition netcdf.h:410
#define NC_UINT
unsigned 4-byte int
Definition netcdf.h:44
#define NC_FILL_BYTE
Default fill value.
Definition netcdf.h:67
void * p
Pointer to VL data.
Definition netcdf.h:748
#define NC_EFILTER
Filter operation failed.
Definition netcdf.h:513
#define NC_INT
signed 4 byte integer
Definition netcdf.h:38
#define NC_FILL_INT
Default fill value.
Definition netcdf.h:70
#define NC_BYTE
signed 1 byte integer
Definition netcdf.h:35
#define NC_QUANTIZE_BITROUND
Use BitRound quantization.
Definition netcdf.h:338
size_t len
Length of VL data (in base type units)
Definition netcdf.h:747
#define NC_VLEN
vlen (variable-length) types
Definition netcdf.h:53
#define NC_FILL_UBYTE
Default fill value.
Definition netcdf.h:73
#define NC_NAT
Not A Type.
Definition netcdf.h:34
#define NC_MAX_UINT
Max or min values for a type.
Definition netcdf.h:102
#define NC_DOUBLE
double precision floating point number
Definition netcdf.h:41
#define NC_UBYTE
unsigned 1 byte int
Definition netcdf.h:42
#define NC_FLOAT
single precision floating point number
Definition netcdf.h:40
#define NC_ENOMEM
Memory allocation (malloc) failure.
Definition netcdf.h:448
#define NC_MAX_INT
Max or min values for a type.
Definition netcdf.h:94
#define NC_COMPOUND
compound types
Definition netcdf.h:56
EXTERNL int nc_copy_data(int ncid, nc_type xtypeid, const void *memory, size_t count, void *copy)
Copy vector of arbitrary type instances.
#define NC_CHUNKED
In HDF5 files you can set storage for each variable to be either contiguous or chunked,...
Definition netcdf.h:304
#define NC_SHORT
signed 2 byte integer
Definition netcdf.h:37
#define NC_QUANTIZE_GRANULARBR
Use Granular BitRound quantization.
Definition netcdf.h:337
#define NC_FILL_UINT64
Default fill value.
Definition netcdf.h:77
#define NC_QUANTIZE_BITGROOM
Use BitGroom quantization.
Definition netcdf.h:336
#define NC_ENUM
enum types
Definition netcdf.h:55
#define NC_INT64
signed 8-byte int
Definition netcdf.h:45
#define NC_GLOBAL
Attribute id to put/get a global attribute.
Definition netcdf.h:254
#define NC_FILL_UINT
Default fill value.
Definition netcdf.h:75
#define NC_FILL_CHAR
Default fill value.
Definition netcdf.h:68
#define NC_UINT64
unsigned 8-byte int
Definition netcdf.h:46
#define NC_ENOTVAR
Variable not found.
Definition netcdf.h:422
#define NC_EINVAL
Invalid Argument.
Definition netcdf.h:378
#define NC_MAX_NAME
Maximum for classic library.
Definition netcdf.h:281
#define NC_NOERR
No Error.
Definition netcdf.h:368
#define NC_FILL_USHORT
Default fill value.
Definition netcdf.h:74
#define NC_FILL_SHORT
Default fill value.
Definition netcdf.h:69
#define NC_FILL_INT64
Default fill value.
Definition netcdf.h:76
#define NC_NOQUANTIZE
No quantization in use.
Definition netcdf.h:335
#define NC_USHORT
unsigned 2-byte int
Definition netcdf.h:43
#define NC_OPAQUE
opaque types
Definition netcdf.h:54
#define NC_ENOPAR
Parallel operation on file opened for non-parallel access.
Definition netcdf.h:494
#define NC_STRING
string
Definition netcdf.h:47
#define NC_EBADCHUNK
Bad chunksize.
Definition netcdf.h:507
#define NC_FILL_FLOAT
Default fill value.
Definition netcdf.h:71
#define NC_ENOFILTER
Filter not defined on variable.
Definition netcdf.h:517
#define NC_FILL_DOUBLE
Default fill value.
Definition netcdf.h:72
#define NC_CHAR
ISO/ASCII character.
Definition netcdf.h:36
#define NC_ERANGE
Math result not representable.
Definition netcdf.h:447
#define NC_FILL_STRING
Default fill value.
Definition netcdf.h:78
int nc_type
The nc_type type is just an int.
Definition netcdf.h:25
This is the type of arrays of vlens.
Definition netcdf.h:746
#define NC_COLLECTIVE
Use with nc_var_par_access() to set parallel access to collective.
Definition netcdf_par.h:28
#define NC_INDEPENDENT
Use with nc_var_par_access() to set parallel access to independent.
Definition netcdf_par.h:26