Actual source code: hdf5io.c

petsc-master 2020-11-28
Report Typos and Errors
  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) */