NetCDF 4.10.0
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,(size_t)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 = (int)(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, (size_t)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 = (int)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 {
267 int xtype = var->type_info->hdr.id;
268 if((retval = NC_copy_data(h5->controller,xtype,var->fill_value,1,fill_valuep))) return retval;
269 }
270 else
271 {
272 if ((retval = nc4_get_default_fill_value(var->type_info, fill_valuep)))
273 return retval;
274 }
275 }
276
277 /* Does the user want the endianness of this variable? */
278 if (endiannessp)
279 *endiannessp = var->endianness;
280
281 return NC_NOERR;
282}
283
300int
301nc_inq_var_chunking_ints(int ncid, int varid, int *storagep, int *chunksizesp)
302{
303 NC_VAR_INFO_T *var;
304 size_t *cs = NULL;
305 int i, retval;
306
307 /* Get pointer to the var. */
308 if ((retval = nc4_find_grp_h5_var(ncid, varid, NULL, NULL, &var)))
309 return retval;
310 assert(var);
311
312 /* Allocate space for the size_t copy of the chunksizes array. */
313 if (var->ndims)
314 if (!(cs = malloc(var->ndims * sizeof(size_t))))
315 return NC_ENOMEM;
316
317 /* Call the netcdf-4 version directly. */
318 retval = NC4_inq_var_all(ncid, varid, NULL, NULL, NULL, NULL, NULL,
319 NULL, NULL, NULL, NULL, storagep, cs, NULL,
320 NULL, NULL, NULL, NULL, NULL);
321
322 /* Copy from size_t array. */
323 if (!retval && chunksizesp && var->storage == NC_CHUNKED)
324 {
325 for (i = 0; i < var->ndims; i++)
326 {
327 chunksizesp[i] = (int)cs[i];
328 if (cs[i] > NC_MAX_INT)
329 retval = NC_ERANGE;
330 }
331 }
332
333 if (var->ndims)
334 free(cs);
335 return retval;
336}
337
350int
351NC4_inq_varid(int ncid, const char *name, int *varidp)
352{
353 NC *nc;
354 NC_GRP_INFO_T *grp;
355 NC_VAR_INFO_T *var;
356 char norm_name[NC_MAX_NAME + 1];
357 int retval;
358
359 if (!name)
360 return NC_EINVAL;
361 if (!varidp)
362 return NC_NOERR;
363
364 LOG((2, "%s: ncid 0x%x name %s", __func__, ncid, name));
365
366 /* Find info for this file and group, and set pointer to each. */
367 if ((retval = nc4_find_nc_grp_h5(ncid, &nc, &grp, NULL)))
368 return retval;
369
370 /* Normalize name. */
371 if ((retval = nc4_normalize_name(name, norm_name)))
372 return retval;
373
374 /* Find var of this name. */
375 var = (NC_VAR_INFO_T*)ncindexlookup(grp->vars,norm_name);
376 if(var)
377 {
378 *varidp = var->hdr.id;
379 return NC_NOERR;
380 }
381 return NC_ENOTVAR;
382}
383
402int
403NC4_var_par_access(int ncid, int varid, int par_access)
404{
405#ifndef USE_PARALLEL4
406 NC_UNUSED(ncid);
407 NC_UNUSED(varid);
408 NC_UNUSED(par_access);
409 return NC_ENOPAR;
410#else
411 NC *nc;
412 NC_GRP_INFO_T *grp;
413 NC_FILE_INFO_T *h5;
414 NC_VAR_INFO_T *var;
415 int retval;
416
417 LOG((1, "%s: ncid 0x%x varid %d par_access %d", __func__, ncid,
418 varid, par_access));
419
420 if (par_access != NC_INDEPENDENT && par_access != NC_COLLECTIVE)
421 return NC_EINVAL;
422
423 /* Find info for this file and group, and set pointer to each. */
424 if ((retval = nc4_find_nc_grp_h5(ncid, &nc, &grp, &h5)))
425 return retval;
426
427 /* This function only for files opened with nc_open_par or nc_create_par. */
428 if (!h5->parallel)
429 return NC_ENOPAR;
430
431 /* Find the var, and set its preference. */
432 var = (NC_VAR_INFO_T*)ncindexith(grp->vars,varid);
433 if (!var) return NC_ENOTVAR;
434 assert(var->hdr.id == varid);
435
436 /* If zlib, shuffle, or fletcher32 filters are in use, then access
437 * must be collective. Fail an attempt to set such a variable to
438 * independent access. */
439 if (nclistlength((NClist*)var->filters) > 0 &&
440 par_access == NC_INDEPENDENT)
441 return NC_EINVAL;
442
443 if (par_access)
444 var->parallel_access = NC_COLLECTIVE;
445 else
446 var->parallel_access = NC_INDEPENDENT;
447 return NC_NOERR;
448#endif /* USE_PARALLEL4 */
449}
450
483int
484nc4_convert_type(const void *src, void *dest, const nc_type src_type,
485 const nc_type dest_type, const size_t len, int *range_error,
486 const void *fill_value, int strict_nc3, int quantize_mode,
487 int nsd)
488{
489 /* These vars are used with quantize feature. */
490 const double bit_per_dgt = M_LN10 / M_LN2; /* 3.32 [frc] Bits per decimal digit of precision = log2(10) */
491 const double dgt_per_bit= M_LN2 / M_LN10; /* 0.301 [frc] Decimal digits per bit of precision = log10(2) */
492 double mnt; /* [frc] Mantissa, 0.5 <= mnt < 1.0 */
493 double mnt_fabs; /* [frc] fabs(mantissa) */
494 double mnt_log10_fabs; /* [frc] log10(fabs(mantissa))) */
495 double mss_val_cmp_dbl; /* Missing value for comparison to double precision values */
496 double val_dbl; /* [frc] Copy of input value to avoid indirection */
497 float mss_val_cmp_flt; /* Missing value for comparison to single precision values */
498 float val_flt; /* [frc] Copy of input value to avoid indirection */
499 int bit_xpl_nbr_zro; /* [nbr] Number of explicit bits to zero */
500 int dgt_nbr; /* [nbr] Number of digits before decimal point */
501 int qnt_pwr; /* [nbr] Power of two in quantization mask: qnt_msk = 2^qnt_pwr */
502 int xpn_bs2; /* [nbr] Binary exponent xpn_bs2 in val = sign(val) * 2^xpn_bs2 * mnt, 0.5 < mnt <= 1.0 */
503 size_t idx;
504 unsigned int *u32_ptr;
505 unsigned int msk_f32_u32_zro;
506 unsigned int msk_f32_u32_one;
507 unsigned int msk_f32_u32_hshv;
508 unsigned long long int *u64_ptr;
509 unsigned long long int msk_f64_u64_zro;
510 unsigned long long int msk_f64_u64_one;
511 unsigned long long int msk_f64_u64_hshv;
512 unsigned short prc_bnr_xpl_rqr; /* [nbr] Explicitly represented binary digits required to retain */
513 ptr_unn op1; /* I/O [frc] Values to quantize */
514
515 char *cp, *cp1;
516 float *fp, *fp1;
517 double *dp, *dp1;
518 int *ip, *ip1;
519 short *sp, *sp1;
520 signed char *bp, *bp1;
521 unsigned char *ubp, *ubp1;
522 unsigned short *usp, *usp1;
523 unsigned int *uip, *uip1;
524 long long *lip, *lip1;
525 unsigned long long *ulip, *ulip1;
526 size_t count = 0;
527
528 *range_error = 0;
529 LOG((3, "%s: len %d src_type %d dest_type %d", __func__, len, src_type,
530 dest_type));
531
532 /* If quantize is in use, set up some values. Quantize can only be
533 * used when the destination type is NC_FLOAT or NC_DOUBLE. */
534 if (quantize_mode != NC_NOQUANTIZE)
535 {
536 assert(dest_type == NC_FLOAT || dest_type == NC_DOUBLE);
537
538 /* Parameters shared by all quantization codecs */
539 if (dest_type == NC_FLOAT)
540 {
541 /* Determine the fill value. */
542 if (fill_value)
543 mss_val_cmp_flt = *(float *)fill_value;
544 else
545 mss_val_cmp_flt = NC_FILL_FLOAT;
546
547 }
548 else
549 {
550
551 /* Determine the fill value. */
552 if (fill_value)
553 mss_val_cmp_dbl = *(double *)fill_value;
554 else
555 mss_val_cmp_dbl = NC_FILL_DOUBLE;
556
557 }
558
559 /* Set parameters used by BitGroom and BitRound here, outside value loop.
560 Equivalent parameters used by GranularBR are set inside value loop,
561 since keep bits and thus masks can change for every value. */
562 if (quantize_mode == NC_QUANTIZE_BITGROOM ||
563 quantize_mode == NC_QUANTIZE_BITROUND )
564 {
565
566 if (quantize_mode == NC_QUANTIZE_BITGROOM){
567
568 /* BitGroom interprets nsd as number of significant decimal digits
569 * Must convert that to number of significant bits to preserve
570 * How many bits to preserve? Being conservative, we round up the
571 * exact binary digits of precision. Add one because the first bit
572 * is implicit not explicit but corner cases prevent our taking
573 * advantage of this. */
574 prc_bnr_xpl_rqr = (unsigned short)ceil(nsd * bit_per_dgt) + 1;
575
576 }else if (quantize_mode == NC_QUANTIZE_BITROUND){
577
578 /* BitRound interprets nsd as number of significant binary digits (bits) */
579 prc_bnr_xpl_rqr = (unsigned short)nsd;
580
581 }
582
583 if (dest_type == NC_FLOAT)
584 {
585
586 bit_xpl_nbr_zro = BIT_XPL_NBR_SGN_FLT - prc_bnr_xpl_rqr;
587
588 /* Create mask */
589 msk_f32_u32_zro = 0U; /* Zero all bits */
590 msk_f32_u32_zro = ~msk_f32_u32_zro; /* Turn all bits to ones */
591
592 /* BitShave mask for AND: Left shift zeros into bits to be
593 * rounded, leave ones in untouched bits. */
594 msk_f32_u32_zro <<= bit_xpl_nbr_zro;
595
596 /* BitSet mask for OR: Put ones into bits to be set, zeros in
597 * untouched bits. */
598 msk_f32_u32_one = ~msk_f32_u32_zro;
599
600 /* BitRound mask for ADD: Set one bit: the MSB of LSBs */
601 msk_f32_u32_hshv=msk_f32_u32_one & (msk_f32_u32_zro >> 1);
602
603 }
604 else
605 {
606
607 bit_xpl_nbr_zro = BIT_XPL_NBR_SGN_DBL - prc_bnr_xpl_rqr;
608 /* Create mask. */
609 msk_f64_u64_zro = 0UL; /* Zero all bits. */
610 msk_f64_u64_zro = ~msk_f64_u64_zro; /* Turn all bits to ones. */
611
612 /* BitShave mask for AND: Left shift zeros into bits to be
613 * rounded, leave ones in untouched bits. */
614 msk_f64_u64_zro <<= bit_xpl_nbr_zro;
615
616 /* BitSet mask for OR: Put ones into bits to be set, zeros in
617 * untouched bits. */
618 msk_f64_u64_one =~ msk_f64_u64_zro;
619
620 /* BitRound mask for ADD: Set one bit: the MSB of LSBs */
621 msk_f64_u64_hshv = msk_f64_u64_one & (msk_f64_u64_zro >> 1);
622
623 }
624
625 }
626
627 } /* endif quantize */
628
629 /* OK, this is ugly. If you can think of anything better, I'm open
630 to suggestions!
631
632 Note that we don't use a default fill value for type
633 NC_BYTE. This is because Lord Voldemort cast a nofilleramous spell
634 at Harry Potter, but it bounced off his scar and hit the netcdf-4
635 code.
636 */
637 switch (src_type)
638 {
639 case NC_CHAR:
640 switch (dest_type)
641 {
642 case NC_CHAR:
643 for (cp = (char *)src, cp1 = dest; count < len; count++)
644 *cp1++ = *cp++;
645 break;
646 default:
647 LOG((0, "%s: Unknown destination type.", __func__));
648 }
649 break;
650
651 case NC_BYTE:
652 switch (dest_type)
653 {
654 case NC_BYTE:
655 for (bp = (signed char *)src, bp1 = dest; count < len; count++)
656 *bp1++ = *bp++;
657 break;
658 case NC_UBYTE:
659 for (bp = (signed char *)src, ubp = dest; count < len; count++)
660 {
661 if (*bp < 0)
662 (*range_error)++;
663 *ubp++ = (unsigned char)*bp++;
664 }
665 break;
666 case NC_SHORT:
667 for (bp = (signed char *)src, sp = dest; count < len; count++)
668 *sp++ = *bp++;
669 break;
670 case NC_USHORT:
671 for (bp = (signed char *)src, usp = dest; count < len; count++)
672 {
673 if (*bp < 0)
674 (*range_error)++;
675 *usp++ = (unsigned short)*bp++;
676 }
677 break;
678 case NC_INT:
679 for (bp = (signed char *)src, ip = dest; count < len; count++)
680 *ip++ = *bp++;
681 break;
682 case NC_UINT:
683 for (bp = (signed char *)src, uip = dest; count < len; count++)
684 {
685 if (*bp < 0)
686 (*range_error)++;
687 *uip++ = (unsigned int)*bp++;
688 }
689 break;
690 case NC_INT64:
691 for (bp = (signed char *)src, lip = dest; count < len; count++)
692 *lip++ = *bp++;
693 break;
694 case NC_UINT64:
695 for (bp = (signed char *)src, ulip = dest; count < len; count++)
696 {
697 if (*bp < 0)
698 (*range_error)++;
699 *ulip++ = (unsigned long long)*bp++;
700 }
701 break;
702 case NC_FLOAT:
703 for (bp = (signed char *)src, fp = dest; count < len; count++)
704 *fp++ = *bp++;
705 break;
706 case NC_DOUBLE:
707 for (bp = (signed char *)src, dp = dest; count < len; count++)
708 *dp++ = *bp++;
709 break;
710 default:
711 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
712 __func__, src_type, dest_type));
713 return NC_EBADTYPE;
714 }
715 break;
716
717 case NC_UBYTE:
718 switch (dest_type)
719 {
720 case NC_BYTE:
721 for (ubp = (unsigned char *)src, bp = dest; count < len; count++)
722 {
723 if (!strict_nc3 && *ubp > X_SCHAR_MAX)
724 (*range_error)++;
725 *bp++ = (signed char)*ubp++;
726 }
727 break;
728 case NC_SHORT:
729 for (ubp = (unsigned char *)src, sp = dest; count < len; count++)
730 *sp++ = *ubp++;
731 break;
732 case NC_UBYTE:
733 for (ubp = (unsigned char *)src, ubp1 = dest; count < len; count++)
734 *ubp1++ = *ubp++;
735 break;
736 case NC_USHORT:
737 for (ubp = (unsigned char *)src, usp = dest; count < len; count++)
738 *usp++ = *ubp++;
739 break;
740 case NC_INT:
741 for (ubp = (unsigned char *)src, ip = dest; count < len; count++)
742 *ip++ = *ubp++;
743 break;
744 case NC_UINT:
745 for (ubp = (unsigned char *)src, uip = dest; count < len; count++)
746 *uip++ = *ubp++;
747 break;
748 case NC_INT64:
749 for (ubp = (unsigned char *)src, lip = dest; count < len; count++)
750 *lip++ = *ubp++;
751 break;
752 case NC_UINT64:
753 for (ubp = (unsigned char *)src, ulip = dest; count < len; count++)
754 *ulip++ = *ubp++;
755 break;
756 case NC_FLOAT:
757 for (ubp = (unsigned char *)src, fp = dest; count < len; count++)
758 *fp++ = *ubp++;
759 break;
760 case NC_DOUBLE:
761 for (ubp = (unsigned char *)src, dp = dest; count < len; count++)
762 *dp++ = *ubp++;
763 break;
764 default:
765 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
766 __func__, src_type, dest_type));
767 return NC_EBADTYPE;
768 }
769 break;
770
771 case NC_SHORT:
772 switch (dest_type)
773 {
774 case NC_UBYTE:
775 for (sp = (short *)src, ubp = dest; count < len; count++)
776 {
777 if (*sp > X_UCHAR_MAX || *sp < 0)
778 (*range_error)++;
779 *ubp++ = (unsigned char)*sp++;
780 }
781 break;
782 case NC_BYTE:
783 for (sp = (short *)src, bp = dest; count < len; count++)
784 {
785 if (*sp > X_SCHAR_MAX || *sp < X_SCHAR_MIN)
786 (*range_error)++;
787 *bp++ = (signed char)*sp++;
788 }
789 break;
790 case NC_SHORT:
791 for (sp = (short *)src, sp1 = dest; count < len; count++)
792 *sp1++ = *sp++;
793 break;
794 case NC_USHORT:
795 for (sp = (short *)src, usp = dest; count < len; count++)
796 {
797 if (*sp < 0)
798 (*range_error)++;
799 *usp++ = (unsigned short)*sp++;
800 }
801 break;
802 case NC_INT:
803 for (sp = (short *)src, ip = dest; count < len; count++)
804 *ip++ = *sp++;
805 break;
806 case NC_UINT:
807 for (sp = (short *)src, uip = dest; count < len; count++)
808 {
809 if (*sp < 0)
810 (*range_error)++;
811 *uip++ = (unsigned int)*sp++;
812 }
813 break;
814 case NC_INT64:
815 for (sp = (short *)src, lip = dest; count < len; count++)
816 *lip++ = *sp++;
817 break;
818 case NC_UINT64:
819 for (sp = (short *)src, ulip = dest; count < len; count++)
820 {
821 if (*sp < 0)
822 (*range_error)++;
823 *ulip++ = (unsigned long long)*sp++;
824 }
825 break;
826 case NC_FLOAT:
827 for (sp = (short *)src, fp = dest; count < len; count++)
828 *fp++ = *sp++;
829 break;
830 case NC_DOUBLE:
831 for (sp = (short *)src, dp = dest; count < len; count++)
832 *dp++ = *sp++;
833 break;
834 default:
835 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
836 __func__, src_type, dest_type));
837 return NC_EBADTYPE;
838 }
839 break;
840
841 case NC_USHORT:
842 switch (dest_type)
843 {
844 case NC_UBYTE:
845 for (usp = (unsigned short *)src, ubp = dest; count < len; count++)
846 {
847 if (*usp > X_UCHAR_MAX)
848 (*range_error)++;
849 *ubp++ = (unsigned char)*usp++;
850 }
851 break;
852 case NC_BYTE:
853 for (usp = (unsigned short *)src, bp = dest; count < len; count++)
854 {
855 if (*usp > X_SCHAR_MAX)
856 (*range_error)++;
857 *bp++ = (signed char)*usp++;
858 }
859 break;
860 case NC_SHORT:
861 for (usp = (unsigned short *)src, sp = dest; count < len; count++)
862 {
863 if (*usp > X_SHORT_MAX)
864 (*range_error)++;
865 *sp++ = (signed short)*usp++;
866 }
867 break;
868 case NC_USHORT:
869 for (usp = (unsigned short *)src, usp1 = dest; count < len; count++)
870 *usp1++ = *usp++;
871 break;
872 case NC_INT:
873 for (usp = (unsigned short *)src, ip = dest; count < len; count++)
874 *ip++ = *usp++;
875 break;
876 case NC_UINT:
877 for (usp = (unsigned short *)src, uip = dest; count < len; count++)
878 *uip++ = *usp++;
879 break;
880 case NC_INT64:
881 for (usp = (unsigned short *)src, lip = dest; count < len; count++)
882 *lip++ = *usp++;
883 break;
884 case NC_UINT64:
885 for (usp = (unsigned short *)src, ulip = dest; count < len; count++)
886 *ulip++ = *usp++;
887 break;
888 case NC_FLOAT:
889 for (usp = (unsigned short *)src, fp = dest; count < len; count++)
890 *fp++ = *usp++;
891 break;
892 case NC_DOUBLE:
893 for (usp = (unsigned short *)src, dp = dest; count < len; count++)
894 *dp++ = *usp++;
895 break;
896 default:
897 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
898 __func__, src_type, dest_type));
899 return NC_EBADTYPE;
900 }
901 break;
902
903 case NC_INT:
904 switch (dest_type)
905 {
906 case NC_UBYTE:
907 for (ip = (int *)src, ubp = dest; count < len; count++)
908 {
909 if (*ip > X_UCHAR_MAX || *ip < 0)
910 (*range_error)++;
911 *ubp++ = (unsigned char)*ip++;
912 }
913 break;
914 case NC_BYTE:
915 for (ip = (int *)src, bp = dest; count < len; count++)
916 {
917 if (*ip > X_SCHAR_MAX || *ip < X_SCHAR_MIN)
918 (*range_error)++;
919 *bp++ = (signed char)*ip++;
920 }
921 break;
922 case NC_SHORT:
923 for (ip = (int *)src, sp = dest; count < len; count++)
924 {
925 if (*ip > X_SHORT_MAX || *ip < X_SHORT_MIN)
926 (*range_error)++;
927 *sp++ = (short)*ip++;
928 }
929 break;
930 case NC_USHORT:
931 for (ip = (int *)src, usp = dest; count < len; count++)
932 {
933 if (*ip > X_USHORT_MAX || *ip < 0)
934 (*range_error)++;
935 *usp++ = (unsigned short)*ip++;
936 }
937 break;
938 case NC_INT: /* src is int */
939 for (ip = (int *)src, ip1 = dest; count < len; count++)
940 {
941 if (*ip > X_INT_MAX || *ip < X_INT_MIN)
942 (*range_error)++;
943 *ip1++ = *ip++;
944 }
945 break;
946 case NC_UINT:
947 for (ip = (int *)src, uip = dest; count < len; count++)
948 {
949 if (*ip > X_UINT_MAX || *ip < 0)
950 (*range_error)++;
951 *uip++ = (unsigned int)*ip++;
952 }
953 break;
954 case NC_INT64:
955 for (ip = (int *)src, lip = dest; count < len; count++)
956 *lip++ = *ip++;
957 break;
958 case NC_UINT64:
959 for (ip = (int *)src, ulip = dest; count < len; count++)
960 {
961 if (*ip < 0)
962 (*range_error)++;
963 *ulip++ = (unsigned long long)*ip++;
964 }
965 break;
966 case NC_FLOAT:
967 for (ip = (int *)src, fp = dest; count < len; count++)
968 *fp++ = (float)*ip++;
969 break;
970 case NC_DOUBLE:
971 for (ip = (int *)src, dp = dest; count < len; count++)
972 *dp++ = (double)*ip++;
973 break;
974 default:
975 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
976 __func__, src_type, dest_type));
977 return NC_EBADTYPE;
978 }
979 break;
980
981 case NC_UINT:
982 switch (dest_type)
983 {
984 case NC_UBYTE:
985 for (uip = (unsigned int *)src, ubp = dest; count < len; count++)
986 {
987 if (*uip > X_UCHAR_MAX)
988 (*range_error)++;
989 *ubp++ = (unsigned char)*uip++;
990 }
991 break;
992 case NC_BYTE:
993 for (uip = (unsigned int *)src, bp = dest; count < len; count++)
994 {
995 if (*uip > X_SCHAR_MAX)
996 (*range_error)++;
997 *bp++ = (signed char)*uip++;
998 }
999 break;
1000 case NC_SHORT:
1001 for (uip = (unsigned int *)src, sp = dest; count < len; count++)
1002 {
1003 if (*uip > X_SHORT_MAX)
1004 (*range_error)++;
1005 *sp++ = (signed short)*uip++;
1006 }
1007 break;
1008 case NC_USHORT:
1009 for (uip = (unsigned int *)src, usp = dest; count < len; count++)
1010 {
1011 if (*uip > X_USHORT_MAX)
1012 (*range_error)++;
1013 *usp++ = (unsigned short)*uip++;
1014 }
1015 break;
1016 case NC_INT:
1017 for (uip = (unsigned int *)src, ip = dest; count < len; count++)
1018 {
1019 if (*uip > X_INT_MAX)
1020 (*range_error)++;
1021 *ip++ = (int)*uip++;
1022 }
1023 break;
1024 case NC_UINT:
1025 for (uip = (unsigned int *)src, uip1 = dest; count < len; count++)
1026 {
1027 if (*uip > X_UINT_MAX)
1028 (*range_error)++;
1029 *uip1++ = *uip++;
1030 }
1031 break;
1032 case NC_INT64:
1033 for (uip = (unsigned int *)src, lip = dest; count < len; count++)
1034 *lip++ = *uip++;
1035 break;
1036 case NC_UINT64:
1037 for (uip = (unsigned int *)src, ulip = dest; count < len; count++)
1038 *ulip++ = *uip++;
1039 break;
1040 case NC_FLOAT:
1041 for (uip = (unsigned int *)src, fp = dest; count < len; count++)
1042 *fp++ = (float)*uip++;
1043 break;
1044 case NC_DOUBLE:
1045 for (uip = (unsigned int *)src, dp = dest; count < len; count++)
1046 *dp++ = *uip++;
1047 break;
1048 default:
1049 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1050 __func__, src_type, dest_type));
1051 return NC_EBADTYPE;
1052 }
1053 break;
1054
1055 case NC_INT64:
1056 switch (dest_type)
1057 {
1058 case NC_UBYTE:
1059 for (lip = (long long *)src, ubp = dest; count < len; count++)
1060 {
1061 if (*lip > X_UCHAR_MAX || *lip < 0)
1062 (*range_error)++;
1063 *ubp++ = (unsigned char)*lip++;
1064 }
1065 break;
1066 case NC_BYTE:
1067 for (lip = (long long *)src, bp = dest; count < len; count++)
1068 {
1069 if (*lip > X_SCHAR_MAX || *lip < X_SCHAR_MIN)
1070 (*range_error)++;
1071 *bp++ = (signed char)*lip++;
1072 }
1073 break;
1074 case NC_SHORT:
1075 for (lip = (long long *)src, sp = dest; count < len; count++)
1076 {
1077 if (*lip > X_SHORT_MAX || *lip < X_SHORT_MIN)
1078 (*range_error)++;
1079 *sp++ = (short)*lip++;
1080 }
1081 break;
1082 case NC_USHORT:
1083 for (lip = (long long *)src, usp = dest; count < len; count++)
1084 {
1085 if (*lip > X_USHORT_MAX || *lip < 0)
1086 (*range_error)++;
1087 *usp++ = (unsigned short)*lip++;
1088 }
1089 break;
1090 case NC_UINT:
1091 for (lip = (long long *)src, uip = dest; count < len; count++)
1092 {
1093 if (*lip > X_UINT_MAX || *lip < 0)
1094 (*range_error)++;
1095 *uip++ = (unsigned int)*lip++;
1096 }
1097 break;
1098 case NC_INT:
1099 for (lip = (long long *)src, ip = dest; count < len; count++)
1100 {
1101 if (*lip > X_INT_MAX || *lip < X_INT_MIN)
1102 (*range_error)++;
1103 *ip++ = (int)*lip++;
1104 }
1105 break;
1106 case NC_INT64:
1107 for (lip = (long long *)src, lip1 = dest; count < len; count++)
1108 *lip1++ = *lip++;
1109 break;
1110 case NC_UINT64:
1111 for (lip = (long long *)src, ulip = dest; count < len; count++)
1112 {
1113 if (*lip < 0)
1114 (*range_error)++;
1115 *ulip++ = (unsigned long long)*lip++;
1116 }
1117 break;
1118 case NC_FLOAT:
1119 for (lip = (long long *)src, fp = dest; count < len; count++)
1120 *fp++ = (float)*lip++;
1121 break;
1122 case NC_DOUBLE:
1123 for (lip = (long long *)src, dp = dest; count < len; count++)
1124 *dp++ = (double)*lip++;
1125 break;
1126 default:
1127 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1128 __func__, src_type, dest_type));
1129 return NC_EBADTYPE;
1130 }
1131 break;
1132
1133 case NC_UINT64:
1134 switch (dest_type)
1135 {
1136 case NC_UBYTE:
1137 for (ulip = (unsigned long long *)src, ubp = dest; count < len; count++)
1138 {
1139 if (*ulip > X_UCHAR_MAX)
1140 (*range_error)++;
1141 *ubp++ = (unsigned char)*ulip++;
1142 }
1143 break;
1144 case NC_BYTE:
1145 for (ulip = (unsigned long long *)src, bp = dest; count < len; count++)
1146 {
1147 if (*ulip > X_SCHAR_MAX)
1148 (*range_error)++;
1149 *bp++ = (signed char)*ulip++;
1150 }
1151 break;
1152 case NC_SHORT:
1153 for (ulip = (unsigned long long *)src, sp = dest; count < len; count++)
1154 {
1155 if (*ulip > X_SHORT_MAX)
1156 (*range_error)++;
1157 *sp++ = (short)*ulip++;
1158 }
1159 break;
1160 case NC_USHORT:
1161 for (ulip = (unsigned long long *)src, usp = dest; count < len; count++)
1162 {
1163 if (*ulip > X_USHORT_MAX)
1164 (*range_error)++;
1165 *usp++ = (unsigned short)*ulip++;
1166 }
1167 break;
1168 case NC_UINT:
1169 for (ulip = (unsigned long long *)src, uip = dest; count < len; count++)
1170 {
1171 if (*ulip > X_UINT_MAX)
1172 (*range_error)++;
1173 *uip++ = (unsigned int)*ulip++;
1174 }
1175 break;
1176 case NC_INT:
1177 for (ulip = (unsigned long long *)src, ip = dest; count < len; count++)
1178 {
1179 if (*ulip > X_INT_MAX)
1180 (*range_error)++;
1181 *ip++ = (int)*ulip++;
1182 }
1183 break;
1184 case NC_INT64:
1185 for (ulip = (unsigned long long *)src, lip = dest; count < len; count++)
1186 {
1187 if (*ulip > X_INT64_MAX)
1188 (*range_error)++;
1189 *lip++ = (long long)*ulip++;
1190 }
1191 break;
1192 case NC_UINT64:
1193 for (ulip = (unsigned long long *)src, ulip1 = dest; count < len; count++)
1194 *ulip1++ = *ulip++;
1195 break;
1196 case NC_FLOAT:
1197 for (ulip = (unsigned long long *)src, fp = dest; count < len; count++)
1198 *fp++ = (float)*ulip++;
1199 break;
1200 case NC_DOUBLE:
1201 for (ulip = (unsigned long long *)src, dp = dest; count < len; count++)
1202 *dp++ = (double)*ulip++;
1203 break;
1204 default:
1205 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1206 __func__, src_type, dest_type));
1207 return NC_EBADTYPE;
1208 }
1209 break;
1210
1211 case NC_FLOAT:
1212 switch (dest_type)
1213 {
1214 case NC_UBYTE:
1215 for (fp = (float *)src, ubp = dest; count < len; count++)
1216 {
1217 if (*fp > X_UCHAR_MAX || *fp < 0)
1218 (*range_error)++;
1219 *ubp++ = (unsigned char)*fp++;
1220 }
1221 break;
1222 case NC_BYTE:
1223 for (fp = (float *)src, bp = dest; count < len; count++)
1224 {
1225 if (*fp > (double)X_SCHAR_MAX || *fp < (double)X_SCHAR_MIN)
1226 (*range_error)++;
1227 *bp++ = (signed char)*fp++;
1228 }
1229 break;
1230 case NC_SHORT:
1231 for (fp = (float *)src, sp = dest; count < len; count++)
1232 {
1233 if (*fp > (double)X_SHORT_MAX || *fp < (double)X_SHORT_MIN)
1234 (*range_error)++;
1235 *sp++ = (short)*fp++;
1236 }
1237 break;
1238 case NC_USHORT:
1239 for (fp = (float *)src, usp = dest; count < len; count++)
1240 {
1241 if (*fp > X_USHORT_MAX || *fp < 0)
1242 (*range_error)++;
1243 *usp++ = (unsigned short)*fp++;
1244 }
1245 break;
1246 case NC_UINT:
1247 for (fp = (float *)src, uip = dest; count < len; count++)
1248 {
1249 if (*fp > (float)X_UINT_MAX || *fp < 0)
1250 (*range_error)++;
1251 *uip++ = (unsigned int)*fp++;
1252 }
1253 break;
1254 case NC_INT:
1255 for (fp = (float *)src, ip = dest; count < len; count++)
1256 {
1257 if (*fp > (double)X_INT_MAX || *fp < (double)X_INT_MIN)
1258 (*range_error)++;
1259 *ip++ = (int)*fp++;
1260 }
1261 break;
1262 case NC_INT64:
1263 for (fp = (float *)src, lip = dest; count < len; count++)
1264 {
1265 if (*fp > (float)X_INT64_MAX || *fp <X_INT64_MIN)
1266 (*range_error)++;
1267 *lip++ = (long long)*fp++;
1268 }
1269 break;
1270 case NC_UINT64:
1271 for (fp = (float *)src, ulip = dest; count < len; count++)
1272 {
1273 if (*fp > (float)X_UINT64_MAX || *fp < 0)
1274 (*range_error)++;
1275 *ulip++ = (unsigned long long)*fp++;
1276 }
1277 break;
1278 case NC_FLOAT:
1279 for (fp = (float *)src, fp1 = dest; count < len; count++)
1280 *fp1++ = *fp++;
1281 break;
1282 case NC_DOUBLE:
1283 for (fp = (float *)src, dp = dest; count < len; count++)
1284 *dp++ = *fp++;
1285 break;
1286 default:
1287 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1288 __func__, src_type, dest_type));
1289 return NC_EBADTYPE;
1290 }
1291 break;
1292
1293 case NC_DOUBLE:
1294 switch (dest_type)
1295 {
1296 case NC_UBYTE:
1297 for (dp = (double *)src, ubp = dest; count < len; count++)
1298 {
1299 if (*dp > X_UCHAR_MAX || *dp < 0)
1300 (*range_error)++;
1301 *ubp++ = (unsigned char)*dp++;
1302 }
1303 break;
1304 case NC_BYTE:
1305 for (dp = (double *)src, bp = dest; count < len; count++)
1306 {
1307 if (*dp > X_SCHAR_MAX || *dp < X_SCHAR_MIN)
1308 (*range_error)++;
1309 *bp++ = (signed char)*dp++;
1310 }
1311 break;
1312 case NC_SHORT:
1313 for (dp = (double *)src, sp = dest; count < len; count++)
1314 {
1315 if (*dp > X_SHORT_MAX || *dp < X_SHORT_MIN)
1316 (*range_error)++;
1317 *sp++ = (short)*dp++;
1318 }
1319 break;
1320 case NC_USHORT:
1321 for (dp = (double *)src, usp = dest; count < len; count++)
1322 {
1323 if (*dp > X_USHORT_MAX || *dp < 0)
1324 (*range_error)++;
1325 *usp++ = (unsigned short)*dp++;
1326 }
1327 break;
1328 case NC_UINT:
1329 for (dp = (double *)src, uip = dest; count < len; count++)
1330 {
1331 if (*dp > X_UINT_MAX || *dp < 0)
1332 (*range_error)++;
1333 *uip++ = (unsigned int)*dp++;
1334 }
1335 break;
1336 case NC_INT:
1337 for (dp = (double *)src, ip = dest; count < len; count++)
1338 {
1339 if (*dp > X_INT_MAX || *dp < X_INT_MIN)
1340 (*range_error)++;
1341 *ip++ = (int)*dp++;
1342 }
1343 break;
1344 case NC_INT64:
1345 for (dp = (double *)src, lip = dest; count < len; count++)
1346 {
1347 if (*dp > (double)X_INT64_MAX || *dp < X_INT64_MIN)
1348 (*range_error)++;
1349 *lip++ = (long long)*dp++;
1350 }
1351 break;
1352 case NC_UINT64:
1353 for (dp = (double *)src, ulip = dest; count < len; count++)
1354 {
1355 if (*dp > (double)X_UINT64_MAX || *dp < 0)
1356 (*range_error)++;
1357 *ulip++ = (unsigned long long)*dp++;
1358 }
1359 break;
1360 case NC_FLOAT:
1361 for (dp = (double *)src, fp = dest; count < len; count++)
1362 {
1363 if (isgreater(*dp, X_FLOAT_MAX) || isless(*dp, X_FLOAT_MIN))
1364 (*range_error)++;
1365 *fp++ = (float)*dp++;
1366 }
1367 break;
1368 case NC_DOUBLE:
1369 for (dp = (double *)src, dp1 = dest; count < len; count++)
1370 *dp1++ = *dp++;
1371 break;
1372 default:
1373 LOG((0, "%s: unexpected dest type. src_type %d, dest_type %d",
1374 __func__, src_type, dest_type));
1375 return NC_EBADTYPE;
1376 }
1377 break;
1378
1379 default:
1380 LOG((0, "%s: unexpected src type. src_type %d, dest_type %d",
1381 __func__, src_type, dest_type));
1382 return NC_EBADTYPE;
1383 }
1384
1385 /* If quantize is in use, determine masks, copy the data, do the
1386 * quantization. */
1387 if (quantize_mode == NC_QUANTIZE_BITGROOM)
1388 {
1389 if (dest_type == NC_FLOAT)
1390 {
1391 /* BitGroom: alternately shave and set LSBs */
1392 op1.fp = (float *)dest;
1393 u32_ptr = op1.ui32p;
1394 /* Do not quantize _FillValue, +/- zero, or NaN */
1395 for (idx = 0L; idx < len; idx += 2L)
1396 if ((val_flt=op1.fp[idx]) != mss_val_cmp_flt && val_flt != 0.0f && !isnan(val_flt))
1397 u32_ptr[idx] &= msk_f32_u32_zro;
1398 for (idx = 1L; idx < len; idx += 2L)
1399 if ((val_flt=op1.fp[idx]) != mss_val_cmp_flt && val_flt != 0.0f && !isnan(val_flt))
1400 u32_ptr[idx] |= msk_f32_u32_one;
1401 }
1402 else
1403 {
1404 /* BitGroom: alternately shave and set LSBs. */
1405 op1.dp = (double *)dest;
1406 u64_ptr = op1.ui64p;
1407 /* Do not quantize _FillValue, +/- zero, or NaN */
1408 for (idx = 0L; idx < len; idx += 2L)
1409 if ((val_dbl=op1.dp[idx]) != mss_val_cmp_dbl && val_dbl != 0.0 && !isnan(val_dbl))
1410 u64_ptr[idx] &= msk_f64_u64_zro;
1411 for (idx = 1L; idx < len; idx += 2L)
1412 if ((val_dbl=op1.dp[idx]) != mss_val_cmp_dbl && val_dbl != 0.0 && !isnan(val_dbl))
1413 u64_ptr[idx] |= msk_f64_u64_one;
1414 }
1415 } /* endif BitGroom */
1416
1417 if (quantize_mode == NC_QUANTIZE_BITROUND)
1418 {
1419 if (dest_type == NC_FLOAT)
1420 {
1421 /* BitRound: Quantize to user-specified NSB with IEEE-rounding */
1422 op1.fp = (float *)dest;
1423 u32_ptr = op1.ui32p;
1424 for (idx = 0L; idx < len; idx++){
1425 /* Do not quantize _FillValue, +/- zero, or NaN */
1426 if ((val_flt=op1.fp[idx]) != mss_val_cmp_flt && val_flt != 0.0f && !isnan(val_flt)){
1427 u32_ptr[idx] += msk_f32_u32_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1428 u32_ptr[idx] &= msk_f32_u32_zro; /* Shave it */
1429 }
1430 }
1431 }
1432 else
1433 {
1434 /* BitRound: Quantize to user-specified NSB with IEEE-rounding */
1435 op1.dp = (double *)dest;
1436 u64_ptr = op1.ui64p;
1437 for (idx = 0L; idx < len; idx++){
1438 /* Do not quantize _FillValue, +/- zero, or NaN */
1439 if ((val_dbl=op1.dp[idx]) != mss_val_cmp_dbl && val_dbl != 0.0 && !isnan(val_dbl)){
1440 u64_ptr[idx] += msk_f64_u64_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1441 u64_ptr[idx] &= msk_f64_u64_zro; /* Shave it */
1442 }
1443 }
1444 }
1445 } /* endif BitRound */
1446
1447 if (quantize_mode == NC_QUANTIZE_GRANULARBR)
1448 {
1449 if (dest_type == NC_FLOAT)
1450 {
1451 /* Granular BitRound */
1452 op1.fp = (float *)dest;
1453 u32_ptr = op1.ui32p;
1454 for (idx = 0L; idx < len; idx++)
1455 {
1456 /* Do not quantize _FillValue, +/- zero, or NaN */
1457 if((val_dbl=op1.fp[idx]) != mss_val_cmp_flt && val_dbl != 0.0 && !isnan(val_dbl))
1458 {
1459 mnt = frexp(val_dbl, &xpn_bs2); /* DGG19 p. 4102 (8) */
1460 mnt_fabs = fabs(mnt);
1461 mnt_log10_fabs = log10(mnt_fabs);
1462 /* 20211003 Continuous determination of dgt_nbr improves CR by ~10% */
1463 dgt_nbr = (int)floor(xpn_bs2 * dgt_per_bit + mnt_log10_fabs) + 1; /* DGG19 p. 4102 (8.67) */
1464 qnt_pwr = (int)floor(bit_per_dgt * (dgt_nbr - nsd)); /* DGG19 p. 4101 (7) */
1465 prc_bnr_xpl_rqr = mnt_fabs == 0.0 ? 0 : (unsigned short)abs((int)floor(xpn_bs2 - bit_per_dgt*mnt_log10_fabs) - qnt_pwr); /* Protect against mnt = -0.0 */
1466 prc_bnr_xpl_rqr--; /* 20211003 Reduce formula result by 1 bit: Passes all tests, improves CR by ~10% */
1467
1468 bit_xpl_nbr_zro = BIT_XPL_NBR_SGN_FLT - prc_bnr_xpl_rqr;
1469 msk_f32_u32_zro = 0U; /* Zero all bits */
1470 msk_f32_u32_zro = ~msk_f32_u32_zro; /* Turn all bits to ones */
1471 /* Bit Shave mask for AND: Left shift zeros into bits to be rounded, leave ones in untouched bits */
1472 msk_f32_u32_zro <<= bit_xpl_nbr_zro;
1473 /* Bit Set mask for OR: Put ones into bits to be set, zeros in untouched bits */
1474 msk_f32_u32_one = ~msk_f32_u32_zro;
1475 msk_f32_u32_hshv = msk_f32_u32_one & (msk_f32_u32_zro >> 1); /* Set one bit: the MSB of LSBs */
1476 u32_ptr[idx] += msk_f32_u32_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1477 u32_ptr[idx] &= msk_f32_u32_zro; /* Shave it */
1478
1479 } /* !mss_val_cmp_flt */
1480
1481 }
1482 }
1483 else
1484 {
1485 /* Granular BitRound */
1486 op1.dp = (double *)dest;
1487 u64_ptr = op1.ui64p;
1488 for (idx = 0L; idx < len; idx++)
1489 {
1490 /* Do not quantize _FillValue, +/- zero, or NaN */
1491 if((val_dbl=op1.dp[idx]) != mss_val_cmp_dbl && val_dbl != 0.0 && !isnan(val_dbl))
1492 {
1493 mnt = frexp(val_dbl, &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 : (unsigned short)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_DBL - prc_bnr_xpl_rqr;
1503 msk_f64_u64_zro = 0ULL; /* Zero all bits */
1504 msk_f64_u64_zro = ~msk_f64_u64_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_f64_u64_zro <<= bit_xpl_nbr_zro;
1507 /* Bit Set mask for OR: Put ones into bits to be set, zeros in untouched bits */
1508 msk_f64_u64_one = ~msk_f64_u64_zro;
1509 msk_f64_u64_hshv = msk_f64_u64_one & (msk_f64_u64_zro >> 1); /* Set one bit: the MSB of LSBs */
1510 u64_ptr[idx] += msk_f64_u64_hshv; /* Add 1 to the MSB of LSBs, carry 1 to mantissa or even exponent */
1511 u64_ptr[idx] &= msk_f64_u64_zro; /* Shave it */
1512
1513 } /* !mss_val_cmp_dbl */
1514
1515 }
1516 }
1517 } /* endif GranularBR */
1518
1519 return NC_NOERR;
1520}
1521
1533int
1534nc4_get_fill_value(NC_FILE_INFO_T *h5, NC_VAR_INFO_T *var, void **fillp)
1535{
1536 size_t size;
1537 int retval;
1538
1539 /* Find out how much space we need for this type's fill value. */
1540 if (var->type_info->nc_type_class == NC_VLEN)
1541 size = sizeof(nc_vlen_t);
1542 else if (var->type_info->nc_type_class == NC_STRING)
1543 size = sizeof(char *);
1544 else
1545 {
1546 if ((retval = nc4_get_typelen_mem(h5, var->type_info->hdr.id, &size)))
1547 return retval;
1548 }
1549 assert(size);
1550
1551 /* Allocate the space. */
1552 if (!((*fillp) = calloc(1, size)))
1553 return NC_ENOMEM;
1554
1555 /* If the user has set a fill_value for this var, use, otherwise
1556 * find the default fill value. */
1557 if (var->fill_value)
1558 {
1559 LOG((4, "Found a fill value for var %s", var->hdr.name));
1560 if (var->type_info->nc_type_class == NC_VLEN)
1561 {
1562 nc_vlen_t *in_vlen = (nc_vlen_t *)(var->fill_value), *fv_vlen = (nc_vlen_t *)(*fillp);
1563 size_t basetypesize = 0;
1564
1565 if((retval=nc4_get_typelen_mem(h5, var->type_info->u.v.base_nc_typeid, &basetypesize)))
1566 return retval;
1567
1568 fv_vlen->len = in_vlen->len;
1569 if (!(fv_vlen->p = malloc(basetypesize * in_vlen->len)))
1570 {
1571 free(*fillp);
1572 *fillp = NULL;
1573 return NC_ENOMEM;
1574 }
1575 memcpy(fv_vlen->p, in_vlen->p, in_vlen->len * basetypesize);
1576 }
1577 else if (var->type_info->nc_type_class == NC_STRING)
1578 {
1579 if (*(char **)var->fill_value)
1580 if (!(**(char ***)fillp = strdup(*(char **)var->fill_value)))
1581 {
1582 free(*fillp);
1583 *fillp = NULL;
1584 return NC_ENOMEM;
1585 }
1586 }
1587 else
1588 memcpy((*fillp), var->fill_value, size);
1589 }
1590 else
1591 {
1592 if (nc4_get_default_fill_value(var->type_info, *fillp))
1593 {
1594 /* Note: release memory, but don't return error on failure */
1595 free(*fillp);
1596 *fillp = NULL;
1597 }
1598 }
1599
1600 return NC_NOERR;
1601}
1602
1615int
1616nc4_get_typelen_mem(NC_FILE_INFO_T *h5, nc_type xtype, size_t *len)
1617{
1618 NC_TYPE_INFO_T *type;
1619 int retval;
1620
1621 LOG((4, "%s xtype: %d", __func__, xtype));
1622 assert(len);
1623
1624 /* If this is an atomic type, the answer is easy. */
1625 switch (xtype)
1626 {
1627 case NC_BYTE:
1628 case NC_CHAR:
1629 case NC_UBYTE:
1630 *len = sizeof(char);
1631 return NC_NOERR;
1632 case NC_SHORT:
1633 case NC_USHORT:
1634 *len = sizeof(short);
1635 return NC_NOERR;
1636 case NC_INT:
1637 case NC_UINT:
1638 *len = sizeof(int);
1639 return NC_NOERR;
1640 case NC_FLOAT:
1641 *len = sizeof(float);
1642 return NC_NOERR;
1643 case NC_DOUBLE:
1644 *len = sizeof(double);
1645 return NC_NOERR;
1646 case NC_INT64:
1647 case NC_UINT64:
1648 *len = sizeof(long long);
1649 return NC_NOERR;
1650 case NC_STRING:
1651 *len = sizeof(char *);
1652 return NC_NOERR;
1653 }
1654
1655 /* See if var is compound type. */
1656 if ((retval = nc4_find_type(h5, xtype, &type)))
1657 return retval;
1658
1659 if (!type)
1660 return NC_EBADTYPE;
1661
1662 *len = type->size;
1663
1664 LOG((5, "type->size: %d", type->size));
1665
1666 return NC_NOERR;
1667}
1668
1669
1683int
1684nc4_check_chunksizes(NC_GRP_INFO_T *grp, NC_VAR_INFO_T *var, const size_t *chunksizes)
1685{
1686 double dprod;
1687 size_t type_len;
1688 int d;
1689 int retval;
1690
1691 if ((retval = nc4_get_typelen_mem(grp->nc4_info, var->type_info->hdr.id, &type_len)))
1692 return retval;
1693 if (var->type_info->nc_type_class == NC_VLEN)
1694 dprod = (double)sizeof(nc_vlen_t);
1695 else
1696 dprod = (double)type_len;
1697 for (d = 0; d < var->ndims; d++)
1698 dprod *= (double)chunksizes[d];
1699
1700 if (dprod > (double) NC_MAX_UINT)
1701 return NC_EBADCHUNK;
1702
1703 return NC_NOERR;
1704}
1705
1717int
1718nc4_find_default_chunksizes2(NC_GRP_INFO_T *grp, NC_VAR_INFO_T *var)
1719{
1720 int d;
1721 size_t type_size;
1722 float num_values = 1, num_unlim = 0;
1723 int retval;
1724 size_t suggested_size;
1725#ifdef LOGGING
1726 double total_chunk_size;
1727#endif
1728
1729 if (var->type_info->nc_type_class == NC_STRING)
1730 type_size = sizeof(char *);
1731 else
1732 type_size = var->type_info->size;
1733
1734#ifdef LOGGING
1735 /* Later this will become the total number of bytes in the default
1736 * chunk. */
1737 total_chunk_size = (double) type_size;
1738#endif
1739
1740 if(var->chunksizes == NULL) {
1741 if((var->chunksizes = calloc(1,sizeof(size_t)*var->ndims)) == NULL)
1742 return NC_ENOMEM;
1743 }
1744
1745 /* How many values in the variable (or one record, if there are
1746 * unlimited dimensions). */
1747 for (d = 0; d < var->ndims; d++)
1748 {
1749 assert(var->dim[d]);
1750 if (! var->dim[d]->unlimited)
1751 num_values *= (float)var->dim[d]->len;
1752 else {
1753 num_unlim++;
1754 var->chunksizes[d] = 1; /* overwritten below, if all dims are unlimited */
1755 }
1756 }
1757 /* Special case to avoid 1D vars with unlim dim taking huge amount
1758 of space (DEFAULT_CHUNK_SIZE bytes). Instead we limit to about
1759 4KB */
1760 if (var->ndims == 1 && num_unlim == 1) {
1761 if (DEFAULT_CHUNK_SIZE / type_size <= 0)
1762 suggested_size = 1;
1763 else if (DEFAULT_CHUNK_SIZE / type_size > DEFAULT_1D_UNLIM_SIZE)
1764 suggested_size = DEFAULT_1D_UNLIM_SIZE;
1765 else
1766 suggested_size = DEFAULT_CHUNK_SIZE / type_size;
1767 var->chunksizes[0] = suggested_size / type_size;
1768 LOG((4, "%s: name %s dim %d DEFAULT_CHUNK_SIZE %d num_values %f type_size %d "
1769 "chunksize %ld", __func__, var->hdr.name, d, DEFAULT_CHUNK_SIZE, num_values, type_size, var->chunksizes[0]));
1770 }
1771 if (var->ndims > 1 && (float)var->ndims == num_unlim) { /* all dims unlimited */
1772 suggested_size = (size_t)pow((double)DEFAULT_CHUNK_SIZE/(double)type_size, 1.0/(double)(var->ndims));
1773 for (d = 0; d < var->ndims; d++)
1774 {
1775 var->chunksizes[d] = suggested_size ? suggested_size : 1;
1776 LOG((4, "%s: name %s dim %d DEFAULT_CHUNK_SIZE %d num_values %f type_size %d "
1777 "chunksize %ld", __func__, var->hdr.name, d, DEFAULT_CHUNK_SIZE, num_values, type_size, var->chunksizes[d]));
1778 }
1779 }
1780
1781 /* Pick a chunk length for each dimension, if one has not already
1782 * been picked above. */
1783 for (d = 0; d < var->ndims; d++)
1784 if (!var->chunksizes[d])
1785 {
1786 suggested_size = (size_t)(pow((double)DEFAULT_CHUNK_SIZE/(num_values * (double)type_size),
1787 1.0/(double)((double)var->ndims - num_unlim)) * (double)var->dim[d]->len - .5);
1788 if (suggested_size > var->dim[d]->len)
1789 suggested_size = var->dim[d]->len;
1790 var->chunksizes[d] = suggested_size ? suggested_size : 1;
1791 LOG((4, "%s: name %s dim %d DEFAULT_CHUNK_SIZE %d num_values %f type_size %d "
1792 "chunksize %ld", __func__, var->hdr.name, d, DEFAULT_CHUNK_SIZE, num_values, type_size, var->chunksizes[d]));
1793 }
1794
1795#ifdef LOGGING
1796 /* Find total chunk size. */
1797 for (d = 0; d < var->ndims; d++)
1798 total_chunk_size *= (double) var->chunksizes[d];
1799 LOG((4, "total_chunk_size %f", total_chunk_size));
1800#endif
1801
1802 /* But did this result in a chunk that is too big? */
1803 retval = nc4_check_chunksizes(grp, var, var->chunksizes);
1804 if (retval)
1805 {
1806 /* Other error? */
1807 if (retval != NC_EBADCHUNK)
1808 return retval;
1809
1810 /* Chunk is too big! Reduce each dimension by half and try again. */
1811 for ( ; retval == NC_EBADCHUNK; retval = nc4_check_chunksizes(grp, var, var->chunksizes))
1812 for (d = 0; d < var->ndims; d++)
1813 var->chunksizes[d] = var->chunksizes[d]/2 ? var->chunksizes[d]/2 : 1;
1814 }
1815
1816 /* Do we have any big data overhangs? They can be dangerous to
1817 * babies, the elderly, or confused campers who have had too much
1818 * beer. */
1819 for (d = 0; d < var->ndims; d++)
1820 {
1821 size_t num_chunks;
1822 size_t overhang;
1823 assert(var->chunksizes[d] > 0);
1824 num_chunks = (var->dim[d]->len + var->chunksizes[d] - 1) / var->chunksizes[d];
1825 if(num_chunks > 0) {
1826 overhang = (num_chunks * var->chunksizes[d]) - var->dim[d]->len;
1827 var->chunksizes[d] -= overhang / num_chunks;
1828 }
1829 }
1830
1831
1832 return NC_NOERR;
1833}
1834
1846int
1847nc4_get_default_fill_value(NC_TYPE_INFO_T* tinfo, void *fill_value)
1848{
1849 if(tinfo->hdr.id > NC_NAT && tinfo->hdr.id <= NC_MAX_ATOMIC_TYPE)
1850 return nc4_get_default_atomic_fill_value(tinfo->hdr.id,fill_value);
1851#ifdef USE_NETCDF4
1852 switch(tinfo->nc_type_class) {
1853 case NC_ENUM:
1854 return nc4_get_default_atomic_fill_value(tinfo->u.e.base_nc_typeid,fill_value);
1855 case NC_OPAQUE:
1856 case NC_VLEN:
1857 case NC_COMPOUND:
1858 if(fill_value)
1859 memset(fill_value,0,tinfo->size);
1860 break;
1861 default: return NC_EBADTYPE;
1862 }
1863#endif
1864 return NC_NOERR;
1865}
1866
1878int
1879nc4_get_default_atomic_fill_value(nc_type xtype, void *fill_value)
1880{
1881 switch (xtype)
1882 {
1883 case NC_CHAR:
1884 *(char *)fill_value = NC_FILL_CHAR;
1885 break;
1886
1887 case NC_STRING:
1888 *(char **)fill_value = strdup(NC_FILL_STRING);
1889 break;
1890
1891 case NC_BYTE:
1892 *(signed char *)fill_value = NC_FILL_BYTE;
1893 break;
1894
1895 case NC_SHORT:
1896 *(short *)fill_value = NC_FILL_SHORT;
1897 break;
1898
1899 case NC_INT:
1900 *(int *)fill_value = NC_FILL_INT;
1901 break;
1902
1903 case NC_UBYTE:
1904 *(unsigned char *)fill_value = NC_FILL_UBYTE;
1905 break;
1906
1907 case NC_USHORT:
1908 *(unsigned short *)fill_value = NC_FILL_USHORT;
1909 break;
1910
1911 case NC_UINT:
1912 *(unsigned int *)fill_value = NC_FILL_UINT;
1913 break;
1914
1915 case NC_INT64:
1916 *(long long *)fill_value = NC_FILL_INT64;
1917 break;
1918
1919 case NC_UINT64:
1920 *(unsigned long long *)fill_value = NC_FILL_UINT64;
1921 break;
1922
1923 case NC_FLOAT:
1924 *(float *)fill_value = NC_FILL_FLOAT;
1925 break;
1926
1927 case NC_DOUBLE:
1928 *(double *)fill_value = NC_FILL_DOUBLE;
1929 break;
1930
1931 default:
1932 return NC_EINVAL;
1933 }
1934 return NC_NOERR;
1935}
1936
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:95
#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:459
#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:852
#define NC_EFILTER
Filter operation failed.
Definition netcdf.h:562
#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:387
size_t len
Length of VL data (in base type units)
Definition netcdf.h:851
#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:497
#define NC_MAX_INT
Max or min values for a type.
Definition netcdf.h:94
#define NC_COMPOUND
compound types
Definition netcdf.h:56
#define NC_CHUNKED
In HDF5 files you can set storage for each variable to be either contiguous or chunked,...
Definition netcdf.h:353
#define NC_SHORT
signed 2 byte integer
Definition netcdf.h:37
#define NC_QUANTIZE_GRANULARBR
Use Granular BitRound quantization.
Definition netcdf.h:386
#define NC_FILL_UINT64
Default fill value.
Definition netcdf.h:77
#define NC_QUANTIZE_BITGROOM
Use BitGroom quantization.
Definition netcdf.h:385
#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:303
#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:471
#define NC_EINVAL
Invalid Argument.
Definition netcdf.h:427
#define NC_MAX_NAME
Maximum for classic library.
Definition netcdf.h:330
#define NC_NOERR
No Error.
Definition netcdf.h:417
#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:384
#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:543
#define NC_STRING
string
Definition netcdf.h:47
#define NC_EBADCHUNK
Bad chunksize.
Definition netcdf.h:556
#define NC_FILL_FLOAT
Default fill value.
Definition netcdf.h:71
#define NC_ENOFILTER
Filter not defined on variable.
Definition netcdf.h:566
#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:496
#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:850
#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