Actual source code: hdf5io.c

      1: #include <petsc/private/viewerhdf5impl.h>
      2: #include <petsclayouthdf5.h>
      3: #include <petscis.h>

      5: struct _n_HDF5ReadCtx {
      6:   const char *name;
      7:   hid_t       file, group, dataset, dataspace;
      8:   int         lenInd, bsInd, complexInd, rdim;
      9:   hsize_t    *dims;
     10:   PetscBool   complexVal, dim2;

     12:   // Needed for compression
     13:   PetscInt  runs;
     14:   PetscInt *cind;
     15: };
     16: typedef struct _n_HDF5ReadCtx *HDF5ReadCtx;

     18: PetscErrorCode PetscViewerHDF5CheckTimestepping_Internal(PetscViewer viewer, const char name[])
     19: {
     20:   PetscViewer_HDF5 *hdf5         = (PetscViewer_HDF5 *)viewer->data;
     21:   PetscBool         timestepping = PETSC_FALSE;

     23:   PetscFunctionBegin;
     24:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, name, "timestepping", PETSC_BOOL, &hdf5->defTimestepping, &timestepping));
     25:   if (timestepping != hdf5->timestepping) {
     26:     const char *group;

     28:     PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &group));
     29:     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]);
     30:   }
     31:   PetscFunctionReturn(PETSC_SUCCESS);
     32: }

     34: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
     35: {
     36:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
     37:   HDF5ReadCtx       h    = NULL;

     39:   PetscFunctionBegin;
     40:   PetscCall(PetscViewerHDF5CheckTimestepping_Internal(viewer, name));
     41:   PetscCall(PetscNew(&h));
     42:   h->name = name;
     43:   PetscCall(PetscViewerHDF5OpenGroup(viewer, NULL, &h->file, &h->group));
     44:   PetscCallHDF5Return(h->dataset, H5Dopen2, (h->group, name, H5P_DEFAULT));
     45:   PetscCallHDF5Return(h->dataspace, H5Dget_space, (h->dataset));
     46:   PetscCall(PetscViewerHDF5ReadAttribute(viewer, name, "complex", PETSC_BOOL, &h->complexVal, &h->complexVal));
     47:   if (!hdf5->horizontal) {
     48:     /* MATLAB stores column vectors horizontally */
     49:     PetscCall(PetscViewerHDF5HasAttribute(viewer, name, "MATLAB_class", &hdf5->horizontal));
     50:   }
     51:   h->runs = 0;
     52:   h->cind = NULL;
     53:   *ctx    = h;
     54:   PetscFunctionReturn(PETSC_SUCCESS);
     55: }

     57: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
     58: {
     59:   HDF5ReadCtx h;

     61:   PetscFunctionBegin;
     62:   h = *ctx;
     63:   PetscCallHDF5(H5Gclose, (h->group));
     64:   PetscCallHDF5(H5Sclose, (h->dataspace));
     65:   PetscCallHDF5(H5Dclose, (h->dataset));
     66:   PetscCall(PetscFree((*ctx)->dims));
     67:   PetscCall(PetscFree((*ctx)->cind));
     68:   PetscCall(PetscFree(*ctx));
     69:   PetscFunctionReturn(PETSC_SUCCESS);
     70: }

     72: // Need forward declaration because we have a cyclic call chain
     73: static PetscErrorCode PetscViewerHDF5Load_Internal(PetscViewer, const char[], PetscBool, PetscLayout, hid_t, void **);

     75: static PetscErrorCode PetscViewerHDF5ReadSizes_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscBool uncompress, PetscBool setup, PetscLayout *map_)
     76: {
     77:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
     78:   PetscInt          bs, N;
     79:   PetscLayout       map;
     80:   PetscBool         compressed;

     82:   PetscFunctionBegin;
     83:   if (!*map_) PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)viewer), map_));
     84:   map = *map_;

     86:   PetscCall(PetscViewerHDF5HasAttribute(viewer, ctx->name, "compressed", &compressed));
     87:   if (compressed && uncompress) {
     88:     hid_t           inttype;
     89:     PetscLayout     cmap;
     90:     PetscInt       *lcind, N = 0;
     91:     PetscMPIInt    *counts, *displs, size, n;
     92:     const PetscInt *range;
     93:     MPI_Comm        comm;

     95: #if defined(PETSC_USE_64BIT_INDICES)
     96:     inttype = H5T_NATIVE_LLONG;
     97: #else
     98:     inttype = H5T_NATIVE_INT;
     99: #endif
    100:     PetscCall(PetscObjectGetComm((PetscObject)viewer, &comm));
    101:     PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)viewer), &cmap));
    102:     cmap->bs = 3;
    103:     PetscCall(PetscViewerHDF5Load_Internal(viewer, ctx->name, PETSC_FALSE, cmap, inttype, (void **)&lcind));
    104:     PetscCheck(!(cmap->n % 3), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Compressed IS must have an even number of entries, not %" PetscInt_FMT, cmap->n);
    105:     for (PetscInt i = 0; i < cmap->n / 3; ++i) N += lcind[i * 3 + 0];
    106:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &N, 1, MPIU_INT, MPIU_SUM, comm));
    107:     ctx->runs = cmap->N / 3;
    108:     PetscCall(PetscMalloc1(cmap->N, &ctx->cind));
    109:     PetscCallMPI(MPI_Comm_size(comm, &size));
    110:     PetscCall(PetscLayoutGetRanges(cmap, &range));
    111:     PetscCall(PetscMalloc2(size, &counts, size, &displs));
    112:     for (PetscInt r = 0; r < size; ++r) {
    113:       PetscCall(PetscMPIIntCast(range[r + 1] - range[r], &counts[r]));
    114:       PetscCall(PetscMPIIntCast(range[r], &displs[r]));
    115:     }
    116:     PetscCall(PetscMPIIntCast(cmap->n, &n));
    117:     PetscCallMPI(MPI_Allgatherv(lcind, n, MPIU_INT, ctx->cind, counts, displs, MPIU_INT, comm));
    118:     PetscCall(PetscFree2(counts, displs));
    119:     PetscCall(PetscFree(lcind));
    120:     PetscCall(PetscLayoutDestroy(&cmap));

    122:     ctx->dim2   = PETSC_FALSE;
    123:     ctx->rdim   = 1;
    124:     ctx->lenInd = 0;
    125:     PetscCall(PetscMalloc1(ctx->rdim, &ctx->dims));
    126:     ctx->dims[0] = N;
    127:     bs           = 1;
    128:     goto layout;
    129:   }

    131:   /* Get actual number of dimensions in dataset */
    132:   PetscCallHDF5Return(ctx->rdim, H5Sget_simple_extent_dims, (ctx->dataspace, NULL, NULL));
    133:   PetscCall(PetscMalloc1(ctx->rdim, &ctx->dims));
    134:   PetscCallHDF5Return(ctx->rdim, H5Sget_simple_extent_dims, (ctx->dataspace, ctx->dims, NULL));

    136:   /*
    137:      Dimensions are in this order:
    138:      [0]        timesteps (optional)
    139:      [lenInd]   entries (numbers or blocks)
    140:      ...
    141:      [bsInd]    entries of blocks (optional)
    142:      [bsInd+1]  real & imaginary part (optional)
    143:       = rdim-1
    144:    */

    146:   /* Get entries dimension index */
    147:   ctx->lenInd = 0;
    148:   if (hdf5->timestepping) ++ctx->lenInd;

    150:   /* Get block dimension index */
    151:   if (ctx->complexVal) {
    152:     ctx->bsInd      = ctx->rdim - 2;
    153:     ctx->complexInd = ctx->rdim - 1;
    154:   } else {
    155:     ctx->bsInd      = ctx->rdim - 1;
    156:     ctx->complexInd = -1;
    157:   }
    158:   PetscCheck(ctx->lenInd <= ctx->bsInd, PetscObjectComm((PetscObject)viewer), PETSC_ERR_PLIB, "Calculated block dimension index = %d < %d = length dimension index.", ctx->bsInd, ctx->lenInd);
    159:   PetscCheck(ctx->bsInd <= ctx->rdim - 1, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Calculated block dimension index = %d > %d = total number of dimensions - 1.", ctx->bsInd, ctx->rdim - 1);
    160:   PetscCheck(!ctx->complexVal || ctx->dims[ctx->complexInd] == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Complex numbers must have exactly 2 parts (%" PRIuHSIZE ")", ctx->dims[ctx->complexInd]);

    162:   if (hdf5->horizontal) {
    163:     /* support horizontal 1D arrays (MATLAB vectors) - swap meaning of blocks and entries */
    164:     int t       = ctx->lenInd;
    165:     ctx->lenInd = ctx->bsInd;
    166:     ctx->bsInd  = t;
    167:   }

    169:   /* Get block size */
    170:   ctx->dim2 = PETSC_FALSE;
    171:   if (ctx->lenInd == ctx->bsInd) {
    172:     bs = 1; /* support vectors stored as 1D array */
    173:   } else {
    174:     bs = (PetscInt)ctx->dims[ctx->bsInd];
    175:     if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
    176:   }

    178: layout:
    179:   /* Get global size */
    180:   PetscCall(PetscIntCast(bs * ctx->dims[ctx->lenInd], &N));

    182:   /* Set global size, blocksize and type if not yet set */
    183:   PetscCall(PetscLayoutSetBlockSize(map, bs));
    184:   if (map->N < 0) PetscCall(PetscLayoutSetSize(map, N));
    185:   else PetscCheck(map->N == N, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Global size of array %s in file is %" PetscInt_FMT ", not %" PetscInt_FMT " as expected", ctx->name, N, map->N);
    186:   if (setup) PetscCall(PetscLayoutSetUp(map));
    187:   PetscFunctionReturn(PETSC_SUCCESS);
    188: }

    190: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
    191: {
    192:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;
    193:   hsize_t          *count, *offset;
    194:   PetscInt          bs, n, low;
    195:   int               i;

    197:   PetscFunctionBegin;
    198:   /* Compute local size and ownership range */
    199:   PetscCall(PetscLayoutSetUp(map));
    200:   PetscCall(PetscLayoutGetBlockSize(map, &bs));
    201:   PetscCall(PetscLayoutGetLocalSize(map, &n));
    202:   PetscCall(PetscLayoutGetRange(map, &low, NULL));

    204:   /* Each process defines a dataset and reads it from the hyperslab in the file */
    205:   PetscCall(PetscMalloc2(ctx->rdim, &count, ctx->rdim, &offset));
    206:   for (i = 0; i < ctx->rdim; i++) {
    207:     /* By default, select all entries with no offset */
    208:     offset[i] = 0;
    209:     count[i]  = ctx->dims[i];
    210:   }
    211:   if (hdf5->timestepping) {
    212:     count[0]  = 1;
    213:     offset[0] = hdf5->timestep;
    214:   }
    215:   {
    216:     PetscCall(PetscHDF5IntCast(n / bs, &count[ctx->lenInd]));
    217:     PetscCall(PetscHDF5IntCast(low / bs, &offset[ctx->lenInd]));
    218:   }
    219:   PetscCallHDF5Return(*memspace, H5Screate_simple, (ctx->rdim, count, NULL));
    220:   PetscCallHDF5(H5Sselect_hyperslab, (ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
    221:   PetscCall(PetscFree2(count, offset));
    222:   PetscFunctionReturn(PETSC_SUCCESS);
    223: }

    225: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
    226: {
    227:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5 *)viewer->data;

    229:   PetscFunctionBegin;
    230:   PetscCallHDF5(H5Dread, (h->dataset, datatype, memspace, h->dataspace, hdf5->dxpl_id, arr));
    231:   PetscFunctionReturn(PETSC_SUCCESS);
    232: }

    234: static PetscErrorCode PetscViewerHDF5Load_Internal(PetscViewer viewer, const char name[], PetscBool uncompress, PetscLayout map, hid_t datatype, void **newarr)
    235: {
    236:   PetscBool   has;
    237:   const char *group;
    238:   HDF5ReadCtx h        = NULL;
    239:   hid_t       memspace = 0;
    240:   size_t      unitsize;
    241:   void       *arr;

    243:   PetscFunctionBegin;
    244:   PetscCall(PetscViewerHDF5GetGroup(viewer, NULL, &group));
    245:   PetscCall(PetscViewerHDF5HasDataset(viewer, name, &has));
    246:   PetscCheck(has, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Object (dataset) \"%s\" not stored in group %s", name, group);
    247:   PetscCall(PetscViewerHDF5ReadInitialize_Private(viewer, name, &h));
    248: #if defined(PETSC_USE_COMPLEX)
    249:   if (!h->complexVal) {
    250:     H5T_class_t clazz = H5Tget_class(datatype);
    251:     PetscCheck(clazz != H5T_FLOAT, 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);
    252:   }
    253: #else
    254:   PetscCheck(!h->complexVal, 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, name);
    255: #endif

    257:   PetscCall(PetscViewerHDF5ReadSizes_Private(viewer, h, uncompress, PETSC_TRUE, &map));
    258:   PetscCall(PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace));

    260:   if (h->runs && uncompress) {
    261:     PetscInt *ind;

    263:     PetscCall(PetscInfo(viewer, "Read compressed object with name %s of size %" PetscInt_FMT ":%" PetscInt_FMT "\n", name, map->n, map->N));
    264:     // Each process stores the whole compression, so skip any leading parts
    265:     PetscCall(PetscMalloc1(map->n, &ind));
    266:     for (PetscInt i = 0, off = 0; i < h->runs; ++i) {
    267:       for (PetscInt j = 0, inc = 0; j < h->cind[i * 3 + 0]; ++j, ++off, inc += h->cind[i * 3 + 1]) {
    268:         if (off >= map->rend) {
    269:           i = h->runs;
    270:           break;
    271:         }
    272:         if (off >= map->rstart) ind[off - map->rstart] = h->cind[i * 3 + 2] + inc;
    273:       }
    274:     }
    275:     *newarr = ind;
    276:     goto cleanup;
    277:   }

    279:   unitsize = H5Tget_size(datatype);
    280:   if (h->complexVal) unitsize *= 2;
    281:   /* unitsize is size_t i.e. always unsigned, so the negative check is pointless? */
    282:   PetscCheck(unitsize > 0 && unitsize <= PetscMax(sizeof(PetscInt), sizeof(PetscScalar)), PETSC_COMM_SELF, PETSC_ERR_LIB, "Sanity check failed: HDF5 function H5Tget_size(datatype) returned suspicious value %zu", unitsize);
    283:   PetscCall(PetscMalloc(map->n * unitsize, &arr));

    285:   PetscCall(PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr));
    286:   *newarr = arr;

    288: cleanup:
    289:   PetscCallHDF5(H5Sclose, (memspace));
    290:   PetscCall(PetscViewerHDF5ReadFinalize_Private(viewer, &h));
    291:   PetscCall(PetscFree(group));
    292:   PetscFunctionReturn(PETSC_SUCCESS);
    293: }

    295: /*@C
    296:   PetscViewerHDF5Load - Read a raw array from the `PETSCVIEWERHDF5` dataset in parallel

    298:   Collective; No Fortran Support

    300:   Input Parameters:
    301: + viewer   - The `PETSCVIEWERHDF5` viewer
    302: . name     - The dataset name
    303: - datatype - The HDF5 datatype of the items in the dataset

    305:   Input/Output Parameter:
    306: . map - The layout which specifies array partitioning, on output the
    307:              set up layout (with global size and blocksize according to dataset)

    309:   Output Parameter:
    310: . newarr - The partitioned array, a memory image of the given dataset

    312:   Level: developer

    314:   Notes:
    315:   This is intended mainly for internal use; users should use higher level routines such as `ISLoad()`, `VecLoad()`, `DMLoad()`.

    317:   The array is partitioned according to the given `PetscLayout` which is converted to an HDF5 hyperslab.

    319:   This name is relative to the current group returned by `PetscViewerHDF5OpenGroup()`.

    321: .seealso: `PetscViewer`, `PETSCVIEWERHDF5`, `PetscViewerHDF5Open()`, `PetscViewerHDF5PushGroup()`, `PetscViewerHDF5OpenGroup()`, `PetscViewerHDF5ReadSizes()`,
    322:           `VecLoad()`, `ISLoad()`, `PetscLayout`
    323: @*/
    324: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char name[], PetscLayout map, hid_t datatype, void **newarr)
    325: {
    326:   PetscFunctionBegin;
    327:   PetscCall(PetscViewerHDF5Load_Internal(viewer, name, PETSC_TRUE, map, datatype, newarr));
    328:   PetscFunctionReturn(PETSC_SUCCESS);
    329: }

    331: /*@
    332:   PetscViewerHDF5ReadSizes - Read block size and global size of a `Vec` or `IS` stored in an HDF5 file.

    334:   Input Parameters:
    335: + viewer - The `PETSCVIEWERHDF5` viewer
    336: - name   - The dataset name

    338:   Output Parameters:
    339: + bs - block size
    340: - N  - global size

    342:   Level: advanced

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

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

    350: .seealso: `PetscViewer`, `PETSCVIEWERHDF5`, `PetscViewerHDF5Open()`, `VecLoad()`, `ISLoad()`, `VecGetSize()`, `ISGetSize()`, `PetscViewerHDF5SetBaseDimension2()`
    351: @*/
    352: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
    353: {
    354:   HDF5ReadCtx h   = NULL;
    355:   PetscLayout map = NULL;

    357:   PetscFunctionBegin;
    359:   PetscCall(PetscViewerHDF5ReadInitialize_Private(viewer, name, &h));
    360:   PetscCall(PetscViewerHDF5ReadSizes_Private(viewer, h, PETSC_TRUE, PETSC_FALSE, &map));
    361:   PetscCall(PetscViewerHDF5ReadFinalize_Private(viewer, &h));
    362:   if (bs) *bs = map->bs;
    363:   if (N) *N = map->N;
    364:   PetscCall(PetscLayoutDestroy(&map));
    365:   PetscFunctionReturn(PETSC_SUCCESS);
    366: }