1 #include <petscis.h> /*I "petscis.h" I*/ 2 #include <petsc-private/isimpl.h> 3 #include <petscviewerhdf5.h> 4 5 #if defined(PETSC_HAVE_HDF5) 6 #undef __FUNCT__ 7 #define __FUNCT__ "ISLoad_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 = 1, bsInd, lenInd, low, timestep; 20 const PetscInt *ind; 21 const char *isname; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 26 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 27 ierr = ISGetBlockSize(is, &bs);CHKERRQ(ierr); 28 /* Create the dataset with default properties and close filespace */ 29 ierr = PetscObjectGetName((PetscObject) is, &isname);CHKERRQ(ierr); 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 if (bs >= 1) ++dim; 43 rdim = H5Sget_simple_extent_dims(filespace, dims, NULL); 44 bsInd = rdim-1; 45 lenInd = timestep >= 0 ? 1 : 0; 46 if (rdim != dim) { 47 if (rdim == dim+1 && bs == -1) bs = dims[bsInd]; 48 else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim); 49 } else if (bs >= 1 && bs != (PetscInt) dims[bsInd]) { 50 ierr = ISSetBlockSize(is, dims[bsInd]);CHKERRQ(ierr); 51 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]); 52 bs = dims[bsInd]; 53 } 54 55 /* Set Vec sizes,blocksize,and type if not already set */ 56 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 57 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 58 if (n < 0 && N < 0) {ierr = PetscLayoutSetSize(is->map, dims[lenInd]*bs);CHKERRQ(ierr);} 59 ierr = PetscLayoutSetUp(is->map);CHKERRQ(ierr); 60 /* If sizes and type already set,check if the vector global size is correct */ 61 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 62 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); 63 64 /* Each process defines a dataset and reads it from the hyperslab in the file */ 65 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 66 dim = 0; 67 if (timestep >= 0) { 68 count[dim] = 1; 69 ++dim; 70 } 71 ierr = PetscHDF5IntCast(n/bs,count + dim);CHKERRQ(ierr); 72 ++dim; 73 if (bs >= 1) { 74 count[dim] = bs; 75 ++dim; 76 } 77 memspace = H5Screate_simple(dim, count, NULL); 78 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Screate_simple()"); 79 80 /* Select hyperslab in the file */ 81 ierr = PetscLayoutGetRange(is->map, &low, NULL);CHKERRQ(ierr); 82 dim = 0; 83 if (timestep >= 0) { 84 offset[dim] = timestep; 85 ++dim; 86 } 87 ierr = PetscHDF5IntCast(low/bs,offset + dim);CHKERRQ(ierr); 88 ++dim; 89 if (bs >= 1) { 90 offset[dim] = 0; 91 ++dim; 92 } 93 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 94 95 /* Create property list for collective dataset read */ 96 plist_id = H5Pcreate(H5P_DATASET_XFER); 97 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Pcreate()"); 98 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 99 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 100 #endif 101 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 102 103 #if defined(PETSC_USE_64BIT_INDICES) 104 inttype = H5T_NATIVE_LLONG; 105 #else 106 inttype = H5T_NATIVE_INT; 107 #endif 108 ierr = PetscMalloc1(n,&ind);CHKERRQ(ierr); 109 status = H5Dread(dset_id, inttype, memspace, filespace, plist_id, (void *) ind);CHKERRQ(status); 110 ierr = ISGeneralSetIndices(is, n, ind, PETSC_OWN_POINTER);CHKERRQ(ierr); 111 112 /* Close/release resources */ 113 if (group != file_id) {status = H5Gclose(group);CHKERRQ(status);} 114 status = H5Pclose(plist_id);CHKERRQ(status); 115 status = H5Sclose(filespace);CHKERRQ(status); 116 status = H5Sclose(memspace);CHKERRQ(status); 117 status = H5Dclose(dset_id);CHKERRQ(status); 118 PetscFunctionReturn(0); 119 } 120 #endif 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "ISLoad_Default" 124 PetscErrorCode ISLoad_Default(IS is, PetscViewer viewer) 125 { 126 PetscBool isbinary; 127 #if defined(PETSC_HAVE_HDF5) 128 PetscBool ishdf5; 129 #endif 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 134 #if defined(PETSC_HAVE_HDF5) 135 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 136 #endif 137 if (isbinary) { 138 SETERRQ(PetscObjectComm((PetscObject) is), PETSC_ERR_SUP, "This should be implemented"); 139 #if defined(PETSC_HAVE_HDF5) 140 } else if (ishdf5) { 141 if (!((PetscObject) is)->name) { 142 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()"); 143 } 144 ierr = ISLoad_HDF5(is, viewer);CHKERRQ(ierr); 145 #endif 146 } 147 PetscFunctionReturn(0); 148 } 149