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 /* 7 This should handle properly the cases where PetscInt is 32 or 64 and hsize_t is 32 or 64. These means properly casting with 8 checks back and forth between the two types of variables. 9 */ 10 PetscErrorCode ISLoad_HDF5(IS is, PetscViewer viewer) 11 { 12 hid_t inttype; /* int type (H5T_NATIVE_INT or H5T_NATIVE_LLONG) */ 13 hid_t file_id, group, dset_id, filespace, memspace, plist_id; 14 int rdim, dim; 15 hsize_t dims[3], count[3], offset[3]; 16 PetscInt n, N, bs, bsInd, lenInd, low, timestep; 17 const PetscInt *ind; 18 const char *isname; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 if (!((PetscObject)is)->name) SETERRQ(PetscObjectComm((PetscObject)is), PETSC_ERR_SUP, "Since HDF5 format gives ASCII name for each object in file; must use ISLoad() after setting name of Vec with PetscObjectSetName()"); 23 ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr); 24 ierr = PetscViewerHDF5GetTimestep(viewer, ×tep);CHKERRQ(ierr); 25 ierr = ISGetBlockSize(is, &bs);CHKERRQ(ierr); 26 /* Create the dataset with default properties and close filespace */ 27 ierr = PetscObjectGetName((PetscObject) is, &isname);CHKERRQ(ierr); 28 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 29 PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT)); 30 #else 31 PetscStackCallHDF5Return(dset_id,H5Dopen,(group, isname)); 32 #endif 33 /* Retrieve the dataspace for the dataset */ 34 PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id)); 35 dim = 0; 36 if (timestep >= 0) ++dim; 37 ++dim; 38 ++dim; 39 PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL)); 40 bsInd = rdim-1; 41 lenInd = timestep >= 0 ? 1 : 0; 42 if (rdim != dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim); 43 else if (bs != (PetscInt) dims[bsInd]) { 44 ierr = ISSetBlockSize(is, dims[bsInd]); 45 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]); 46 bs = dims[bsInd]; 47 } 48 49 /* Set Vec sizes,blocksize,and type if not already set */ 50 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 51 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 52 if (n < 0 && N < 0) {ierr = PetscLayoutSetSize(is->map, dims[lenInd]*bs);CHKERRQ(ierr);} 53 ierr = PetscLayoutSetUp(is->map);CHKERRQ(ierr); 54 /* If sizes and type already set,check if the vector global size is correct */ 55 ierr = ISGetSize(is, &N);CHKERRQ(ierr); 56 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); 57 58 /* Each process defines a dataset and reads it from the hyperslab in the file */ 59 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr); 60 dim = 0; 61 if (timestep >= 0) { 62 count[dim] = 1; 63 ++dim; 64 } 65 ierr = PetscHDF5IntCast(n/bs,count + dim);CHKERRQ(ierr); 66 ++dim; 67 if (bs >= 1) { 68 count[dim] = bs; 69 ++dim; 70 } 71 PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL)); 72 73 /* Select hyperslab in the file */ 74 ierr = PetscLayoutGetRange(is->map, &low, NULL);CHKERRQ(ierr); 75 dim = 0; 76 if (timestep >= 0) { 77 offset[dim] = timestep; 78 ++dim; 79 } 80 ierr = PetscHDF5IntCast(low/bs,offset + dim);CHKERRQ(ierr); 81 ++dim; 82 if (bs >= 1) { 83 offset[dim] = 0; 84 ++dim; 85 } 86 PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL)); 87 88 /* Create property list for collective dataset read */ 89 PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER)); 90 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 91 PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE)); 92 #endif 93 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 94 95 #if defined(PETSC_USE_64BIT_INDICES) 96 inttype = H5T_NATIVE_LLONG; 97 #else 98 inttype = H5T_NATIVE_INT; 99 #endif 100 ierr = PetscMalloc1(n,&ind);CHKERRQ(ierr); 101 PetscStackCallHDF5(H5Dread,(dset_id, inttype, memspace, filespace, plist_id, (void *) ind)); 102 ierr = ISGeneralSetIndices(is, n, ind, PETSC_OWN_POINTER);CHKERRQ(ierr); 103 104 /* Close/release resources */ 105 if (group != file_id) PetscStackCallHDF5(H5Gclose,(group)); 106 PetscStackCallHDF5(H5Pclose,(plist_id)); 107 PetscStackCallHDF5(H5Sclose,(filespace)); 108 PetscStackCallHDF5(H5Sclose,(memspace)); 109 PetscStackCallHDF5(H5Dclose,(dset_id)); 110 PetscFunctionReturn(0); 111 } 112 #endif 113 114 PetscErrorCode ISLoad_Binary(IS is, PetscViewer viewer) 115 { 116 PetscErrorCode ierr; 117 PetscBool isgeneral,skipHeader,useMPIIO; 118 int fd; 119 PetscInt tr[2],N,ln,*idx; 120 MPI_Request request; 121 MPI_Status status; 122 MPI_Comm comm; 123 PetscMPIInt rank,size,tag; 124 125 PetscFunctionBegin; 126 ierr = PetscObjectGetComm((PetscObject)is,&comm);CHKERRQ(ierr); 127 ierr = PetscObjectTypeCompare((PetscObject)is,ISGENERAL,&isgeneral);CHKERRQ(ierr); 128 if (!isgeneral) SETERRQ(comm,PETSC_ERR_ARG_INCOMP,"IS must be of type ISGENERAL to load into it"); 129 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 130 if (skipHeader) SETERRQ(comm,PETSC_ERR_USER, "Currently no support for binary files without headers"); 131 /* force binary viewer to load .info file if it has not yet done so */ 132 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 133 134 ierr = PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);CHKERRQ(ierr); 135 if (tr[0] != IS_FILE_CLASSID) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Not an IS next in file"); 136 137 /* Has IS already had its layout defined */ 138 /* ierr = ISGetLayout(is,&map);CHKERRQ(ierr); */ 139 ierr = PetscLayoutGetSize(is->map,&N);CHKERRQ(ierr); 140 if (N > -1 && N != tr[1]) SETERRQ2(comm,PETSC_ERR_ARG_SIZ,"Size of IS in file %D does not match size of IS provided",tr[1],N); 141 if (N == -1) { 142 N = tr[1]; 143 ierr = PetscLayoutSetSize(is->map,N);CHKERRQ(ierr); 144 ierr = PetscLayoutSetUp(is->map);CHKERRQ(ierr); 145 } 146 ierr = PetscLayoutGetLocalSize(is->map,&ln);CHKERRQ(ierr); 147 ierr = PetscMalloc1(ln,&idx);CHKERRQ(ierr); 148 149 ierr = PetscViewerBinaryGetUseMPIIO(viewer,&useMPIIO);CHKERRQ(ierr); 150 #if defined(PETSC_HAVE_MPIIO) 151 if (useMPIIO) { 152 MPI_File mfdes; 153 MPI_Offset off; 154 PetscMPIInt lsize; 155 PetscInt rstart; 156 157 ierr = PetscMPIIntCast(ln,&lsize);CHKERRQ(ierr); 158 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 159 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 160 ierr = PetscLayoutGetRange(is->map,&rstart,NULL);CHKERRQ(ierr); 161 off += rstart*(MPI_Offset)sizeof(PetscInt); 162 ierr = MPI_File_set_view(mfdes,off,MPIU_INT,MPIU_INT,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr); 163 ierr = MPIU_File_read_all(mfdes,idx,lsize,MPIU_INT,MPI_STATUS_IGNORE);CHKERRQ(ierr); 164 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,N*(MPI_Offset)sizeof(PetscInt));CHKERRQ(ierr); 165 ierr = ISGeneralSetIndices(is,ln,idx,PETSC_OWN_POINTER);CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 #endif 169 170 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 171 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 172 ierr = PetscObjectGetNewTag((PetscObject)viewer,&tag);CHKERRQ(ierr); 173 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 174 175 if (!rank) { 176 ierr = PetscBinaryRead(fd,idx,ln,PETSC_INT);CHKERRQ(ierr); 177 178 if (size > 1) { 179 PetscInt *range,n,i,*idxwork; 180 181 /* read in other chuncks and send to other processors */ 182 /* determine maximum chunck owned by other */ 183 range = is->map->range; 184 n = 1; 185 for (i=1; i<size; i++) n = PetscMax(n,range[i+1] - range[i]); 186 187 ierr = PetscMalloc1(n,&idxwork);CHKERRQ(ierr); 188 for (i=1; i<size; i++) { 189 n = range[i+1] - range[i]; 190 ierr = PetscBinaryRead(fd,idxwork,n,PETSC_INT);CHKERRQ(ierr); 191 ierr = MPI_Isend(idxwork,n,MPIU_INT,i,tag,comm,&request);CHKERRQ(ierr); 192 ierr = MPI_Wait(&request,&status);CHKERRQ(ierr); 193 } 194 ierr = PetscFree(idxwork);CHKERRQ(ierr); 195 } 196 } else { 197 ierr = MPI_Recv(idx,ln,MPIU_INT,0,tag,comm,&status);CHKERRQ(ierr); 198 } 199 ierr = ISGeneralSetIndices(is,ln,idx,PETSC_OWN_POINTER);CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202 203 PetscErrorCode ISLoad_Default(IS is, PetscViewer viewer) 204 { 205 PetscBool isbinary; 206 #if defined(PETSC_HAVE_HDF5) 207 PetscBool ishdf5; 208 #endif 209 PetscErrorCode ierr; 210 211 PetscFunctionBegin; 212 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 213 #if defined(PETSC_HAVE_HDF5) 214 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 215 #endif 216 if (isbinary) { 217 ierr = ISLoad_Binary(is, viewer);CHKERRQ(ierr); 218 #if defined(PETSC_HAVE_HDF5) 219 } else if (ishdf5) { 220 ierr = ISLoad_HDF5(is, viewer);CHKERRQ(ierr); 221 #endif 222 } 223 PetscFunctionReturn(0); 224 } 225