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 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; 22 23 PetscFunctionBegin; 24 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 25 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 26 ierr = ISGetBlockSize(is, &bs);CHKERRQ(ierr); 27 /* Create the dataset with default properties and close filespace */ 28 ierr = PetscObjectGetName((PetscObject) is, &isname);CHKERRQ(ierr); 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 ierr = 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 } 49 50 /* Set Vec sizes,blocksize,and type if not already set */ 51 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 52 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 53 if (n < 0 && N < 0) {ierr = PetscLayoutSetSize(is->map, dims[lenInd]*bs);CHKERRQ(ierr);} 54 ierr = PetscLayoutSetUp(is->map);CHKERRQ(ierr); 55 /* If sizes and type already set,check if the vector global size is correct */ 56 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 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); 58 59 /* Each process defines a dataset and reads it from the hyperslab in the file */ 60 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 61 dim = 0; 62 if (timestep >= 0) { 63 count[dim] = 1; 64 ++dim; 65 } 66 ierr = PetscHDF5IntCast(n/bs,count + dim);CHKERRQ(ierr); 67 ++dim; 68 if (bs >= 1) { 69 count[dim] = bs; 70 ++dim; 71 } 72 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 73 74 /* Select hyperslab in the file */ 75 ierr = PetscLayoutGetRange(is->map, &low, NULL);CHKERRQ(ierr); 76 dim = 0; 77 if (timestep >= 0) { 78 offset[dim] = timestep; 79 ++dim; 80 } 81 ierr = PetscHDF5IntCast(low/bs,offset + dim);CHKERRQ(ierr); 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)); 88 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) */ 95 96 #if defined(PETSC_USE_64BIT_INDICES) 97 inttype = H5T_NATIVE_LLONG; 98 #else 99 inttype = H5T_NATIVE_INT; 100 #endif 101 ierr = PetscMalloc1(n,&ind);CHKERRQ(ierr); 102 PetscStackCallHDF5(H5Dread,(dset_id, inttype, memspace, filespace, plist_id, (void *) ind)); 103 ierr = ISGeneralSetIndices(is, n, ind, PETSC_OWN_POINTER);CHKERRQ(ierr); 104 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 PetscFunctionReturn(0); 112 } 113 #endif 114 115 #undef __FUNCT__ 116 #define __FUNCT__ "ISLoad_Default" 117 PetscErrorCode ISLoad_Default(IS is, PetscViewer viewer) 118 { 119 PetscBool isbinary; 120 #if defined(PETSC_HAVE_HDF5) 121 PetscBool ishdf5; 122 #endif 123 PetscErrorCode ierr; 124 125 PetscFunctionBegin; 126 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 127 #if defined(PETSC_HAVE_HDF5) 128 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 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 ierr = ISLoad_HDF5(is, viewer);CHKERRQ(ierr); 138 #endif 139 } 140 PetscFunctionReturn(0); 141 } 142