1 /*********************************************************************
2 * Copyright 2018, University Corporation for Atmospheric Research
3 * See netcdf/README file for copying and redistribution conditions.
4 * "$Id: nciter.c 400 2010-08-27 21:02:52Z russ $"
5 *********************************************************************/
6
7 #include "config.h" /* for USE_NETCDF4 macro */
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include <netcdf.h>
12 #include "utils.h"
13 #include "nciter.h"
14
15 /* Initialize block iteration for variables, including those that
16 * won't fit in the copy buffer all at once. Returns error if
17 * variable is chunked but size of chunks is too big to fit in bufsize
18 * bytes. */
19 static int
nc_blkio_init(size_t bufsize,size_t value_size,int rank,int chunked,nciter_t * iter)20 nc_blkio_init(size_t bufsize, /* size in bytes of in-memory copy buffer */
21 size_t value_size, /* size in bytes of each variable element */
22 int rank, /* number of dimensions for variable */
23 int chunked, /* 1 if variable is chunked, 0 otherwise */
24 nciter_t *iter /* returned iteration state, don't mess with it */
25 ) {
26 int stat = NC_NOERR;
27 int i;
28 long long prod;
29 size_t *dims = iter->dimsizes;
30
31 iter->rank = rank;
32 iter->first = 1;
33 iter->more = 1;
34 iter->chunked = chunked;
35 prod = value_size;
36 if(iter->chunked == 0) { /* contiguous */
37 iter->right_dim = rank - 1;
38 for(i = rank; i > 0; i--) {
39 if(prod*dims[i-1] <= bufsize) {
40 prod *= dims[i-1];
41 iter->right_dim--;
42 } else {
43 break;
44 }
45 }
46 if (i > 0) { /* variable won't fit in bufsize bytes */
47 iter->rows = bufsize/prod;
48 iter->numrows = dims[iter->right_dim] / iter->rows;
49 iter->leftover = dims[iter->right_dim] - iter->numrows * iter->rows;
50 iter->cur = 1;
51 iter->inc = iter->rows;
52 return stat;
53 }
54 /* else, variable will fit in bufsize bytes of memory. */
55 iter->right_dim = 0;
56 iter->rows = dims[0];
57 iter->inc = 0;
58 return stat;
59 }
60 /* else, handle chunked case */
61 for(i = 0; i < rank; i++) {
62 prod *= iter->chunksizes[i];
63 }
64 if(prod > bufsize) {
65 stat = NC_ENOMEM;
66 fprintf(stderr, "chunksize (= %ld) > copy_buffer size (= %ld)\n",
67 (long)prod, (long)bufsize);
68 }
69 return stat;
70 }
71
72 /* From netCDF type in group igrp, get size in memory needed for each
73 * value. Wouldn't be needed if nc_inq_type() was a netCDF-3 function
74 * too. */
75 static int
inq_value_size(int igrp,nc_type vartype,size_t * value_sizep)76 inq_value_size(int igrp, nc_type vartype, size_t *value_sizep) {
77 int stat = NC_NOERR;
78
79 #ifdef USE_NETCDF4
80 NC_CHECK(nc_inq_type(igrp, vartype, NULL, value_sizep));
81 #else
82 switch(vartype) {
83 case NC_BYTE:
84 *value_sizep = sizeof(signed char);
85 break;
86 case NC_CHAR:
87 *value_sizep = sizeof(char);
88 break;
89 case NC_SHORT:
90 *value_sizep = sizeof(short);
91 break;
92 case NC_INT:
93 *value_sizep = sizeof(int);
94 break;
95 case NC_FLOAT:
96 *value_sizep = sizeof(float);
97 break;
98 case NC_DOUBLE:
99 *value_sizep = sizeof(double);
100 break;
101 default:
102 NC_CHECK(NC_EBADTYPE);
103 break;
104 }
105 #endif /* USE_NETCDF4 */
106 return stat;
107 }
108
109 /*
110 * Updates a vector of size_t, odometer style. Returns 0 if odometer
111 * overflowed, else 1.
112 */
113 static int
up_start(int ndims,const size_t * dims,int incdim,size_t inc,size_t * odom)114 up_start(
115 int ndims, /* Number of dimensions */
116 const size_t *dims, /* The "odometer" limits for each dimension */
117 int incdim, /* the odmometer increment dimension */
118 size_t inc, /* the odometer increment for that dimension */
119 size_t* odom /* The "odometer" vector to be updated */
120 )
121 {
122 int id;
123 int ret = 1;
124
125 if(inc == 0) {
126 return 0;
127 }
128 odom[incdim] += inc;
129 for (id = incdim; id > 0; id--) {
130 if(odom[id] >= dims[id]) {
131 odom[id-1]++;
132 odom[id] -= dims[id];
133 }
134 }
135 if (odom[0] >= dims[0])
136 ret = 0;
137 return ret;
138 }
139
140 /*
141 * Updates a vector of size_t, odometer style, for chunk access.
142 * Returns 0 if odometer overflowed, else 1.
143 */
144 static int
up_start_by_chunks(int ndims,const size_t * dims,const size_t * chunks,size_t * odom)145 up_start_by_chunks(
146 int ndims, /* Number of dimensions */
147 const size_t *dims, /* The "odometer" limits for each dimension */
148 const size_t *chunks, /* the odometer increments for each dimension */
149 size_t* odom /* The "odometer" vector to be updated */
150 )
151 {
152 int incdim = ndims - 1;
153 int id;
154 int ret = 1;
155
156 odom[incdim] += chunks[incdim];
157 for (id = incdim; id > 0; id--) {
158 if(odom[id] >= dims[id]) {
159 odom[id-1] += chunks[id-1];
160 /* odom[id] -= dims[id]; */
161 odom[id] = 0;
162 }
163 }
164 if (odom[0] >= dims[0])
165 ret = 0;
166 return ret;
167 }
168
169 /* Begin public interfaces */
170
171 /* Initialize iteration for a variable. Just a wrapper for
172 * nc_blkio_init() that makes the netCDF calls needed to initialize
173 * lower-level iterator. */
174 int
nc_get_iter(int ncid,int varid,size_t bufsize,nciter_t ** iterpp)175 nc_get_iter(int ncid,
176 int varid,
177 size_t bufsize, /* size in bytes of memory buffer */
178 nciter_t **iterpp /* returned opaque iteration state */)
179 {
180 int stat = NC_NOERR;
181 nciter_t *iterp;
182 nc_type vartype;
183 size_t value_size = 0; /* size in bytes of each variable element */
184 int ndims; /* number of dimensions for variable */
185 int *dimids;
186 long long nvalues = 1;
187 int dim;
188 int chunked = 0;
189
190 /* Caller should free this by calling nc_free_iter(iterp) */
191 iterp = (nciter_t *) emalloc(sizeof(nciter_t));
192 memset((void*)iterp,0,sizeof(nciter_t)); /* make sure it is initialized */
193
194 NC_CHECK(nc_inq_varndims(ncid, varid, &ndims));
195
196 dimids = (int *) emalloc((ndims + 1) * sizeof(int));
197
198 iterp->dimsizes = (size_t *) emalloc((ndims + 1) * sizeof(size_t));
199 iterp->chunksizes = (size_t *) emalloc((ndims + 1) * sizeof(size_t));
200
201 NC_CHECK(nc_inq_vardimid (ncid, varid, dimids));
202 for(dim = 0; dim < ndims; dim++) {
203 size_t len;
204 NC_CHECK(nc_inq_dimlen(ncid, dimids[dim], &len));
205 nvalues *= len;
206 iterp->dimsizes[dim] = len;
207 }
208 NC_CHECK(nc_inq_vartype(ncid, varid, &vartype));
209 NC_CHECK(inq_value_size(ncid, vartype, &value_size));
210 #ifdef USE_NETCDF4
211 {
212 int contig = 1;
213 if(ndims > 0) {
214 NC_CHECK(nc_inq_var_chunking(ncid, varid, &contig, NULL));
215 }
216 if(contig == 0) { /* chunked */
217 NC_CHECK(nc_inq_var_chunking(ncid, varid, &contig, iterp->chunksizes));
218 chunked = 1;
219 }
220 }
221 #endif /* USE_NETCDF4 */
222 NC_CHECK(nc_blkio_init(bufsize, value_size, ndims, chunked, iterp));
223 iterp->to_get = 0;
224 free(dimids);
225 *iterpp = iterp;
226 return stat;
227 }
228
229 /* Iterate on blocks for variables, by updating start and count vector
230 * for next vara call. Assumes nc_get_iter called first. Returns
231 * number of variable values to get, 0 if done, negative number if
232 * error, so use like this:
233 size_t to_get;
234 while((to_get = nc_next_iter(&iter, start, count)) > 0) {
235 ... iteration ...
236 }
237 if(to_get < 0) { ... handle error ... }
238 */
239 size_t
nc_next_iter(nciter_t * iter,size_t * start,size_t * count)240 nc_next_iter(nciter_t *iter, /* returned opaque iteration state */
241 size_t *start, /* returned start vector for next vara call */
242 size_t *count /* returned count vector for next vara call */
243 ) {
244 int i;
245 /* Note: special case for chunked variables is just an
246 * optimization, the contiguous code below is OK even
247 * for chunked variables, but in general will do more I/O ... */
248 if(iter->first) {
249 if(!iter->chunked) { /* contiguous storage */
250 for(i = 0; i < iter->right_dim; i++) {
251 start[i] = 0;
252 count[i] = 1;
253 }
254 start[iter->right_dim] = 0;
255 count[iter->right_dim] = iter->rows;
256 for(i = iter->right_dim + 1; i < iter->rank; i++) {
257 start[i] = 0;
258 count[i] = iter->dimsizes[i];
259 }
260 } else { /* chunked storage */
261 for(i = 0; i < iter->rank; i++) {
262 start[i] = 0;
263 if(iter->chunksizes[i] <= iter->dimsizes[i])
264 count[i] = iter->chunksizes[i];
265 else /* can happen for variables with only unlimited dimensions */
266 count[i] = iter->dimsizes[i];
267 }
268 }
269 iter->first = 0;
270 } else {
271 if(!iter->chunked) { /* contiguous storage */
272 iter->more = up_start(iter->rank, iter->dimsizes, iter->right_dim,
273 iter->inc, start);
274 /* iterate on pieces of variable */
275 if(iter->cur < iter->numrows) {
276 iter->inc = iter->rows;
277 count[iter->right_dim] = iter->rows;
278 iter->cur++;
279 } else {
280 if(iter->leftover > 0) {
281 count[iter->right_dim] = iter->leftover;
282 iter->inc = iter->leftover;
283 iter->cur = 0;
284 }
285 }
286 } else { /* chunked storage */
287 iter->more = up_start_by_chunks(iter->rank, iter->dimsizes,
288 iter->chunksizes, start);
289 /* adjust count to stay in range of dimsizes */
290 for(i = 0; i < iter->rank; i++) {
291 int leftover = iter->dimsizes[i] - start[i];
292 if(iter->chunksizes[i] <= iter->dimsizes[i])
293 count[i] = iter->chunksizes[i];
294 else /* can happen for variables with only unlimited dimensions */
295 count[i] = iter->dimsizes[i];
296 if(leftover < count[i])
297 count[i] = leftover;
298 }
299 }
300 }
301 iter->to_get = 1;
302 for(i = 0; i < iter->rank; i++) {
303 iter->to_get *= count[i];
304 }
305 return iter->more == 0 ? 0 : iter->to_get ;
306 }
307
308 /* Free iterator and its internally allocated memory */
309 int
nc_free_iter(nciter_t * iterp)310 nc_free_iter(nciter_t *iterp) {
311 if(iterp->dimsizes)
312 free(iterp->dimsizes);
313 if(iterp->chunksizes)
314 free(iterp->chunksizes);
315 free(iterp);
316 return NC_NOERR;
317 }
318