1:  #include <petsc/private/viewerimpl.h>
  2: #include <petsclayouthdf5.h>    /*I   "petsclayoutdf5.h"   I*/
  3:  #include <petscis.h>

  5: #if defined(PETSC_HAVE_HDF5)

  7: struct _n_HDF5ReadCtx {
  8:   hid_t file, group, dataset, dataspace, plist;
  9:   PetscInt timestep;
 10:   PetscBool complexVal, dim2, horizontal;
 11: };
 12: typedef struct _n_HDF5ReadCtx* HDF5ReadCtx;

 14: static PetscErrorCode PetscViewerHDF5ReadInitialize_Private(PetscViewer viewer, const char name[], HDF5ReadCtx *ctx)
 15: {
 16:   HDF5ReadCtx    h=NULL;

 20:   PetscNew(&h);
 21:   PetscViewerHDF5OpenGroup(viewer, &h->file, &h->group);
 22:   PetscStackCallHDF5Return(h->dataset,H5Dopen2,(h->group, name, H5P_DEFAULT));
 23:   PetscStackCallHDF5Return(h->dataspace,H5Dget_space,(h->dataset));
 24:   PetscViewerHDF5GetTimestep(viewer, &h->timestep);
 25:   PetscViewerHDF5HasAttribute(viewer,name,"complex",&h->complexVal);
 26:   if (h->complexVal) {PetscViewerHDF5ReadAttribute(viewer,name,"complex",PETSC_BOOL,&h->complexVal);}
 27:   /* MATLAB stores column vectors horizontally */
 28:   PetscViewerHDF5HasAttribute(viewer,name,"MATLAB_class",&h->horizontal);
 29:   /* Create property list for collective dataset read */
 30:   PetscStackCallHDF5Return(h->plist,H5Pcreate,(H5P_DATASET_XFER));
 31: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
 32:   PetscStackCallHDF5(H5Pset_dxpl_mpio,(h->plist, H5FD_MPIO_COLLECTIVE));
 33: #endif
 34:   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
 35:   *ctx = h;
 36:   return(0);
 37: }

 39: static PetscErrorCode PetscViewerHDF5ReadFinalize_Private(PetscViewer viewer, HDF5ReadCtx *ctx)
 40: {
 41:   HDF5ReadCtx    h;

 45:   h = *ctx;
 46:   PetscStackCallHDF5(H5Pclose,(h->plist));
 47:   PetscStackCallHDF5(H5Gclose,(h->group));
 48:   PetscStackCallHDF5(H5Sclose,(h->dataspace));
 49:   PetscStackCallHDF5(H5Dclose,(h->dataset));
 50:   PetscFree(*ctx);
 51:   return(0);
 52: }

 54: static PetscErrorCode PetscViewerHDF5ReadSizes_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout *map_)
 55: {
 56:   int            rdim, dim;
 57:   hsize_t        dims[4];
 58:   PetscInt       bsInd, lenInd, bs, len, N;
 59:   PetscLayout    map;

 63:   if (!(*map_)) {
 64:     PetscLayoutCreate(PetscObjectComm((PetscObject)viewer),map_);
 65:   }
 66:   map = *map_;
 67:   /* calculate expected number of dimensions */
 68:   dim = 0;
 69:   if (ctx->timestep >= 0) ++dim;
 70:   ++dim; /* length in blocks */
 71:   if (ctx->complexVal) ++dim;
 72:   /* get actual number of dimensions in dataset */
 73:   PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(ctx->dataspace, dims, NULL));
 74:   /* calculate expected dimension indices */
 75:   lenInd = 0;
 76:   if (ctx->timestep >= 0) ++lenInd;
 77:   bsInd = lenInd + 1;
 78:   ctx->dim2 = PETSC_FALSE;
 79:   if (rdim == dim) {
 80:     bs = 1; /* support vectors stored as 1D array */
 81:   } else if (rdim == dim+1) {
 82:     bs = (PetscInt) dims[bsInd];
 83:     if (bs == 1) ctx->dim2 = PETSC_TRUE; /* vector with blocksize of 1, still stored as 2D array */
 84:   } else {
 85:     SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Number of dimensions %d not %d as expected", rdim, dim);
 86:   }
 87:   len = dims[lenInd];
 88:   if (ctx->horizontal) {
 89:     if (len != 1) SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Cannot have horizontal array with number of rows > 1. In case of MATLAB MAT-file, vectors must be saved as column vectors.");
 90:     len = bs;
 91:     bs = 1;
 92:   }
 93:   N = (PetscInt) len*bs;

 95:   /* Set Vec sizes,blocksize,and type if not already set */
 96:   if (map->bs < 0) {
 97:     PetscLayoutSetBlockSize(map, bs);
 98:   } 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);
 99:   if (map->N < 0) {
100:     PetscLayoutSetSize(map, N);
101:   } 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);
102:   return(0);
103: }

105: static PetscErrorCode PetscViewerHDF5ReadSelectHyperslab_Private(PetscViewer viewer, HDF5ReadCtx ctx, PetscLayout map, hid_t *memspace)
106: {
107:   hsize_t        count[4], offset[4];
108:   int            dim;
109:   PetscInt       bs, n, low;

113:   /* Compute local size and ownership range */
114:   PetscLayoutSetUp(map);
115:   PetscLayoutGetBlockSize(map, &bs);
116:   PetscLayoutGetLocalSize(map, &n);
117:   PetscLayoutGetRange(map, &low, NULL);

119:   /* Each process defines a dataset and reads it from the hyperslab in the file */
120:   dim  = 0;
121:   if (ctx->timestep >= 0) {
122:     count[dim]  = 1;
123:     offset[dim] = ctx->timestep;
124:     ++dim;
125:   }
126:   if (ctx->horizontal) {
127:     count[dim]  = 1;
128:     offset[dim] = 0;
129:     ++dim;
130:   }
131:   {
132:     PetscHDF5IntCast(n/bs, &count[dim]);
133:     PetscHDF5IntCast(low/bs, &offset[dim]);
134:     ++dim;
135:   }
136:   if (bs > 1 || ctx->dim2) {
137:     if (PetscUnlikely(ctx->horizontal)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "cannot have horizontal array with blocksize > 1");
138:     count[dim]  = bs;
139:     offset[dim] = 0;
140:     ++dim;
141:   }
142:   if (ctx->complexVal) {
143:     count[dim]  = 2;
144:     offset[dim] = 0;
145:     ++dim;
146:   }
147:   PetscStackCallHDF5Return(*memspace,H5Screate_simple,(dim, count, NULL));
148:   PetscStackCallHDF5(H5Sselect_hyperslab,(ctx->dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
149:   return(0);
150: }

152: static PetscErrorCode PetscViewerHDF5ReadArray_Private(PetscViewer viewer, HDF5ReadCtx h, hid_t datatype, hid_t memspace, void *arr)
153: {
155:   PetscStackCallHDF5(H5Dread,(h->dataset, datatype, memspace, h->dataspace, h->plist, arr));
156:   return(0);
157: }

159: PetscErrorCode PetscViewerHDF5Load(PetscViewer viewer, const char *name, PetscLayout map, hid_t datatype, void **newarr)
160: {
161:   HDF5ReadCtx     h=NULL;
162:   hid_t           memspace=0;
163:   size_t          unitsize;
164:   void            *arr;
165:   PetscErrorCode  ierr;

168:   PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
169: #if defined(PETSC_USE_COMPLEX)
170:   if (!h->complexVal) {
171:     H5T_class_t clazz = H5Tget_class(datatype);
172:     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.");
173:   }
174: #else
175:   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.");
176: #endif

178:   PetscViewerHDF5ReadSizes_Private(viewer, h, &map);
179:   PetscLayoutSetUp(map);
180:   PetscViewerHDF5ReadSelectHyperslab_Private(viewer, h, map, &memspace);

182:   unitsize = H5Tget_size(datatype);
183:   if (h->complexVal) unitsize *= 2;
184:   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);
185:   PetscMalloc(map->n*unitsize, &arr);

187:   PetscViewerHDF5ReadArray_Private(viewer, h, datatype, memspace, arr);
188:   PetscStackCallHDF5(H5Sclose,(memspace));
189:   PetscViewerHDF5ReadFinalize_Private(viewer, &h);
190:   *newarr = arr;
191:   return(0);
192: }

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

197:   Input Parameters:
198: + viewer - The HDF5 viewer
199: - name   - The dataset name

201:   Output Parameter:
202: + bs     - block size
203: - N      - global size

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

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

211:   Level: advanced

213: .seealso: PetscViewerHDF5Open(), VecLoad(), ISLoad(), VecGetSize(), ISGetSize(), PetscViewerHDF5SetBaseDimension2()
214: @*/
215: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
216: {
217:   HDF5ReadCtx    h=NULL;
218:   PetscLayout    map=NULL;

223:   PetscViewerHDF5ReadInitialize_Private(viewer, name, &h);
224:   PetscViewerHDF5ReadSizes_Private(viewer, h, &map);
225:   PetscViewerHDF5ReadFinalize_Private(viewer, &h);
226:   if (bs) *bs = map->bs;
227:   if (N) *N = map->N;
228:   PetscLayoutDestroy(&map);
229:   return(0);
230: }

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