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;
19: const char *group;
21: PetscViewerHDF5GetGroup(viewer, &group);
22: PetscViewerHDF5ReadAttribute(viewer,name,"timestepping",PETSC_BOOL,×tepping,×tepping);
24: return 0;
25: }
27: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
28: {
29: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
30: HDF5ReadCtx h=NULL;
32: PetscViewerHDF5CheckTimestepping_Internal(viewer, name);
33: PetscNew(&h);
34: PetscViewerHDF5OpenGroup(viewer, &h->file, &h->group);
35: PetscStackCallHDF5Return(h->dataset,H5Dopen2,(h->group, name, H5P_DEFAULT));
36: PetscStackCallHDF5Return(h->dataspace,H5Dget_space,(h->dataset));
37: PetscViewerHDF5ReadAttribute(viewer,name,"complex",PETSC_BOOL,&h->complexVal,&h->complexVal);
38: if (!hdf5->horizontal) {
39: /* MATLAB stores column vectors horizontally */
40: PetscViewerHDF5HasAttribute(viewer,name,"MATLAB_class",&hdf5->horizontal);
41: }
42: *ctx = h;
43: return 0;
44: }
46: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
47: {
48: HDF5ReadCtx h;
50: h = *ctx;
51: PetscStackCallHDF5(H5Gclose,(h->group));
52: PetscStackCallHDF5(H5Sclose,(h->dataspace));
53: PetscStackCallHDF5(H5Dclose,(h->dataset));
54: PetscFree((*ctx)->dims);
55: PetscFree(*ctx);
56: return 0;
57: }
59: static PetscErrorCode PetscViewerHDF5ReadSizes_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscBool setup, PetscLayout *map_)
60: {
61: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
62: PetscInt bs, len, N;
63: PetscLayout map;
65: if (!(*map_)) {
66: PetscLayoutCreate(PetscObjectComm((PetscObject)viewer),map_);
67: }
68: map = *map_;
70: /* Get actual number of dimensions in dataset */
71: PetscStackCallHDF5Return(ctx->rdim,H5Sget_simple_extent_dims,(ctx->dataspace, NULL, NULL));
72: PetscMalloc1(ctx->rdim, &ctx->dims);
73: PetscStackCallHDF5Return(ctx->rdim,H5Sget_simple_extent_dims,(ctx->dataspace, ctx->dims, NULL));
75: /*
76: Dimensions are in this order:
77: [0] timesteps (optional)
78: [lenInd] entries (numbers or blocks)
79: ...
80: [bsInd] entries of blocks (optional)
81: [bsInd+1] real & imaginary part (optional)
82: = rdim-1
83: */
85: /* Get entries dimension index */
86: ctx->lenInd = 0;
87: if (hdf5->timestepping) ++ctx->lenInd;
89: /* Get block dimension index */
90: if (ctx->complexVal) {
91: ctx->bsInd = ctx->rdim-2;
92: ctx->complexInd = ctx->rdim-1;
93: } else {
94: ctx->bsInd = ctx->rdim-1;
95: ctx->complexInd = -1;
96: }
101: if (hdf5->horizontal) {
102: PetscInt t;
103: /* support horizontal 1D arrays (MATLAB vectors) - swap meaning of blocks and entries */
104: t = ctx->lenInd; ctx->lenInd = ctx->bsInd; ctx->bsInd = t;
105: }
107: /* Get block size */
108: ctx->dim2 = PETSC_FALSE;
109: if (ctx->lenInd == ctx->bsInd) {
110: bs = 1; /* support vectors stored as 1D array */
111: } else {
112: bs = (PetscInt) ctx->dims[ctx->bsInd];
113: if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
114: }
116: /* Get global size */
117: len = ctx->dims[ctx->lenInd];
118: N = (PetscInt) len*bs;
120: /* Set global size, blocksize and type if not yet set */
121: if (map->bs < 0) {
122: PetscLayoutSetBlockSize(map, bs);
124: if (map->N < 0) {
125: PetscLayoutSetSize(map, N);
127: if (setup) PetscLayoutSetUp(map);
128: return 0;
129: }
131: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
132: {
133: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
134: hsize_t *count, *offset;
135: PetscInt bs, n, low;
136: int i;
138: /* Compute local size and ownership range */
139: PetscLayoutSetUp(map);
140: PetscLayoutGetBlockSize(map, &bs);
141: PetscLayoutGetLocalSize(map, &n);
142: PetscLayoutGetRange(map, &low, NULL);
144: /* Each process defines a dataset and reads it from the hyperslab in the file */
145: PetscMalloc2(ctx->rdim, &count, ctx->rdim, &offset);
146: for (i=0; i<ctx->rdim; i++) {
147: /* By default, select all entries with no offset */
148: offset[i] = 0;
149: count[i] = ctx->dims[i];
150: }
151: if (hdf5->timestepping) {
152: count[0] = 1;
153: offset[0] = hdf5->timestep;
154: }
155: {
156: PetscHDF5IntCast(n/bs, &count[ctx->lenInd]);
157: PetscHDF5IntCast(low/bs, &offset[ctx->lenInd]);
158: }
159: PetscStackCallHDF5Return(*memspace,H5Screate_simple,(ctx->rdim, count, NULL));
160: PetscStackCallHDF5(H5Sselect_hyperslab,(ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
161: PetscFree2(count, offset);
162: return 0;
163: }
165: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
166: {
167: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
169: PetscStackCallHDF5(H5Dread,(h->dataset, datatype, memspace, h->dataspace, hdf5->dxpl_id, arr));
170: return 0;
171: }
173: /*@C
174: PetscViewerHDF5Load - Read a raw array from the HDF5 dataset.
176: Input Parameters:
177: + viewer - The HDF5 viewer
178: . name - The dataset name
179: - datatype - The HDF5 datatype of the items in the dataset
181: Input/Output Parameter:
182: . map - The layout which specifies array partitioning, on output the
183: set up layout (with global size and blocksize according to dataset)
185: Output Parameter:
186: . newarr - The partitioned array, a memory image of the given dataset
188: Level: developer
190: Notes:
191: This is intended mainly for internal use; users should use higher level routines such as ISLoad(), VecLoad(), DMLoad().
192: The array is partitioned according to the given PetscLayout which is converted to an HDF5 hyperslab.
193: This name is relative to the current group returned by PetscViewerHDF5OpenGroup().
195: Fortran Notes:
196: This routine is not available in Fortran.
198: .seealso PetscViewerHDF5Open(), PetscViewerHDF5PushGroup(), PetscViewerHDF5OpenGroup(), PetscViewerHDF5ReadSizes(), VecLoad(), ISLoad()
199: @*/
200: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
201: {
202: PetscBool has;
203: const char *group;
204: HDF5ReadCtx h=NULL;
205: hid_t memspace=0;
206: size_t unitsize;
207: void *arr;
209: PetscViewerHDF5GetGroup(viewer, &group);
210: PetscViewerHDF5HasDataset(viewer, name, &has);
212: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
213: #if defined(PETSC_USE_COMPLEX)
214: if (!h->complexVal) {
215: H5T_class_t clazz = H5Tget_class(datatype);
217: }
218: #else
220: #endif
222: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_TRUE, &map);
223: PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace);
225: unitsize = H5Tget_size(datatype);
226: if (h->complexVal) unitsize *= 2;
227: /* unitsize is size_t i.e. always unsigned, so the negative check is pointless? */
229: PetscMalloc(map->n*unitsize, &arr);
231: PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr);
232: PetscStackCallHDF5(H5Sclose,(memspace));
233: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
234: *newarr = arr;
235: return 0;
236: }
238: /*@C
239: PetscViewerHDF5ReadSizes - Read block size and global size of a vector (Vec or IS) stored in an HDF5 file.
241: Input Parameters:
242: + viewer - The HDF5 viewer
243: - name - The dataset name
245: Output Parameters:
246: + bs - block size
247: - N - global size
249: Notes:
250: The dataset is stored as an HDF5 dataspace with 1-4 dimensions in the order
251: 1) # timesteps (optional), 2) # blocks, 3) # elements per block (optional), 4) real and imaginary part (only for complex).
253: The dataset can be stored as a 2D dataspace even if its blocksize is 1; see PetscViewerHDF5SetBaseDimension2().
255: Level: advanced
257: .seealso: PetscViewerHDF5Open(), VecLoad(), ISLoad(), VecGetSize(), ISGetSize(), PetscViewerHDF5SetBaseDimension2()
258: @*/
259: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
260: {
261: HDF5ReadCtx h=NULL;
262: PetscLayout map=NULL;
265: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
266: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_FALSE, &map);
267: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
268: if (bs) *bs = map->bs;
269: if (N) *N = map->N;
270: PetscLayoutDestroy(&map);
271: return 0;
272: }
274: #endif /* defined(PETSC_HAVE_HDF5) */