Actual source code: isio.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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, &timestep);
 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: }