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,&timestepping,&timestepping);
 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: . map      - The layout which specifies array partitioning
191: - datatype - The HDF5 datatype of the items in the dataset

193:   Output Parameter:
194: + map      - The set up layout (with global size and blocksize according to dataset)
195: - newarr   - The partitioned array, a memory image of the given dataset

197:   Level: developer

199:   Notes:
200:   This is intended mainly for internal use; users should use higher level routines such as ISLoad(), VecLoad(), DMLoad().
201:   The array is partitioned according to the given PetscLayout which is converted to an HDF5 hyperslab.
202:   This name is relative to the current group returned by PetscViewerHDF5OpenGroup().

204:   Fortran Notes:
205:   This routine is not available in Fortran.

207: .seealso PetscViewerHDF5Open(), PetscViewerHDF5PushGroup(), PetscViewerHDF5OpenGroup(), PetscViewerHDF5ReadSizes(), VecLoad(), ISLoad()
208: @*/
209: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
210: {
211:   PetscBool       has;
212:   const char      *group;
213:   HDF5ReadCtx     h=NULL;
214:   hid_t           memspace=0;
215:   size_t          unitsize;
216:   void            *arr;
217:   PetscErrorCode  ierr;

220:   PetscViewerHDF5GetGroup(viewer, &group);
221:   PetscViewerHDF5HasDataset(viewer, name, &has);
222:   if (!has) SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Object (dataset) \"%s\" not stored in group %s", name, group ? group : "/");
223:   PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
224: #if defined(PETSC_USE_COMPLEX)
225:   if (!h->complexVal) {
226:     H5T_class_t clazz = H5Tget_class(datatype);
227:     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);
228:   }
229: #else
230:   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);
231: #endif

233:   PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_TRUE, &map);
234:   PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace);

236:   unitsize = H5Tget_size(datatype);
237:   if (h->complexVal) unitsize *= 2;
238:   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);
239:   PetscMalloc(map->n*unitsize, &arr);

241:   PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr);
242:   PetscStackCallHDF5(H5Sclose,(memspace));
243:   PetscViewerHDF5ReadFinalize_Private(viewer, &h);
244:   *newarr = arr;
245:   return(0);
246: }

248: /*@C
249:  PetscViewerHDF5ReadSizes - Read block size and global size of a vector (Vec or IS) stored in an HDF5 file.

251:   Input Parameters:
252: + viewer - The HDF5 viewer
253: - name   - The dataset name

255:   Output Parameter:
256: + bs     - block size
257: - N      - global size

259:   Notes:
260:   The dataset is stored as an HDF5 dataspace with 1-4 dimensions in the order
261:   1) # timesteps (optional), 2) # blocks, 3) # elements per block (optional), 4) real and imaginary part (only for complex).

263:   The dataset can be stored as a 2D dataspace even if its blocksize is 1; see PetscViewerHDF5SetBaseDimension2().

265:   Level: advanced

267: .seealso: PetscViewerHDF5Open(), VecLoad(), ISLoad(), VecGetSize(), ISGetSize(), PetscViewerHDF5SetBaseDimension2()
268: @*/
269: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
270: {
271:   HDF5ReadCtx    h=NULL;
272:   PetscLayout    map=NULL;

277:   PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
278:   PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_FALSE, &map);
279:   PetscViewerHDF5ReadFinalize_Private(viewer, &h);
280:   if (bs) *bs = map->bs;
281:   if (N) *N = map->N;
282:   PetscLayoutDestroy(&map);
283:   return(0);
284: }

286: #endif /* defined(PETSC_HAVE_HDF5) */