Actual source code: isio.c
petsc-3.7.7 2017-09-25
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: 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;
24: 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()");
25: PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
26: PetscViewerHDF5GetTimestep(viewer, ×tep);
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: PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, isname, H5P_DEFAULT));
32: #else
33: PetscStackCallHDF5Return(dset_id,H5Dopen,(group, isname));
34: #endif
35: /* Retrieve the dataspace for the dataset */
36: PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
37: dim = 0;
38: if (timestep >= 0) ++dim;
39: ++dim;
40: ++dim;
41: PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL));
42: bsInd = rdim-1;
43: lenInd = timestep >= 0 ? 1 : 0;
44: if (rdim != dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim);
45: else if (bs != (PetscInt) dims[bsInd]) {
46: ISSetBlockSize(is, dims[bsInd]);
47: 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]);
48: bs = dims[bsInd];
49: }
51: /* Set Vec sizes,blocksize,and type if not already set */
52: ISGetLocalSize(is, &n);
53: ISGetSize(is, &N);
54: if (n < 0 && N < 0) {PetscLayoutSetSize(is->map, dims[lenInd]*bs);}
55: PetscLayoutSetUp(is->map);
56: /* If sizes and type already set,check if the vector global size is correct */
57: ISGetSize(is, &N);
58: 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);
60: /* Each process defines a dataset and reads it from the hyperslab in the file */
61: ISGetLocalSize(is, &n);
62: dim = 0;
63: if (timestep >= 0) {
64: count[dim] = 1;
65: ++dim;
66: }
67: PetscHDF5IntCast(n/bs,count + dim);
68: ++dim;
69: if (bs >= 1) {
70: count[dim] = bs;
71: ++dim;
72: }
73: PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
75: /* Select hyperslab in the file */
76: PetscLayoutGetRange(is->map, &low, NULL);
77: dim = 0;
78: if (timestep >= 0) {
79: offset[dim] = timestep;
80: ++dim;
81: }
82: PetscHDF5IntCast(low/bs,offset + dim);
83: ++dim;
84: if (bs >= 1) {
85: offset[dim] = 0;
86: ++dim;
87: }
88: PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
90: /* Create property list for collective dataset read */
91: PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER));
92: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
93: PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE));
94: #endif
95: /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
97: #if defined(PETSC_USE_64BIT_INDICES)
98: inttype = H5T_NATIVE_LLONG;
99: #else
100: inttype = H5T_NATIVE_INT;
101: #endif
102: PetscMalloc1(n,&ind);
103: PetscStackCallHDF5(H5Dread,(dset_id, inttype, memspace, filespace, plist_id, (void *) ind));
104: ISGeneralSetIndices(is, n, ind, PETSC_OWN_POINTER);
106: /* Close/release resources */
107: if (group != file_id) PetscStackCallHDF5(H5Gclose,(group));
108: PetscStackCallHDF5(H5Pclose,(plist_id));
109: PetscStackCallHDF5(H5Sclose,(filespace));
110: PetscStackCallHDF5(H5Sclose,(memspace));
111: PetscStackCallHDF5(H5Dclose,(dset_id));
112: return(0);
113: }
114: #endif
118: PetscErrorCode ISLoad_Binary(IS is, PetscViewer viewer)
119: {
121: PetscBool isgeneral,skipheader;
122: int fd;
123: PetscInt tr[2],N,ln,*idx;
124: MPI_Request request;
125: MPI_Status status;
126: MPI_Comm comm;
127: PetscMPIInt rank,size,tag;
130: PetscObjectGetComm((PetscObject)is,&comm);
131: MPI_Comm_rank(comm,&rank);
132: MPI_Comm_size(comm,&size);
133: PetscObjectGetNewTag((PetscObject)viewer,&tag);
134: /* force binary viewer to load .info file if it has not yet done so */
135: PetscViewerSetUp(viewer);
136: PetscObjectTypeCompare((PetscObject)is,ISGENERAL,&isgeneral);
137: if (!isgeneral) SETERRQ(comm,PETSC_ERR_ARG_INCOMP,"IS must be of type ISGENERAL to load into it");
138: PetscViewerBinaryGetDescriptor(viewer,&fd);
139: PetscViewerBinaryGetSkipHeader(viewer,&skipheader);
140: if (skipheader) SETERRQ(comm,PETSC_ERR_USER, "Currently no support for binary files without headers");
142: PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
143: if (tr[0] != IS_FILE_CLASSID) SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Not an IS next in file");
145: /* Has IS already had its layout defined */
146: PetscLayoutGetSize(is->map,&N);
147: 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);
148: if (N == -1) {
149: N = tr[1];
150: PetscLayoutSetSize(is->map,N);
151: PetscLayoutSetUp(is->map);
152: }
153: PetscLayoutGetLocalSize(is->map,&ln);
154: PetscMalloc1(ln,&idx);
155: if (!rank) {
156: PetscBinaryRead(fd,idx,ln,PETSC_INT);
158: if (size > 1) {
159: PetscInt *range,n,i,*idxwork;
161: /* read in other chuncks and send to other processors */
162: /* determine maximum chunck owned by other */
163: range = is->map->range;
164: n = 1;
165: for (i=1; i<size; i++) n = PetscMax(n,range[i+1] - range[i]);
167: PetscMalloc1(n,&idxwork);
168: for (i=1; i<size; i++) {
169: n = range[i+1] - range[i];
170: PetscBinaryRead(fd,idxwork,n,PETSC_INT);
171: MPI_Isend(idxwork,n,MPIU_INT,i,tag,comm,&request);
172: MPI_Wait(&request,&status);
173: }
174: PetscFree(idxwork);
175: }
176: } else {
177: MPI_Recv(idx,ln,MPIU_INT,0,tag,comm,&status);
178: }
179: ISGeneralSetIndices(is,ln,idx,PETSC_OWN_POINTER);
180: return(0);
181: }
185: PetscErrorCode ISLoad_Default(IS is, PetscViewer viewer)
186: {
187: PetscBool isbinary;
188: #if defined(PETSC_HAVE_HDF5)
189: PetscBool ishdf5;
190: #endif
194: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
195: #if defined(PETSC_HAVE_HDF5)
196: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
197: #endif
198: if (isbinary) {
199: ISLoad_Binary(is, viewer);
200: #if defined(PETSC_HAVE_HDF5)
201: } else if (ishdf5) {
202: ISLoad_HDF5(is, viewer);
203: #endif
204: }
205: return(0);
206: }