Actual source code: hdf5io.c
petsc-3.13.6 2020-09-29
1: #include <petsc/private/viewerimpl.h>
2: #include <petsc/private/viewerhdf5impl.h>
3: #include <petsclayouthdf5.h>
4: #include <petscis.h>
6: #if defined(PETSC_HAVE_HDF5)
8: struct _n_HDF5ReadCtx {
9: hid_t file, group, dataset, dataspace;
10: PetscInt timestep;
11: int lenInd, bsInd, rdim;
12: hsize_t *dims;
13: PetscBool complexVal, dim2, horizontal;
14: };
15: typedef struct _n_HDF5ReadCtx* HDF5ReadCtx;
17: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
18: {
19: HDF5ReadCtx h=NULL;
23: PetscNew(&h);
24: PetscViewerHDF5OpenGroup(viewer, &h->file, &h->group);
25: PetscStackCallHDF5Return(h->dataset,H5Dopen2,(h->group, name, H5P_DEFAULT));
26: PetscStackCallHDF5Return(h->dataspace,H5Dget_space,(h->dataset));
27: PetscViewerHDF5GetTimestep(viewer, &h->timestep);
28: PetscViewerHDF5HasAttribute(viewer,name,"complex",&h->complexVal);
29: if (h->complexVal) {PetscViewerHDF5ReadAttribute(viewer,name,"complex",PETSC_BOOL,&h->complexVal);}
30: /* MATLAB stores column vectors horizontally */
31: PetscViewerHDF5HasAttribute(viewer,name,"MATLAB_class",&h->horizontal);
32: *ctx = h;
33: return(0);
34: }
36: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
37: {
38: HDF5ReadCtx h;
42: h = *ctx;
43: PetscStackCallHDF5(H5Gclose,(h->group));
44: PetscStackCallHDF5(H5Sclose,(h->dataspace));
45: PetscStackCallHDF5(H5Dclose,(h->dataset));
46: PetscFree((*ctx)->dims);
47: PetscFree(*ctx);
48: return(0);
49: }
51: static PetscErrorCode PetscViewerHDF5ReadSizes_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscBool setup, PetscLayout *map_)
52: {
53: PetscInt bs, len, N;
54: PetscLayout map;
58: if (!(*map_)) {
59: PetscLayoutCreate(PetscObjectComm((PetscObject)viewer),map_);
60: }
61: map = *map_;
63: /* Get actual number of dimensions in dataset */
64: PetscStackCallHDF5Return(ctx->rdim,H5Sget_simple_extent_dims,(ctx->dataspace, NULL, NULL));
65: PetscMalloc1(ctx->rdim, &ctx->dims);
66: PetscStackCallHDF5Return(ctx->rdim,H5Sget_simple_extent_dims,(ctx->dataspace, ctx->dims, NULL));
68: /*
69: Dimensions are in this order:
70: [0] timesteps (optional)
71: [lenInd] entries (numbers or blocks)
72: ...
73: [bsInd] entries of blocks (optional)
74: [bsInd+1] real & imaginary part (optional)
75: = rdim-1
76: */
78: /* Get entries dimension index */
79: ctx->lenInd = 0;
80: if (ctx->timestep >= 0) ++ctx->lenInd;
82: /* Get block dimension index */
83: ctx->bsInd = ctx->rdim-1;
84: if (ctx->complexVal) --ctx->bsInd;
85: 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);
86: 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);
88: /* Get block size */
89: ctx->dim2 = PETSC_FALSE;
90: if (ctx->lenInd == ctx->bsInd) {
91: bs = 1; /* support vectors stored as 1D array */
92: } else {
93: bs = (PetscInt) ctx->dims[ctx->bsInd];
94: if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
95: }
97: /* Get global size */
98: len = ctx->dims[ctx->lenInd];
99: if (ctx->horizontal) {
100: PetscInt t;
101: /* support horizontal 1D arrays (MATLAB vectors) - swap meaning of blocks and entries */
102: if (ctx->complexVal) SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Complex and horizontal at the same time not allowed.");
103: t = len; len = bs; bs = t;
104: t = ctx->lenInd; ctx->lenInd = ctx->bsInd; ctx->bsInd = t;
105: }
106: N = (PetscInt) len*bs;
108: /* Set global size, blocksize and type if not yet set */
109: if (map->bs < 0) {
110: PetscLayoutSetBlockSize(map, bs);
111: } 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);
112: if (map->N < 0) {
113: PetscLayoutSetSize(map, N);
114: } 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);
115: if (setup) {PetscLayoutSetUp(map);}
116: return(0);
117: }
119: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
120: {
121: hsize_t *count, *offset;
122: PetscInt bs, n, low;
123: int i;
127: /* Compute local size and ownership range */
128: PetscLayoutSetUp(map);
129: PetscLayoutGetBlockSize(map, &bs);
130: PetscLayoutGetLocalSize(map, &n);
131: PetscLayoutGetRange(map, &low, NULL);
133: /* Each process defines a dataset and reads it from the hyperslab in the file */
134: PetscMalloc2(ctx->rdim, &count, ctx->rdim, &offset);
135: for (i=0; i<ctx->rdim; i++) {
136: /* By default, select all entries with no offset */
137: offset[i] = 0;
138: count[i] = ctx->dims[i];
139: }
140: if (ctx->timestep >= 0) {
141: count[0] = 1;
142: offset[0] = ctx->timestep;
143: }
144: {
145: PetscHDF5IntCast(n/bs, &count[ctx->lenInd]);
146: PetscHDF5IntCast(low/bs, &offset[ctx->lenInd]);
147: }
148: if (ctx->complexVal && count[ctx->bsInd+1] != 2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Complex numbers must have exactly 2 parts (%D)",count[ctx->bsInd+1]);
149: PetscStackCallHDF5Return(*memspace,H5Screate_simple,(ctx->rdim, count, NULL));
150: PetscStackCallHDF5(H5Sselect_hyperslab,(ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
151: PetscFree2(count, offset);
152: return(0);
153: }
155: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
156: {
157: PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
160: PetscStackCallHDF5(H5Dread,(h->dataset, datatype, memspace, h->dataspace, hdf5->dxpl_id, arr));
161: return(0);
162: }
164: /*@C
165: PetscViewerHDF5Load - Read a raw array from the HDF5 dataset.
167: Input Parameters:
168: + viewer - The HDF5 viewer
169: . name - The dataset name
170: . map - The layout which specifies array partitioning
171: - datatype - The HDF5 datatype of the items in the dataset
173: Output Parameter:
174: + map - The set up layout (with global size and blocksize according to dataset)
175: - newarr - The partitioned array, a memory image of the given dataset
177: Level: developer
179: Notes:
180: This is intended mainly for internal use; users should use higher level routines such as ISLoad(), VecLoad(), DMLoad().
181: The array is partitioned according to the given PetscLayout which is converted to an HDF5 hyperslab.
182: This name is relative to the current group returned by PetscViewerHDF5OpenGroup().
184: Fortran Notes:
185: This routine is not available in Fortran.
187: .seealso PetscViewerHDF5Open(), PetscViewerHDF5PushGroup(), PetscViewerHDF5OpenGroup(), PetscViewerHDF5ReadSizes(), VecLoad(), ISLoad()
188: @*/
189: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
190: {
191: HDF5ReadCtx h=NULL;
192: hid_t memspace=0;
193: size_t unitsize;
194: void *arr;
195: PetscErrorCode ierr;
198: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
199: #if defined(PETSC_USE_COMPLEX)
200: if (!h->complexVal) {
201: H5T_class_t clazz = H5Tget_class(datatype);
202: if (clazz == H5T_FLOAT) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"File contains real numbers but PETSc is configured for complex. The conversion is not yet implemented. Configure with --with-scalar-type=real.");
203: }
204: #else
205: if (h->complexVal) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"File contains complex numbers but PETSc not configured for them. Configure with --with-scalar-type=complex.");
206: #endif
208: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_TRUE, &map);
209: PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace);
211: unitsize = H5Tget_size(datatype);
212: if (h->complexVal) unitsize *= 2;
213: 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);
214: PetscMalloc(map->n*unitsize, &arr);
216: PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr);
217: PetscStackCallHDF5(H5Sclose,(memspace));
218: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
219: *newarr = arr;
220: return(0);
221: }
223: /*@C
224: PetscViewerHDF5ReadSizes - Read block size and global size of a vector (Vec or IS) stored in an HDF5 file.
226: Input Parameters:
227: + viewer - The HDF5 viewer
228: - name - The dataset name
230: Output Parameter:
231: + bs - block size
232: - N - global size
234: Notes:
235: The dataset is stored as an HDF5 dataspace with 1-4 dimensions in the order
236: 1) # timesteps (optional), 2) # blocks, 3) # elements per block (optional), 4) real and imaginary part (only for complex).
238: The dataset can be stored as a 2D dataspace even if its blocksize is 1; see PetscViewerHDF5SetBaseDimension2().
240: Level: advanced
242: .seealso: PetscViewerHDF5Open(), VecLoad(), ISLoad(), VecGetSize(), ISGetSize(), PetscViewerHDF5SetBaseDimension2()
243: @*/
244: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
245: {
246: HDF5ReadCtx h=NULL;
247: PetscLayout map=NULL;
252: PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
253: PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_FALSE, &map);
254: PetscViewerHDF5ReadFinalize_Private(viewer, &h);
255: if (bs) *bs = map->bs;
256: if (N) *N = map->N;
257: PetscLayoutDestroy(&map);
258: return(0);
259: }
261: #endif /* defined(PETSC_HAVE_HDF5) */