Actual source code: hdf5io.c
1: #include <petsc/private/viewerhdf5impl.h>
2: #include <petsclayouthdf5.h>
3: #include <petscis.h>
5: #if defined(PETSC_HAVE_HDF5)
7: struct _n_HDF5ReadCtx {
8: hid_t file, group, dataset, dataspace;
9: int lenInd, bsInd, complexInd, rdim;
10: hsize_t *dims;
11: PetscBool complexVal, dim2;
12: };
13: typedef struct _n_HDF5ReadCtx *HDF5ReadCtx;
15: PetscErrorCode PetscViewerHDF5CheckTimestepping_Internal(PetscViewer viewer, const char name[])
16: {
17: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
18: PetscBool timestepping = PETSC_FALSE;
20: PetscViewerHDF5ReadAttribute(viewer, name, "timestepping", PETSC_BOOL, ×tepping, ×tepping);
21: if (timestepping != hdf5->timestepping) {
22: char *group;
24: PetscViewerHDF5GetGroup(viewer, NULL, &group);
25: SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Dataset %s/%s stored with timesteps? %s Timestepping pushed? %s", group, name, PetscBools[timestepping], PetscBools[hdf5->timestepping]);
26: }
27: return 0;
28: }
30: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
31: {
32: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
33: HDF5ReadCtx h = NULL;
35: PetscViewerHDF5CheckTimestepping_Internal(viewer, name);
36: PetscNew(&h);
37: PetscViewerHDF5OpenGroup(viewer, NULL, &h->file, &h->group);
38: PetscCallHDF5Return(h->dataset, H5Dopen2, (h->group, name, H5P_DEFAULT));
39: PetscCallHDF5Return(h->dataspace, H5Dget_space, (h->dataset));
40: PetscViewerHDF5ReadAttribute(viewer, name, "complex", PETSC_BOOL, &h->complexVal, &h->complexVal);
41: if (!hdf5->horizontal) {
42: /* MATLAB stores column vectors horizontally */
43: PetscViewerHDF5HasAttribute(viewer, name, "MATLAB_class", &hdf5->horizontal);
44: }
45: *ctx = h;
46: return 0;
47: }
49: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
50: {
51: HDF5ReadCtx h;
53: h = *ctx;
54: PetscCallHDF5(H5Gclose, (h->group));
55: PetscCallHDF5(H5Sclose, (h->dataspace));
56: PetscCallHDF5(H5Dclose, (h->dataset));
57: PetscFree((*ctx)->dims);
58: PetscFree(*ctx);
59: return 0;
60: }
62: static PetscErrorCode PetscViewerHDF5ReadSizes_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscBool setup, PetscLayout *map_)
63: {
64: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
65: PetscInt bs, len, N;
66: PetscLayout map;
68: if (!(*map_)) PetscLayoutCreate(PetscObjectComm((PetscObject)viewer), map_);
69: map = *map_;
71: /* Get actual number of dimensions in dataset */
72: PetscCallHDF5Return(ctx->rdim, H5Sget_simple_extent_dims, (ctx->dataspace, NULL, NULL));
73: PetscMalloc1(ctx->rdim, &ctx->dims);
74: PetscCallHDF5Return(ctx->rdim, H5Sget_simple_extent_dims, (ctx->dataspace, ctx->dims, NULL));
76: /*
77: Dimensions are in this order:
78: [0] timesteps (optional)
79: [lenInd] entries (numbers or blocks)
80: ...
81: [bsInd] entries of blocks (optional)
82: [bsInd+1] real & imaginary part (optional)
83: = rdim-1
84: */
86: /* Get entries dimension index */
87: ctx->lenInd = 0;
88: if (hdf5->timestepping) ++ctx->lenInd;
90: /* Get block dimension index */
91: if (ctx->complexVal) {
92: ctx->bsInd = ctx->rdim - 2;
93: ctx->complexInd = ctx->rdim - 1;
94: } else {
95: ctx->bsInd = ctx->rdim - 1;
96: ctx->complexInd = -1;
97: }
102: if (hdf5->horizontal) {
103: PetscInt t;
104: /* support horizontal 1D arrays (MATLAB vectors) - swap meaning of blocks and entries */
105: t = ctx->lenInd;
106: ctx->lenInd = ctx->bsInd;
107: ctx->bsInd = t;
108: }
110: /* Get block size */
111: ctx->dim2 = PETSC_FALSE;
112: if (ctx->lenInd == ctx->bsInd) {
113: bs = 1; /* support vectors stored as 1D array */
114: } else {
115: bs = (PetscInt)ctx->dims[ctx->bsInd];
116: if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
117: }
119: /* Get global size */
120: len = ctx->dims[ctx->lenInd];
121: N = (PetscInt)len * bs;
123: /* Set global size, blocksize and type if not yet set */
124: if (map->bs < 0) {
125: PetscLayoutSetBlockSize(map, bs);
127: if (map->N < 0) {
128: PetscLayoutSetSize(map, N);
130: if (setup) PetscLayoutSetUp(map);
131: return 0;
132: }
134: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
135: {
136: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
137: hsize_t *count, *offset;
138: PetscInt bs, n, low;
139: int i;
141: /* Compute local size and ownership range */
142: PetscLayoutSetUp(map);
143: PetscLayoutGetBlockSize(map, &bs);
144: PetscLayoutGetLocalSize(map, &n);
145: PetscLayoutGetRange(map, &low, NULL);
147: /* Each process defines a dataset and reads it from the hyperslab in the file */
148: PetscMalloc2(ctx->rdim, &count, ctx->rdim, &offset);
149: for (i = 0; i < ctx->rdim; i++) {
150: /* By default, select all entries with no offset */
151: offset[i] = 0;
152: count[i] = ctx->dims[i];
153: }
154: if (hdf5->timestepping) {
155: count[0] = 1;
156: offset[0] = hdf5->timestep;
157: }
158: {
159: PetscHDF5IntCast(n / bs, &count[ctx->lenInd]);
160: PetscHDF5IntCast(low / bs, &offset[ctx->lenInd]);
161: }
162: PetscCallHDF5Return(*memspace, H5Screate_simple, (ctx->rdim, count, NULL));
163: PetscCallHDF5(H5Sselect_hyperslab, (ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
164: PetscFree2(count, offset);
165: return 0;
166: }
168: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
169: {
170: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
172: PetscCallHDF5(H5Dread, (h->dataset, datatype, memspace, h->dataspace, hdf5->dxpl_id, arr));
173: return 0;
174: }
176: /*@C
177: PetscViewerHDF5Load - Read a raw array from the `PETSCVIEWERHDF5` dataset.
179: Input Parameters:
180: + viewer - The HDF5 viewer
181: . name - The dataset name
182: - datatype - The HDF5 datatype of the items in the dataset
184: Input/Output Parameter:
185: . map - The layout which specifies array partitioning, on output the
186: set up layout (with global size and blocksize according to dataset)
188: Output Parameter:
189: . newarr - The partitioned array, a memory image of the given dataset
191: Level: developer
193: Notes:
194: This is intended mainly for internal use; users should use higher level routines such as `ISLoad()`, `VecLoad()`, `DMLoad()`.
196: The array is partitioned according to the given `PetscLayout` which is converted to an HDF5 hyperslab.
198: This name is relative to the current group returned by `PetscViewerHDF5OpenGroup()`.
200: Fortran Note:
201: This routine is not available in Fortran.
203: .seealso: `PetscViewer`, `PETSCVIEWERHDF5`, `PetscViewerHDF5Open()`, `PetscViewerHDF5PushGroup()`, `PetscViewerHDF5OpenGroup()`, `PetscViewerHDF5ReadSizes()`, `VecLoad()`, `ISLoad()`
204: @*/
205: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
206: {
207: PetscBool has;
208: char *group;
209: HDF5ReadCtx h = NULL;
210: hid_t memspace = 0;
211: size_t unitsize;
212: void *arr;
214: PetscViewerHDF5GetGroup(viewer, NULL, &group);
215: PetscViewerHDF5HasDataset(viewer, name, &has);
217: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
218: #if defined(PETSC_USE_COMPLEX)
219: if (!h->complexVal) {
220: H5T_class_t clazz = H5Tget_class(datatype);
222: }
223: #else
225: #endif
227: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_TRUE, &map);
228: PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace);
230: unitsize = H5Tget_size(datatype);
231: if (h->complexVal) unitsize *= 2;
232: /* unitsize is size_t i.e. always unsigned, so the negative check is pointless? */
234: PetscMalloc(map->n * unitsize, &arr);
236: PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr);
237: PetscCallHDF5(H5Sclose, (memspace));
238: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
239: PetscFree(group);
240: *newarr = arr;
241: return 0;
242: }
244: /*@C
245: PetscViewerHDF5ReadSizes - Read block size and global size of a `Vec` or `IS` stored in an HDF5 file.
247: Input Parameters:
248: + viewer - The HDF5 viewer
249: - name - The dataset name
251: Output Parameters:
252: + bs - block size
253: - N - global size
255: Level: advanced
257: Notes:
258: The dataset is stored as an HDF5 dataspace with 1-4 dimensions in the order
259: 1) # timesteps (optional), 2) # blocks, 3) # elements per block (optional), 4) real and imaginary part (only for complex).
261: The dataset can be stored as a 2D dataspace even if its blocksize is 1; see `PetscViewerHDF5SetBaseDimension2()`.
263: .seealso: `PetscViewer`, `PETSCVIEWERHDF5`, `PetscViewerHDF5Open()`, `VecLoad()`, `ISLoad()`, `VecGetSize()`, `ISGetSize()`, `PetscViewerHDF5SetBaseDimension2()`
264: @*/
265: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
266: {
267: HDF5ReadCtx h = NULL;
268: PetscLayout map = NULL;
271: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
272: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_FALSE, &map);
273: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
274: if (bs) *bs = map->bs;
275: if (N) *N = map->N;
276: PetscLayoutDestroy(&map);
277: return 0;
278: }
280: #endif /* defined(PETSC_HAVE_HDF5) */