1: #include <petscis.h>         /*I  "petscis.h"  I*/
  2: #include <petsc/private/isimpl.h>
  3: #include <petscviewerhdf5.h>

  5: #if defined(PETSC_HAVE_HDF5)
  8: /*
  9:      This should handle properly the cases where PetscInt is 32 or 64 and hsize_t is 32 or 64. These means properly casting with
 10:    checks back and forth between the two types of variables.
 11: */
 12: PetscErrorCode ISLoad_HDF5(IS is, PetscViewer viewer)
 13: {
 14:   hid_t           inttype;    /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */
 15:   hid_t           file_id, group, dset_id, filespace, memspace, plist_id;
 16:   int             rdim, dim;
 17:   hsize_t         dims[3], count[3], offset[3];
 18:   PetscInt        n, N, bs, bsInd, lenInd, low, timestep;
 19:   const PetscInt *ind;
 20:   const char     *isname;
 21:   PetscErrorCode  ierr;

 24:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
 25:   PetscViewerHDF5GetTimestep(viewer, &timestep);
 26:   ISGetBlockSize(is, &bs);
 27:   /* Create the dataset with default properties and close filespace */
 28:   PetscObjectGetName((PetscObject) is, &isname);
 29: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
 30:   PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
 31: #else
 32:   PetscStackCallHDF5Return(dset_id,H5Dopen,(group, isname));
 33: #endif
 34:   /* Retrieve the dataspace for the dataset */
 35:   PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
 36:   dim = 0;
 37:   if (timestep >= 0) ++dim;
 38:   ++dim;
 39:   ++dim;
 40:   PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL));
 41:   bsInd = rdim-1;
 42:   lenInd = timestep >= 0 ? 1 : 0;
 43:   if (rdim != dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim);
 44:   else if (bs != (PetscInt) dims[bsInd]) {
 45:     ISSetBlockSize(is, dims[bsInd]);
 46:     if (ierr) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Block size %d specified for IS does not match blocksize in file %d",bs,dims[bsInd]);
 47:     bs = dims[bsInd];
 48:   }

 50:   /* Set Vec sizes,blocksize,and type if not already set */
 51:   ISGetLocalSize(is, &n);
 52:   ISGetSize(is, &N);
 53:   if (n < 0 && N < 0) {PetscLayoutSetSize(is->map, dims[lenInd]*bs);}
 54:   PetscLayoutSetUp(is->map);
 55:   /* If sizes and type already set,check if the vector global size is correct */
 56:   ISGetSize(is, &N);
 57:   if (N/bs != (PetscInt) dims[lenInd]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "IS in file different length (%d) then input vector (%d)", (PetscInt) dims[lenInd], N/bs);

 59:   /* Each process defines a dataset and reads it from the hyperslab in the file */
 60:   ISGetLocalSize(is, &n);
 61:   dim  = 0;
 62:   if (timestep >= 0) {
 63:     count[dim] = 1;
 64:     ++dim;
 65:   }
 66:   PetscHDF5IntCast(n/bs,count + dim);
 67:   ++dim;
 68:   if (bs >= 1) {
 69:     count[dim] = bs;
 70:     ++dim;
 71:   }
 72:   PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));

 74:   /* Select hyperslab in the file */
 75:   PetscLayoutGetRange(is->map, &low, NULL);
 76:   dim  = 0;
 77:   if (timestep >= 0) {
 78:     offset[dim] = timestep;
 79:     ++dim;
 80:   }
 81:   PetscHDF5IntCast(low/bs,offset + dim);
 82:   ++dim;
 83:   if (bs >= 1) {
 84:     offset[dim] = 0;
 85:     ++dim;
 86:   }
 87:   PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));

 89:   /* Create property list for collective dataset read */
 90:   PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER));
 91: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
 92:   PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE));
 93: #endif
 94:   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */

 96: #if defined(PETSC_USE_64BIT_INDICES)
 97:   inttype = H5T_NATIVE_LLONG;
 98: #else
 99:   inttype = H5T_NATIVE_INT;
100: #endif
101:   PetscMalloc1(n,&ind);
102:   PetscStackCallHDF5(H5Dread,(dset_id, inttype, memspace, filespace, plist_id, (void *) ind));
103:   ISGeneralSetIndices(is, n, ind, PETSC_OWN_POINTER);

105:   /* Close/release resources */
106:   if (group != file_id) PetscStackCallHDF5(H5Gclose,(group));
107:   PetscStackCallHDF5(H5Pclose,(plist_id));
108:   PetscStackCallHDF5(H5Sclose,(filespace));
109:   PetscStackCallHDF5(H5Sclose,(memspace));
110:   PetscStackCallHDF5(H5Dclose,(dset_id));
111:   return(0);
112: }
113: #endif

117: PetscErrorCode ISLoad_Default(IS is, PetscViewer viewer)
118: {
119:   PetscBool      isbinary;
120: #if defined(PETSC_HAVE_HDF5)
121:   PetscBool      ishdf5;
122: #endif

126:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
127: #if defined(PETSC_HAVE_HDF5)
128:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
129: #endif
130:   if (isbinary) {
131:     SETERRQ(PetscObjectComm((PetscObject) is), PETSC_ERR_SUP, "This should be implemented");
132: #if defined(PETSC_HAVE_HDF5)
133:   } else if (ishdf5) {
134:     if (!((PetscObject) is)->name) {
135:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Since HDF5 format gives ASCII name for each object in file; must use ISLoad() after setting name of Vec with PetscObjectSetName()");
136:     }
137:     ISLoad_HDF5(is, viewer);
138: #endif
139:   }
140:   return(0);
141: }