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: PetscFunctionBegin;
21: PetscCall(PetscViewerHDF5ReadAttribute(viewer, name, "timestepping", PETSC_BOOL, &hdf5->defTimestepping, ×tepping));
22: if (timestepping != hdf5->timestepping) {
23: char *group;
25: PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &group));
26: 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]);
27: }
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
32: {
33: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
34: HDF5ReadCtx h = NULL;
36: PetscFunctionBegin;
37: PetscCall(PetscViewerHDF5CheckTimestepping_Internal(viewer, name));
38: PetscCall(PetscNew(&h));
39: PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &h->file, &h->group));
40: PetscCallHDF5Return(h->dataset, H5Dopen2, (h->group, name, H5P_DEFAULT));
41: PetscCallHDF5Return(h->dataspace, H5Dget_space, (h->dataset));
42: PetscCall(PetscViewerHDF5ReadAttribute(viewer, name, "complex", PETSC_BOOL, &h->complexVal, &h->complexVal));
43: if (!hdf5->horizontal) {
44: /* MATLAB stores column vectors horizontally */
45: PetscCall(PetscViewerHDF5HasAttribute(viewer, name, "MATLAB_class", &hdf5->horizontal));
46: }
47: *ctx = h;
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
52: {
53: HDF5ReadCtx h;
55: PetscFunctionBegin;
56: h = *ctx;
57: PetscCallHDF5(H5Gclose, (h->group));
58: PetscCallHDF5(H5Sclose, (h->dataspace));
59: PetscCallHDF5(H5Dclose, (h->dataset));
60: PetscCall(PetscFree((*ctx)->dims));
61: PetscCall(PetscFree(*ctx));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode PetscViewerHDF5ReadSizes_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscBool setup, PetscLayout *map_)
66: {
67: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
68: PetscInt bs, len, N;
69: PetscLayout map;
71: PetscFunctionBegin;
72: if (!*map_) PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)viewer), map_));
73: map = *map_;
75: /* Get actual number of dimensions in dataset */
76: PetscCallHDF5Return(ctx->rdim, H5Sget_simple_extent_dims, (ctx->dataspace, NULL, NULL));
77: PetscCall(PetscMalloc1(ctx->rdim, &ctx->dims));
78: PetscCallHDF5Return(ctx->rdim, H5Sget_simple_extent_dims, (ctx->dataspace, ctx->dims, NULL));
80: /*
81: Dimensions are in this order:
82: [0] timesteps (optional)
83: [lenInd] entries (numbers or blocks)
84: ...
85: [bsInd] entries of blocks (optional)
86: [bsInd+1] real & imaginary part (optional)
87: = rdim-1
88: */
90: /* Get entries dimension index */
91: ctx->lenInd = 0;
92: if (hdf5->timestepping) ++ctx->lenInd;
94: /* Get block dimension index */
95: if (ctx->complexVal) {
96: ctx->bsInd = ctx->rdim - 2;
97: ctx->complexInd = ctx->rdim - 1;
98: } else {
99: ctx->bsInd = ctx->rdim - 1;
100: ctx->complexInd = -1;
101: }
102: PetscCheck(ctx->lenInd <= ctx->bsInd, PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Calculated block dimension index = %d < %d = length dimension index.", ctx->bsInd, ctx->lenInd);
103: PetscCheck(ctx->bsInd <= ctx->rdim - 1, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Calculated block dimension index = %d > %d = total number of dimensions - 1.", ctx->bsInd, ctx->rdim - 1);
104: PetscCheck(!ctx->complexVal || ctx->dims[ctx->complexInd] == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Complex numbers must have exactly 2 parts (%" PRIuHSIZE ")", ctx->dims[ctx->complexInd]);
106: if (hdf5->horizontal) {
107: PetscInt t;
108: /* support horizontal 1D arrays (MATLAB vectors) - swap meaning of blocks and entries */
109: t = ctx->lenInd;
110: ctx->lenInd = ctx->bsInd;
111: ctx->bsInd = t;
112: }
114: /* Get block size */
115: ctx->dim2 = PETSC_FALSE;
116: if (ctx->lenInd == ctx->bsInd) {
117: bs = 1; /* support vectors stored as 1D array */
118: } else {
119: bs = (PetscInt)ctx->dims[ctx->bsInd];
120: if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
121: }
123: /* Get global size */
124: len = ctx->dims[ctx->lenInd];
125: N = (PetscInt)len * bs;
127: /* Set global size, blocksize and type if not yet set */
128: if (map->bs < 0) {
129: PetscCall(PetscLayoutSetBlockSize(map, bs));
130: } else PetscCheck(map->bs == bs, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Block size of array in file is %" PetscInt_FMT ", not %" PetscInt_FMT " as expected", bs, map->bs);
131: if (map->N < 0) {
132: PetscCall(PetscLayoutSetSize(map, N));
133: } else PetscCheck(map->N == N, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Global size of array in file is %" PetscInt_FMT ", not %" PetscInt_FMT " as expected", N, map->N);
134: if (setup) PetscCall(PetscLayoutSetUp(map));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
139: {
140: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
141: hsize_t *count, *offset;
142: PetscInt bs, n, low;
143: int i;
145: PetscFunctionBegin;
146: /* Compute local size and ownership range */
147: PetscCall(PetscLayoutSetUp(map));
148: PetscCall(PetscLayoutGetBlockSize(map, &bs));
149: PetscCall(PetscLayoutGetLocalSize(map, &n));
150: PetscCall(PetscLayoutGetRange(map, &low, NULL));
152: /* Each process defines a dataset and reads it from the hyperslab in the file */
153: PetscCall(PetscMalloc2(ctx->rdim, &count, ctx->rdim, &offset));
154: for (i = 0; i < ctx->rdim; i++) {
155: /* By default, select all entries with no offset */
156: offset[i] = 0;
157: count[i] = ctx->dims[i];
158: }
159: if (hdf5->timestepping) {
160: count[0] = 1;
161: offset[0] = hdf5->timestep;
162: }
163: {
164: PetscCall(PetscHDF5IntCast(n / bs, &count[ctx->lenInd]));
165: PetscCall(PetscHDF5IntCast(low / bs, &offset[ctx->lenInd]));
166: }
167: PetscCallHDF5Return(*memspace, H5Screate_simple, (ctx->rdim, count, NULL));
168: PetscCallHDF5(H5Sselect_hyperslab, (ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
169: PetscCall(PetscFree2(count, offset));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
174: {
175: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
177: PetscFunctionBegin;
178: PetscCallHDF5(H5Dread, (h->dataset, datatype, memspace, h->dataspace, hdf5->dxpl_id, arr));
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: /*@C
183: PetscViewerHDF5Load - Read a raw array from the `PETSCVIEWERHDF5` dataset in parallel
185: Collective; No Fortran Support
187: Input Parameters:
188: + viewer - The `PETSCVIEWERHDF5` viewer
189: . name - The dataset name
190: - datatype - The HDF5 datatype of the items in the dataset
192: Input/Output Parameter:
193: . map - The layout which specifies array partitioning, on output the
194: set up layout (with global size and blocksize according to dataset)
196: Output Parameter:
197: . newarr - The partitioned array, a memory image of the given dataset
199: Level: developer
201: Notes:
202: This is intended mainly for internal use; users should use higher level routines such as `ISLoad()`, `VecLoad()`, `DMLoad()`.
204: The array is partitioned according to the given `PetscLayout` which is converted to an HDF5 hyperslab.
206: This name is relative to the current group returned by `PetscViewerHDF5OpenGroup()`.
208: .seealso: `PetscViewer`, `PETSCVIEWERHDF5`, `PetscViewerHDF5Open()`, `PetscViewerHDF5PushGroup()`, `PetscViewerHDF5OpenGroup()`, `PetscViewerHDF5ReadSizes()`,
209: `VecLoad()`, `ISLoad()`, `PetscLayout`
210: @*/
211: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
212: {
213: PetscBool has;
214: char *group;
215: HDF5ReadCtx h = NULL;
216: hid_t memspace = 0;
217: size_t unitsize;
218: void *arr;
220: PetscFunctionBegin;
221: PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &group));
222: PetscCall(PetscViewerHDF5HasDataset(viewer, name, &has));
223: PetscCheck(has, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Object (dataset) \"%s\" not stored in group %s", name, group);
224: PetscCall(PetscViewerHDF5ReadInitialize_Private(viewer, name, &h));
225: #if defined(PETSC_USE_COMPLEX)
226: if (!h->complexVal) {
227: H5T_class_t clazz = H5Tget_class(datatype);
228: PetscCheck(clazz != H5T_FLOAT, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Dataset %s/%s is marked as real but PETSc is configured for complex scalars. The conversion is not yet implemented. Configure with --with-scalar-type=real to read this dataset", group ? group : "", name);
229: }
230: #else
231: PetscCheck(!h->complexVal, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Dataset %s/%s is marked as complex but PETSc is configured for real scalars. Configure with --with-scalar-type=complex to read this dataset", group, name);
232: #endif
234: PetscCall(PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_TRUE, &map));
235: PetscCall(PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace));
237: unitsize = H5Tget_size(datatype);
238: if (h->complexVal) unitsize *= 2;
239: /* unitsize is size_t i.e. always unsigned, so the negative check is pointless? */
240: PetscCheck(unitsize > 0 && unitsize <= PetscMax(sizeof(PetscInt), sizeof(PetscScalar)), PETSC_COMM_SELF, PETSC_ERR_LIB, "Sanity check failed: HDF5 function H5Tget_size(datatype) returned suspicious value %zu", unitsize);
241: PetscCall(PetscMalloc(map->n * unitsize, &arr));
243: PetscCall(PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr));
244: PetscCallHDF5(H5Sclose, (memspace));
245: PetscCall(PetscViewerHDF5ReadFinalize_Private(viewer, &h));
246: PetscCall(PetscFree(group));
247: *newarr = arr;
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: /*@C
252: PetscViewerHDF5ReadSizes - Read block size and global size of a `Vec` or `IS` stored in an HDF5 file.
254: Input Parameters:
255: + viewer - The `PETSCVIEWERHDF5` viewer
256: - name - The dataset name
258: Output Parameters:
259: + bs - block size
260: - N - global size
262: Level: advanced
264: Notes:
265: The dataset is stored as an HDF5 dataspace with 1-4 dimensions in the order
266: 1) # timesteps (optional), 2) # blocks, 3) # elements per block (optional), 4) real and imaginary part (only for complex).
268: The dataset can be stored as a 2D dataspace even if its blocksize is 1; see `PetscViewerHDF5SetBaseDimension2()`.
270: .seealso: `PetscViewer`, `PETSCVIEWERHDF5`, `PetscViewerHDF5Open()`, `VecLoad()`, `ISLoad()`, `VecGetSize()`, `ISGetSize()`, `PetscViewerHDF5SetBaseDimension2()`
271: @*/
272: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
273: {
274: HDF5ReadCtx h = NULL;
275: PetscLayout map = NULL;
277: PetscFunctionBegin;
279: PetscCall(PetscViewerHDF5ReadInitialize_Private(viewer, name, &h));
280: PetscCall(PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_FALSE, &map));
281: PetscCall(PetscViewerHDF5ReadFinalize_Private(viewer, &h));
282: if (bs) *bs = map->bs;
283: if (N) *N = map->N;
284: PetscCall(PetscLayoutDestroy(&map));
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: #endif /* defined(PETSC_HAVE_HDF5) */