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:   hsize_t         rdim, dim;
 17:   hsize_t         dims[3], count[3], offset[3];
 18:   herr_t          status;
 19:   PetscInt        n, N, bs, bsInd, lenInd, low, timestep;
 20:   const PetscInt *ind;
 21:   const char     *isname;
 22:   PetscErrorCode  ierr;

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

 53:   /* Set Vec sizes,blocksize,and type if not already set */
 54:   ISGetLocalSize(is, &n);
 55:   ISGetSize(is, &N);
 56:   if (n < 0 && N < 0) {PetscLayoutSetSize(is->map, dims[lenInd]*bs);}
 57:   PetscLayoutSetUp(is->map);
 58:   /* If sizes and type already set,check if the vector global size is correct */
 59:   ISGetSize(is, &N);
 60:   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);

 62:   /* Each process defines a dataset and reads it from the hyperslab in the file */
 63:   ISGetLocalSize(is, &n);
 64:   dim  = 0;
 65:   if (timestep >= 0) {
 66:     count[dim] = 1;
 67:     ++dim;
 68:   }
 69:   PetscHDF5IntCast(n/bs,count + dim);
 70:   ++dim;
 71:   if (bs >= 1) {
 72:     count[dim] = bs;
 73:     ++dim;
 74:   }
 75:   memspace = H5Screate_simple(dim, count, NULL);
 76:   if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Screate_simple()");

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

 93:   /* Create property list for collective dataset read */
 94:   plist_id = H5Pcreate(H5P_DATASET_XFER);
 95:   if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Pcreate()");
 96: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
 97:   status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status);
 98: #endif
 99:   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */

101: #if defined(PETSC_USE_64BIT_INDICES)
102:   inttype = H5T_NATIVE_LLONG;
103: #else
104:   inttype = H5T_NATIVE_INT;
105: #endif
106:   PetscMalloc1(n,&ind);
107:   status = H5Dread(dset_id, inttype, memspace, filespace, plist_id, (void *) ind);CHKERRQ(status);
108:   ISGeneralSetIndices(is, n, ind, PETSC_OWN_POINTER);

110:   /* Close/release resources */
111:   if (group != file_id) {status = H5Gclose(group);CHKERRQ(status);}
112:   status = H5Pclose(plist_id);CHKERRQ(status);
113:   status = H5Sclose(filespace);CHKERRQ(status);
114:   status = H5Sclose(memspace);CHKERRQ(status);
115:   status = H5Dclose(dset_id);CHKERRQ(status);
116:   return(0);
117: }
118: #endif

122: PetscErrorCode ISLoad_Default(IS is, PetscViewer viewer)
123: {
124:   PetscBool      isbinary;
125: #if defined(PETSC_HAVE_HDF5)
126:   PetscBool      ishdf5;
127: #endif

131:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
132: #if defined(PETSC_HAVE_HDF5)
133:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
134: #endif
135:   if (isbinary) {
136:     SETERRQ(PetscObjectComm((PetscObject) is), PETSC_ERR_SUP, "This should be implemented");
137: #if defined(PETSC_HAVE_HDF5)
138:   } else if (ishdf5) {
139:     if (!((PetscObject) is)->name) {
140:       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()");
141:     }
142:     ISLoad_HDF5(is, viewer);
143: #endif
144:   }
145:   return(0);
146: }