1 /* Copyright (c) 1999-2008 Massachusetts Institute of Technology
2 *
3 * Permission is hereby granted, free of charge, to any person obtaining
4 * a copy of this software and associated documentation files (the
5 * "Software"), to deal in the Software without restriction, including
6 * without limitation the rights to use, copy, modify, merge, publish,
7 * distribute, sublicense, and/or sell copies of the Software, and to
8 * permit persons to whom the Software is furnished to do so, subject to
9 * the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be
12 * included in all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
17 * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
18 * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
19 * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
20 * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
21 */
22
23 #include <stdlib.h>
24 #include <string.h>
25
26 /* don't use new HDF5 1.8 API (which isn't even fully documented yet, grrr) */
27 #define H5_USE_16_API 1
28
29 #include <hdf5.h>
30
31 #include "arrayh5.h"
32
33 #define CHECK(cond, msg) { if (!(cond)) { fprintf(stderr, "arrayh5 error: %s\n", msg); exit(EXIT_FAILURE); } }
34
35 #define CHK_MALLOC(p, t, n) CHECK(p = (t *) malloc(sizeof(t) * (n)), "out of memory")
36
37 /* Normally, HDF5 prints out all sorts of error messages, e.g. if a dataset
38 can't be found, in addition to returning an error code. The following
39 macro can be wrapped around code to temporarily suppress error messages. */
40
41 #define SUPPRESS_HDF5_ERRORS(statements) { \
42 H5E_auto_t xxxxx_err_func; \
43 void *xxxxx_err_func_data; \
44 H5Eget_auto(&xxxxx_err_func, &xxxxx_err_func_data); \
45 H5Eset_auto(NULL, NULL); \
46 { statements; } \
47 H5Eset_auto(xxxxx_err_func, xxxxx_err_func_data); \
48 }
49
arrayh5_create_withdata(int rank,const int * dims,double * data)50 arrayh5 arrayh5_create_withdata(int rank, const int *dims, double *data)
51 {
52 arrayh5 a;
53 int i;
54
55 CHECK(rank >= 0, "non-positive rank");
56 a.rank = rank;
57
58 CHK_MALLOC(a.dims, int, rank);
59
60 a.N = 1;
61 for (i = 0; i < rank; ++i) {
62 a.dims[i] = dims[i];
63 a.N *= dims[i];
64 }
65
66 if (data)
67 a.data = data;
68 else {
69 CHK_MALLOC(a.data, double, a.N);
70 }
71 return a;
72 }
73
arrayh5_create(int rank,const int * dims)74 arrayh5 arrayh5_create(int rank, const int *dims)
75 {
76 return arrayh5_create_withdata(rank, dims, NULL);
77 }
78
arrayh5_clone(arrayh5 a)79 arrayh5 arrayh5_clone(arrayh5 a)
80 {
81 arrayh5 b = arrayh5_create(a.rank, a.dims);
82 if (a.data) memcpy(b.data, a.data, sizeof(double) * a.N);
83 return b;
84 }
85
arrayh5_destroy(arrayh5 a)86 void arrayh5_destroy(arrayh5 a)
87 {
88 free(a.dims);
89 free(a.data);
90 }
91
arrayh5_conformant(arrayh5 a,arrayh5 b)92 int arrayh5_conformant(arrayh5 a, arrayh5 b)
93 {
94 int i;
95
96 if (a.rank != b.rank)
97 return 0;
98 for (i = 0; i < a.rank; ++i)
99 if (a.dims[i] != b.dims[i])
100 return 0;
101 return 1;
102 }
103
rtranspose(int curdim,int rank,const int * dims,int curindex,int curindex_t,const double * data,double * data_t)104 static void rtranspose(int curdim, int rank, const int *dims,
105 int curindex, int curindex_t,
106 const double *data, double *data_t)
107 {
108 int prod_before = 1, prod_after = 1;
109 int i;
110
111 if (rank == 0) {
112 *data_t = *data;
113 return;
114 }
115
116 for (i = 0; i < curdim; ++i)
117 prod_before *= dims[i];
118 for (i = curdim + 1; i < rank; ++i)
119 prod_after *= dims[i];
120
121 if (curdim == rank - 1) {
122 for (i = 0; i < dims[curdim]; ++i)
123 data_t[curindex_t + i * prod_before] = data[curindex + i];
124 }
125 else {
126 for (i = 0; i < dims[curdim]; ++i)
127 rtranspose(curdim + 1, rank, dims,
128 curindex + i * prod_after,
129 curindex_t + i * prod_before,
130 data, data_t);
131 }
132 }
133
arrayh5_transpose(arrayh5 * a)134 void arrayh5_transpose(arrayh5 *a)
135 {
136 double *data_t;
137 int i;
138
139 CHK_MALLOC(data_t, double, a->N);
140 rtranspose(0, a->rank, a->dims, 0, 0, a->data, data_t);
141 free(a->data);
142 a->data = data_t;
143
144 for (i = 0; i < a->rank - 1 - i; ++i) {
145 int dummy = a->dims[i];
146 a->dims[i] = a->dims[a->rank - 1 - i];
147 a->dims[a->rank - 1 - i] = dummy;
148 }
149 }
150
arrayh5_getrange(arrayh5 a,double * min,double * max)151 void arrayh5_getrange(arrayh5 a, double *min, double *max)
152 {
153 int i;
154
155 CHECK(a.N > 0, "no elements in array");
156 *min = *max = a.data[0];
157 for (i = 1; i < a.N; ++i) {
158 if (a.data[i] < *min)
159 *min = a.data[i];
160 if (a.data[i] > *max)
161 *max = a.data[i];
162 }
163 }
164
find_dataset(hid_t group_id,const char * name,void * d)165 static herr_t find_dataset(hid_t group_id, const char *name, void *d)
166 {
167 char **dname = (char **) d;
168 H5G_stat_t info;
169
170 H5Gget_objinfo(group_id, name, 1, &info);
171 if (info.type == H5G_DATASET) {
172 CHK_MALLOC(*dname, char, strlen(name) + 1);
173 strcpy(*dname, name);
174 return 1;
175 }
176 return 0;
177 }
178
179 typedef enum { NO_ERROR = 0, OPEN_FAILED, NO_DATA, READ_FAILED, SLICE_FAILED,
180 INVALID_SLICE, INVALID_RANK, OPEN_DATA_FAILED } arrayh5_err;
181
182 const char arrayh5_read_strerror[][100] = {
183 "no error",
184 "error opening HD5 file",
185 "couldn't find data set in HDF5 file",
186 "error reading data from HDF5",
187 "error reading data slice from HDF5",
188 "invalid slice of HDF5 data",
189 "non-positive rank in HDF file",
190 "error opening data set in HDF file",
191 };
192
arrayh5_read(arrayh5 * a,const char * fname,const char * datapath,char ** dataname,int nslicedims_,const int * slicedim_,const int * islice_,const int * center_slice)193 int arrayh5_read(arrayh5 *a, const char *fname, const char *datapath,
194 char **dataname,
195 int nslicedims_, const int *slicedim_, const int *islice_,
196 const int *center_slice)
197 {
198 hid_t file_id = -1, data_id = -1, space_id = -1;
199 char *dname = NULL;
200 int err = NO_ERROR;
201 hsize_t i, rank, *dims_copy, *maxdims, *slicedim = 0;
202 int *islice = 0;
203 int *dims = 0;
204 hsize_t nslicedims = (hsize_t) nslicedims_;
205
206 CHECK(a, "NULL array passed to arrayh5_read");
207 a->dims = NULL;
208 a->data = NULL;
209
210 file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
211 if (file_id < 0) {
212 err = OPEN_FAILED;
213 goto done;
214 }
215
216 if (datapath && datapath[0]) {
217 CHK_MALLOC(dname, char, strlen(datapath) + 1);
218 strcpy(dname, datapath);
219 }
220 else {
221 if (H5Giterate(file_id, "/", NULL, find_dataset, &dname) <= 0) {
222 err = NO_DATA;
223 goto done;
224 }
225 }
226
227 data_id = H5Dopen(file_id, dname);
228 if (data_id < 0) {
229 err = OPEN_DATA_FAILED;
230 goto done;
231 }
232
233 space_id = H5Dget_space(data_id);
234 rank = H5Sget_simple_extent_ndims(space_id);
235 if (rank <= 0) {
236 err = INVALID_RANK;
237 goto done;
238 }
239
240 CHK_MALLOC(dims, int, rank);
241 CHK_MALLOC(dims_copy, hsize_t, rank);
242 CHK_MALLOC(maxdims, hsize_t, rank);
243
244 H5Sget_simple_extent_dims(space_id, dims_copy, maxdims);
245 for (i = 0; i < rank; ++i)
246 dims[i] = dims_copy[i];
247
248 free(maxdims);
249 free(dims_copy);
250
251 for (i = 0; i < nslicedims && slicedim_[i] == NO_SLICE_DIM; ++i)
252 ;
253
254 if (i == nslicedims) { /* no slices */
255 *a = arrayh5_create(rank, dims);
256
257 if (H5Dread(data_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
258 H5P_DEFAULT, (void *) a->data) < 0) {
259 err = READ_FAILED;
260 goto done;
261 }
262 }
263 else if (nslicedims > 0) {
264 int j, rank2 = rank;
265 hsize_t *start;
266 hsize_t *count;
267 hid_t mem_space_id;
268 herr_t readerr;
269
270 CHK_MALLOC(slicedim, hsize_t, nslicedims);
271 CHK_MALLOC(islice, int, nslicedims);
272
273 for (i = j = 0; i < nslicedims; ++i)
274 if (slicedim_[i] != NO_SLICE_DIM) {
275 if (slicedim_[i] == LAST_SLICE_DIM)
276 slicedim[j] = rank - 1;
277 else
278 slicedim[j] = (hsize_t) slicedim_[i];
279 if (slicedim[j] >= rank) {
280 err = INVALID_SLICE;
281 goto done;
282 }
283 islice[j] = islice_[i];
284 if (center_slice[i])
285 islice[j] += dims[slicedim[j]] / 2;
286 if (islice[j] < 0 || islice[j] >= dims[slicedim[j]]) {
287 err = INVALID_SLICE;
288 goto done;
289 }
290 j++;
291 }
292 nslicedims = j;
293
294 CHK_MALLOC(start, hsize_t, rank);
295 CHK_MALLOC(count, hsize_t, rank);
296
297 for (i = 0; i < rank; ++i) {
298 count[i] = dims[i];
299 start[i] = 0;
300 }
301 for (i = 0; i < nslicedims; ++i) {
302 start[slicedim[i]] = islice[i];
303 count[slicedim[i]] = 1;
304 }
305
306 H5Sselect_hyperslab(space_id, H5S_SELECT_SET,
307 start, NULL, count, NULL);
308
309 for (i = j = 0; i < rank; ++i)
310 if (count[i] > 1)
311 dims[j++] = count[i];
312 rank2 = j;
313
314 *a = arrayh5_create(rank2, dims);
315
316 mem_space_id = H5Screate_simple(rank, count, NULL);
317 H5Sselect_all(mem_space_id);
318
319 readerr = H5Dread(data_id, H5T_NATIVE_DOUBLE,
320 mem_space_id, space_id,
321 H5P_DEFAULT, (void *) a->data);
322
323 H5Sclose(mem_space_id);
324 free(count);
325 free(start);
326
327 if (readerr < 0) {
328 err = SLICE_FAILED;
329 goto done;
330 }
331 }
332 else {
333 err = INVALID_SLICE;
334 goto done;
335 }
336
337 done:
338 if (err != NO_ERROR)
339 arrayh5_destroy(*a);
340 free(islice);
341 free(slicedim);
342 free(dims);
343 if (space_id >= 0)
344 H5Sclose(space_id);
345 if (data_id >= 0)
346 H5Dclose(data_id);
347 if (dataname)
348 *dataname = dname;
349 else
350 free(dname);
351 if (file_id >= 0)
352 H5Fclose(file_id);
353
354 return err;
355 }
356
dataset_exists(hid_t id,const char * name)357 static int dataset_exists(hid_t id, const char *name)
358 {
359 hid_t data_id;
360 SUPPRESS_HDF5_ERRORS(data_id = H5Dopen(id, name));
361 if (data_id >= 0)
362 H5Dclose(data_id);
363 return (data_id >= 0);
364 }
365
arrayh5_write(arrayh5 a,char * filename,char * dataname,short append_data)366 void arrayh5_write(arrayh5 a, char *filename, char *dataname,
367 short append_data)
368 {
369 int i;
370 hid_t file_id, space_id, type_id, data_id;
371 hsize_t *dims_copy;
372
373 if (append_data)
374 file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
375 else
376 file_id = H5Fcreate(filename, H5F_ACC_TRUNC,
377 H5P_DEFAULT, H5P_DEFAULT);
378 CHECK(file_id >= 0, "error opening HDF5 output file");
379
380 if (dataset_exists(file_id, dataname))
381 H5Gunlink(file_id, dataname); /* delete it */
382
383 CHECK(a.rank > 0, "non-positive rank");
384 CHK_MALLOC(dims_copy, hsize_t, a.rank);
385 for (i = 0; i < a.rank; ++i)
386 dims_copy[i] = a.dims[i];
387 space_id = H5Screate_simple(a.rank, dims_copy, NULL);
388 free(dims_copy);
389
390 type_id = H5T_NATIVE_DOUBLE;
391 data_id = H5Dcreate(file_id, dataname, type_id, space_id, H5P_DEFAULT);
392 H5Sclose(space_id);
393
394 H5Dwrite(data_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, a.data);
395
396 H5Dclose(data_id);
397 H5Fclose(file_id);
398 }
399
arrayh5_read_rank(const char * fname,const char * datapath,int * rank)400 int arrayh5_read_rank(const char *fname, const char *datapath, int *rank)
401 {
402 hid_t file_id = -1, data_id = -1, space_id = -1;
403 char *dname = NULL;
404 int err = NO_ERROR;
405
406 file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
407 if (file_id < 0) {
408 err = OPEN_FAILED;
409 goto done;
410 }
411
412 if (datapath && datapath[0]) {
413 CHK_MALLOC(dname, char, strlen(datapath) + 1);
414 strcpy(dname, datapath);
415 }
416 else {
417 if (H5Giterate(file_id, "/", NULL, find_dataset, &dname) <= 0) {
418 err = NO_DATA;
419 goto done;
420 }
421 }
422
423 data_id = H5Dopen(file_id, dname);
424 if (data_id < 0) {
425 err = OPEN_DATA_FAILED;
426 goto done;
427 }
428
429 space_id = H5Dget_space(data_id);
430 *rank = H5Sget_simple_extent_ndims(space_id);
431
432 done:
433 if (space_id >= 0)
434 H5Sclose(space_id);
435 if (data_id >= 0)
436 H5Dclose(data_id);
437 free(dname);
438 if (file_id >= 0)
439 H5Fclose(file_id);
440
441 return err;
442 }
443