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;
20: PetscErrorCode ierr;
23: PetscViewerHDF5GetGroup(viewer, &group);
24: PetscViewerHDF5ReadAttribute(viewer,name,"timestepping",PETSC_BOOL,×tepping,×tepping);
25: if (timestepping != hdf5->timestepping) SETERRQ4(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: return(0);
27: }
29: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
30: {
31: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
32: HDF5ReadCtx h=NULL;
33: PetscErrorCode ierr;
36: PetscViewerHDF5CheckTimestepping_Internal(viewer, name);
37: PetscNew(&h);
38: PetscViewerHDF5OpenGroup(viewer, &h->file, &h->group);
39: PetscStackCallHDF5Return(h->dataset,H5Dopen2,(h->group, name, H5P_DEFAULT));
40: PetscStackCallHDF5Return(h->dataspace,H5Dget_space,(h->dataset));
41: PetscViewerHDF5ReadAttribute(viewer,name,"complex",PETSC_BOOL,&h->complexVal,&h->complexVal);
42: if (!hdf5->horizontal) {
43: /* MATLAB stores column vectors horizontally */
44: PetscViewerHDF5HasAttribute(viewer,name,"MATLAB_class",&hdf5->horizontal);
45: }
46: *ctx = h;
47: return(0);
48: }
50: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
51: {
52: HDF5ReadCtx h;
56: h = *ctx;
57: PetscStackCallHDF5(H5Gclose,(h->group));
58: PetscStackCallHDF5(H5Sclose,(h->dataspace));
59: PetscStackCallHDF5(H5Dclose,(h->dataset));
60: PetscFree((*ctx)->dims);
61: PetscFree(*ctx);
62: return(0);
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;
70: PetscErrorCode ierr;
73: if (!(*map_)) {
74: PetscLayoutCreate(PetscObjectComm((PetscObject)viewer),map_);
75: }
76: map = *map_;
78: /* Get actual number of dimensions in dataset */
79: PetscStackCallHDF5Return(ctx->rdim,H5Sget_simple_extent_dims,(ctx->dataspace, NULL, NULL));
80: PetscMalloc1(ctx->rdim, &ctx->dims);
81: PetscStackCallHDF5Return(ctx->rdim,H5Sget_simple_extent_dims,(ctx->dataspace, ctx->dims, NULL));
83: /*
84: Dimensions are in this order:
85: [0] timesteps (optional)
86: [lenInd] entries (numbers or blocks)
87: ...
88: [bsInd] entries of blocks (optional)
89: [bsInd+1] real & imaginary part (optional)
90: = rdim-1
91: */
93: /* Get entries dimension index */
94: ctx->lenInd = 0;
95: if (hdf5->timestepping) ++ctx->lenInd;
97: /* Get block dimension index */
98: if (ctx->complexVal) {
99: ctx->bsInd = ctx->rdim-2;
100: ctx->complexInd = ctx->rdim-1;
101: } else {
102: ctx->bsInd = ctx->rdim-1;
103: ctx->complexInd = -1;
104: }
105: if (ctx->lenInd > ctx->bsInd) SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Calculated block dimension index = %D < %D = length dimension index.",ctx->bsInd,ctx->lenInd);
106: if (ctx->bsInd > ctx->rdim - 1) SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Calculated block dimension index = %D > %D = total number of dimensions - 1.",ctx->bsInd,ctx->rdim-1);
107: if (ctx->complexVal && ctx->dims[ctx->complexInd] != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Complex numbers must have exactly 2 parts (%D)",ctx->dims[ctx->complexInd]);
109: if (hdf5->horizontal) {
110: PetscInt t;
111: /* support horizontal 1D arrays (MATLAB vectors) - swap meaning of blocks and entries */
112: t = ctx->lenInd; ctx->lenInd = ctx->bsInd; ctx->bsInd = t;
113: }
115: /* Get block size */
116: ctx->dim2 = PETSC_FALSE;
117: if (ctx->lenInd == ctx->bsInd) {
118: bs = 1; /* support vectors stored as 1D array */
119: } else {
120: bs = (PetscInt) ctx->dims[ctx->bsInd];
121: if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
122: }
124: /* Get global size */
125: len = ctx->dims[ctx->lenInd];
126: N = (PetscInt) len*bs;
128: /* Set global size, blocksize and type if not yet set */
129: if (map->bs < 0) {
130: PetscLayoutSetBlockSize(map, bs);
131: } else if (map->bs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Block size of array in file is %D, not %D as expected",bs,map->bs);
132: if (map->N < 0) {
133: PetscLayoutSetSize(map, N);
134: } else if (map->N != N) SETERRQ2(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Global size of array in file is %D, not %D as expected",N,map->N);
135: if (setup) {PetscLayoutSetUp(map);}
136: return(0);
137: }
139: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
140: {
141: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
142: hsize_t *count, *offset;
143: PetscInt bs, n, low;
144: int i;
145: PetscErrorCode ierr;
148: /* Compute local size and ownership range */
149: PetscLayoutSetUp(map);
150: PetscLayoutGetBlockSize(map, &bs);
151: PetscLayoutGetLocalSize(map, &n);
152: PetscLayoutGetRange(map, &low, NULL);
154: /* Each process defines a dataset and reads it from the hyperslab in the file */
155: PetscMalloc2(ctx->rdim, &count, ctx->rdim, &offset);
156: for (i=0; i<ctx->rdim; i++) {
157: /* By default, select all entries with no offset */
158: offset[i] = 0;
159: count[i] = ctx->dims[i];
160: }
161: if (hdf5->timestepping) {
162: count[0] = 1;
163: offset[0] = hdf5->timestep;
164: }
165: {
166: PetscHDF5IntCast(n/bs, &count[ctx->lenInd]);
167: PetscHDF5IntCast(low/bs, &offset[ctx->lenInd]);
168: }
169: PetscStackCallHDF5Return(*memspace,H5Screate_simple,(ctx->rdim, count, NULL));
170: PetscStackCallHDF5(H5Sselect_hyperslab,(ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
171: PetscFree2(count, offset);
172: return(0);
173: }
175: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
176: {
177: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
180: PetscStackCallHDF5(H5Dread,(h->dataset, datatype, memspace, h->dataspace, hdf5->dxpl_id, arr));
181: return(0);
182: }
184: /*@C
185: PetscViewerHDF5Load - Read a raw array from the HDF5 dataset.
187: Input Parameters:
188: + viewer - The HDF5 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().
203: The array is partitioned according to the given PetscLayout which is converted to an HDF5 hyperslab.
204: This name is relative to the current group returned by PetscViewerHDF5OpenGroup().
206: Fortran Notes:
207: This routine is not available in Fortran.
209: .seealso PetscViewerHDF5Open(), PetscViewerHDF5PushGroup(), PetscViewerHDF5OpenGroup(), PetscViewerHDF5ReadSizes(), VecLoad(), ISLoad()
210: @*/
211: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
212: {
213: PetscBool has;
214: const char *group;
215: HDF5ReadCtx h=NULL;
216: hid_t memspace=0;
217: size_t unitsize;
218: void *arr;
219: PetscErrorCode ierr;
222: PetscViewerHDF5GetGroup(viewer, &group);
223: PetscViewerHDF5HasDataset(viewer, name, &has);
224: if (!has) SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Object (dataset) \"%s\" not stored in group %s", name, group ? group : "/");
225: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
226: #if defined(PETSC_USE_COMPLEX)
227: if (!h->complexVal) {
228: H5T_class_t clazz = H5Tget_class(datatype);
229: if (clazz == H5T_FLOAT) SETERRQ2(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);
230: }
231: #else
232: if (h->complexVal) SETERRQ2(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 ? group : "",name);
233: #endif
235: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_TRUE, &map);
236: PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace);
238: unitsize = H5Tget_size(datatype);
239: if (h->complexVal) unitsize *= 2;
240: if (unitsize <= 0 || unitsize > PetscMax(sizeof(PetscInt),sizeof(PetscScalar))) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Sanity check failed: HDF5 function H5Tget_size(datatype) returned suspicious value %D",unitsize);
241: PetscMalloc(map->n*unitsize, &arr);
243: PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr);
244: PetscStackCallHDF5(H5Sclose,(memspace));
245: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
246: *newarr = arr;
247: return(0);
248: }
250: /*@C
251: PetscViewerHDF5ReadSizes - Read block size and global size of a vector (Vec or IS) stored in an HDF5 file.
253: Input Parameters:
254: + viewer - The HDF5 viewer
255: - name - The dataset name
257: Output Parameters:
258: + bs - block size
259: - N - global size
261: Notes:
262: The dataset is stored as an HDF5 dataspace with 1-4 dimensions in the order
263: 1) # timesteps (optional), 2) # blocks, 3) # elements per block (optional), 4) real and imaginary part (only for complex).
265: The dataset can be stored as a 2D dataspace even if its blocksize is 1; see PetscViewerHDF5SetBaseDimension2().
267: Level: advanced
269: .seealso: PetscViewerHDF5Open(), VecLoad(), ISLoad(), VecGetSize(), ISGetSize(), PetscViewerHDF5SetBaseDimension2()
270: @*/
271: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
272: {
273: HDF5ReadCtx h=NULL;
274: PetscLayout map=NULL;
279: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
280: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_FALSE, &map);
281: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
282: if (bs) *bs = map->bs;
283: if (N) *N = map->N;
284: PetscLayoutDestroy(&map);
285: return(0);
286: }
288: #endif /* defined(PETSC_HAVE_HDF5) */