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, 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 ++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 ierr = 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 } 52 53 /* Set Vec sizes,blocksize,and type if not already set */ 54 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 55 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 56 if (n < 0 && N < 0) {ierr = PetscLayoutSetSize(is->map, dims[lenInd]*bs);CHKERRQ(ierr);} 57 ierr = PetscLayoutSetUp(is->map);CHKERRQ(ierr); 58 /* If sizes and type already set,check if the vector global size is correct */ 59 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 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); 61 62 /* Each process defines a dataset and reads it from the hyperslab in the file */ 63 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 64 dim = 0; 65 if (timestep >= 0) { 66 count[dim] = 1; 67 ++dim; 68 } 69 ierr = PetscHDF5IntCast(n/bs,count + dim);CHKERRQ(ierr); 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()"); 77 78 /* Select hyperslab in the file */ 79 ierr = PetscLayoutGetRange(is->map, &low, NULL);CHKERRQ(ierr); 80 dim = 0; 81 if (timestep >= 0) { 82 offset[dim] = timestep; 83 ++dim; 84 } 85 ierr = PetscHDF5IntCast(low/bs,offset + dim);CHKERRQ(ierr); 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); 92 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) */ 100 101 #if defined(PETSC_USE_64BIT_INDICES) 102 inttype = H5T_NATIVE_LLONG; 103 #else 104 inttype = H5T_NATIVE_INT; 105 #endif 106 ierr = PetscMalloc1(n,&ind);CHKERRQ(ierr); 107 status = H5Dread(dset_id, inttype, memspace, filespace, plist_id, (void *) ind);CHKERRQ(status); 108 ierr = ISGeneralSetIndices(is, n, ind, PETSC_OWN_POINTER);CHKERRQ(ierr); 109 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 PetscFunctionReturn(0); 117 } 118 #endif 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "ISLoad_Default" 122 PetscErrorCode ISLoad_Default(IS is, PetscViewer viewer) 123 { 124 PetscBool isbinary; 125 #if defined(PETSC_HAVE_HDF5) 126 PetscBool ishdf5; 127 #endif 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 132 #if defined(PETSC_HAVE_HDF5) 133 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 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 ierr = ISLoad_HDF5(is, viewer);CHKERRQ(ierr); 143 #endif 144 } 145 PetscFunctionReturn(0); 146 } 147