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, &timestepping, &timestepping);
 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) */