Actual source code: vecio.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  2: /*
  3:    This file contains simple binary input routines for vectors.  The
  4:    analogous output routines are within each vector implementation's
  5:    VecView (with viewer types PETSCVIEWERBINARY)
  6:  */

  8: #include <petscsys.h>
  9: #include <petscvec.h>         /*I  "petscvec.h"  I*/
 10: #include <petsc-private/vecimpl.h>
 11: #include <petscmat.h> /* so that MAT_FILE_CLASSID is defined */
 12: #include <petscviewerhdf5.h>

 16: static PetscErrorCode PetscViewerBinaryReadVecHeader_Private(PetscViewer viewer,PetscInt *rows)
 17: {
 19:   MPI_Comm       comm;
 20:   PetscInt       tr[2],type;

 23:   PetscObjectGetComm((PetscObject)viewer,&comm);
 24:   /* Read vector header */
 25:   PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);
 26:   type = tr[0];
 27:   if (type != VEC_FILE_CLASSID) {
 28:     PetscLogEventEnd(VEC_Load,viewer,0,0,0);
 29:     if (type == MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix is next in file, not a vector as you requested");
 30:     else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a vector next in file");
 31:   }
 32:   *rows = tr[1];
 33:   return(0);
 34: }

 36: #if defined(PETSC_HAVE_MPIIO)
 39: static PetscErrorCode VecLoad_Binary_MPIIO(Vec vec, PetscViewer viewer)
 40: {
 42:   PetscMPIInt    lsize;
 43:   PetscScalar    *avec;
 44:   MPI_File       mfdes;
 45:   MPI_Offset     off;

 48:   VecGetArray(vec,&avec);
 49:   PetscMPIIntCast(vec->map->n,&lsize);

 51:   PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
 52:   PetscViewerBinaryGetMPIIOOffset(viewer,&off);
 53:   off += vec->map->rstart*sizeof(PetscScalar);
 54:   MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
 55:   MPIU_File_read_all(mfdes,avec,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
 56:   PetscViewerBinaryAddMPIIOOffset(viewer,vec->map->N*sizeof(PetscScalar));

 58:   VecRestoreArray(vec,&avec);
 59:   VecAssemblyBegin(vec);
 60:   VecAssemblyEnd(vec);
 61:   return(0);
 62: }
 63: #endif

 67: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
 68: {
 69:   PetscMPIInt    size,rank,tag;
 70:   int            fd;
 71:   PetscInt       i,rows = 0,n,*range,N,bs;
 73:   PetscBool      flag,skipheader;
 74:   PetscScalar    *avec,*avecwork;
 75:   MPI_Comm       comm;
 76:   MPI_Request    request;
 77:   MPI_Status     status;
 78: #if defined(PETSC_HAVE_MPIIO)
 79:   PetscBool      useMPIIO;
 80: #endif

 83:   PetscObjectGetComm((PetscObject)viewer,&comm);
 84:   MPI_Comm_rank(comm,&rank);
 85:   MPI_Comm_size(comm,&size);

 87:   PetscViewerBinaryGetDescriptor(viewer,&fd);
 88:   PetscViewerBinaryGetSkipHeader(viewer,&skipheader);
 89:   if (!skipheader) {
 90:     PetscViewerBinaryReadVecHeader_Private(viewer,&rows);
 91:   } else {
 92:     VecType vtype;
 93:     VecGetType(vec,&vtype);
 94:     if (!vtype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER, "Vector binary file header was skipped, thus the user must specify the type and length of input vector");
 95:     VecGetSize(vec,&N);
 96:     if (!N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER, "Vector binary file header was skipped, thus the user must specify the length of input vector");
 97:     rows = N;
 98:   }
 99:   /* Set Vec sizes,blocksize,and type if not already set. Block size first so that local sizes will be compatible. */
100:   PetscOptionsGetInt(((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flag);
101:   if (flag) {
102:     VecSetBlockSize(vec, bs);
103:   }
104:   if (vec->map->n < 0 && vec->map->N < 0) {
105:     VecSetSizes(vec,PETSC_DECIDE,rows);
106:   }

108:   /* If sizes and type already set,check if the vector global size is correct */
109:   VecGetSize(vec, &N);
110:   if (N != rows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%d) then input vector (%d)", rows, N);

112: #if defined(PETSC_HAVE_MPIIO)
113:   PetscViewerBinaryGetMPIIO(viewer,&useMPIIO);
114:   if (useMPIIO) {
115:     VecLoad_Binary_MPIIO(vec, viewer);
116:     return(0);
117:   }
118: #endif

120:   VecGetLocalSize(vec,&n);
121:   PetscObjectGetNewTag((PetscObject)viewer,&tag);
122:   VecGetArray(vec,&avec);
123:   if (!rank) {
124:     PetscBinaryRead(fd,avec,n,PETSC_SCALAR);

126:     if (size > 1) {
127:       /* read in other chuncks and send to other processors */
128:       /* determine maximum chunck owned by other */
129:       range = vec->map->range;
130:       n = 1;
131:       for (i=1; i<size; i++) n = PetscMax(n,range[i+1] - range[i]);

133:       PetscMalloc1(n,&avecwork);
134:       for (i=1; i<size; i++) {
135:         n    = range[i+1] - range[i];
136:         PetscBinaryRead(fd,avecwork,n,PETSC_SCALAR);
137:         MPI_Isend(avecwork,n,MPIU_SCALAR,i,tag,comm,&request);
138:         MPI_Wait(&request,&status);
139:       }
140:       PetscFree(avecwork);
141:     }
142:   } else {
143:     MPI_Recv(avec,n,MPIU_SCALAR,0,tag,comm,&status);
144:   }

146:   VecRestoreArray(vec,&avec);
147:   VecAssemblyBegin(vec);
148:   VecAssemblyEnd(vec);
149:   return(0);
150: }

152: #if defined(PETSC_HAVE_HDF5)
155: PetscErrorCode PetscViewerHDF5OpenGroup(PetscViewer viewer, hid_t *fileId, hid_t *groupId)
156: {
157:   hid_t          file_id, group;
158:   htri_t         found;
159:   const char     *groupName = NULL;

163:   PetscViewerHDF5GetFileId(viewer, &file_id);
164:   PetscViewerHDF5GetGroup(viewer, &groupName);
165:   /* Open group */
166:   if (groupName) {
167:     PetscBool root;

169:     PetscStrcmp(groupName, "/", &root);
170:     found = H5Lexists(file_id, groupName, H5P_DEFAULT);
171:     if (!root && (found <= 0)) {
172: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
173:       group = H5Gcreate2(file_id, groupName, 0, H5P_DEFAULT, H5P_DEFAULT);
174: #else /* deprecated HDF5 1.6 API */
175:       group = H5Gcreate(file_id, groupName, 0);
176: #endif
177:       if (group < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Could not create group %s", groupName);
178:       H5Gclose(group);
179:     }
180: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
181:     group = H5Gopen2(file_id, groupName, H5P_DEFAULT);
182: #else
183:     group = H5Gopen(file_id, groupName);
184: #endif
185:     if (group < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Could not open group %s", groupName);
186:   } else group = file_id;

188:   *fileId  = file_id;
189:   *groupId = group;
190:   return(0);
191: }

195: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
196: {
197:   hid_t          file_id, group, dset_id, filespace;
198:   hsize_t        rdim, dim;
199:   hsize_t        dims[4];
200:   herr_t         status;
201:   PetscInt       bsInd, lenInd, timestep;

205:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
206:   PetscViewerHDF5GetTimestep(viewer, &timestep);
207: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
208:   dset_id = H5Dopen2(group, name, H5P_DEFAULT);
209: #else
210:   dset_id = H5Dopen(group, name);
211: #endif
212:   if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Could not H5Dopen() with PetscObject named %s", name);
213:   filespace = H5Dget_space(dset_id);
214:   if (filespace == -1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Could not H5Dget_space()");
215:   dim = 0;
216:   if (timestep >= 0) ++dim;
217:   ++dim; /* length in blocks */
218:   ++dim; /* block size */
219: #if defined(PETSC_USE_COMPLEX)
220:   ++dim;
221: #endif
222:   rdim = H5Sget_simple_extent_dims(filespace, dims, NULL);
223: #if defined(PETSC_USE_COMPLEX)
224:   bsInd = rdim-2;
225: #else
226:   bsInd = rdim-1;
227: #endif
228:   lenInd = timestep >= 0 ? 1 : 0;
229:   if (rdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected", rdim, dim);
230:   /* Close/release resources */
231:   status = H5Sclose(filespace);CHKERRQ(status);
232:   status = H5Dclose(dset_id);CHKERRQ(status);
233:   if (group != file_id) {status = H5Gclose(group);CHKERRQ(status);}
234:   if (bs) *bs = (PetscInt) dims[bsInd];
235:   if (N)  *N  = (PetscInt) dims[lenInd]*dims[bsInd];
236:   return(0);
237: }

241: /*
242:      This should handle properly the cases where PetscInt is 32 or 64 and hsize_t is 32 or 64. These means properly casting with
243:    checks back and forth between the two types of variables.
244: */
245: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
246: {
247:   hid_t          file_id, group, dset_id, filespace, memspace, plist_id;
248:   hsize_t        rdim, dim;
249:   hsize_t        dims[4], count[4], offset[4];
250:   herr_t         status;
251:   PetscInt       n, N, bs = 1, bsInd, lenInd, low, timestep;
252:   PetscScalar    *x;
253:   hid_t          scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
254:   const char     *vecname;

258: #if defined(PETSC_USE_REAL_SINGLE)
259:   scalartype = H5T_NATIVE_FLOAT;
260: #elif defined(PETSC_USE_REAL___FLOAT128)
261: #error "HDF5 output with 128 bit floats not supported."
262: #else
263:   scalartype = H5T_NATIVE_DOUBLE;
264: #endif

266:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
267:   PetscViewerHDF5GetTimestep(viewer, &timestep);
268:   VecGetBlockSize(xin,&bs);
269:   /* Create the dataset with default properties and close filespace */
270:   PetscObjectGetName((PetscObject)xin,&vecname);
271: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
272:   dset_id = H5Dopen2(group, vecname, H5P_DEFAULT);
273: #else
274:   dset_id = H5Dopen(group, vecname);
275: #endif
276:   if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Dopen() with Vec named %s",vecname);
277:   /* Retrieve the dataspace for the dataset */
278:   filespace = H5Dget_space(dset_id);
279:   if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Dget_space()");
280:   dim = 0;
281:   if (timestep >= 0) ++dim;
282:   ++dim;
283:   if (bs >= 1) ++dim;
284: #if defined(PETSC_USE_COMPLEX)
285:   ++dim;
286: #endif
287:   rdim = H5Sget_simple_extent_dims(filespace, dims, NULL);
288: #if defined(PETSC_USE_COMPLEX)
289:   bsInd = rdim-2;
290: #else
291:   bsInd = rdim-1;
292: #endif
293:   lenInd = timestep >= 0 ? 1 : 0;
294:   if (rdim != dim) {
295:     if (rdim == dim+1 && bs == -1) bs = dims[bsInd];
296:     else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim);
297:   } else if (bs >= 1 && bs != (PetscInt) dims[bsInd]) {
298:     VecSetBlockSize(xin, dims[bsInd]);
299:     if (ierr) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Block size %d specified for vector does not match blocksize in file %d",bs,dims[bsInd]);
300:     bs = dims[bsInd];
301:   }

303:   /* Set Vec sizes,blocksize,and type if not already set */
304:   if ((xin)->map->n < 0 && (xin)->map->N < 0) {
305:     VecSetSizes(xin, PETSC_DECIDE, dims[lenInd]*bs);
306:   }
307:   /* If sizes and type already set,check if the vector global size is correct */
308:   VecGetSize(xin, &N);
309:   if (N/bs != (PetscInt) dims[lenInd]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%d) then input vector (%d)", (PetscInt) dims[lenInd], N/bs);

311:   /* Each process defines a dataset and reads it from the hyperslab in the file */
312:   VecGetLocalSize(xin, &n);
313:   dim  = 0;
314:   if (timestep >= 0) {
315:     count[dim] = 1;
316:     ++dim;
317:   }
318:   PetscHDF5IntCast(n/bs,count + dim);
319:   ++dim;
320:   if (bs >= 1) {
321:     count[dim] = bs;
322:     ++dim;
323:   }
324: #if defined(PETSC_USE_COMPLEX)
325:   count[dim] = 2;
326:   ++dim;
327: #endif
328:   memspace = H5Screate_simple(dim, count, NULL);
329:   if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Screate_simple()");

331:   /* Select hyperslab in the file */
332:   VecGetOwnershipRange(xin, &low, NULL);
333:   dim  = 0;
334:   if (timestep >= 0) {
335:     offset[dim] = timestep;
336:     ++dim;
337:   }
338:   PetscHDF5IntCast(low/bs,offset + dim);
339:   ++dim;
340:   if (bs >= 1) {
341:     offset[dim] = 0;
342:     ++dim;
343:   }
344: #if defined(PETSC_USE_COMPLEX)
345:   offset[dim] = 0;
346:   ++dim;
347: #endif
348:   status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status);

350:   /* Create property list for collective dataset read */
351:   plist_id = H5Pcreate(H5P_DATASET_XFER);
352:   if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Could not H5Pcreate()");
353: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
354:   status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status);
355: #endif
356:   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */

358:   VecGetArray(xin, &x);
359:   status = H5Dread(dset_id, scalartype, memspace, filespace, plist_id, x);CHKERRQ(status);
360:   VecRestoreArray(xin, &x);

362:   /* Close/release resources */
363:   if (group != file_id) {
364:     status = H5Gclose(group);CHKERRQ(status);
365:   }
366:   status = H5Pclose(plist_id);CHKERRQ(status);
367:   status = H5Sclose(filespace);CHKERRQ(status);
368:   status = H5Sclose(memspace);CHKERRQ(status);
369:   status = H5Dclose(dset_id);CHKERRQ(status);

371:   VecAssemblyBegin(xin);
372:   VecAssemblyEnd(xin);
373:   return(0);
374: }
375: #endif


380: PetscErrorCode  VecLoad_Default(Vec newvec, PetscViewer viewer)
381: {
383:   PetscBool      isbinary;
384: #if defined(PETSC_HAVE_HDF5)
385:   PetscBool      ishdf5;
386: #endif

389:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
390: #if defined(PETSC_HAVE_HDF5)
391:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
392: #endif

394: #if defined(PETSC_HAVE_HDF5)
395:   if (ishdf5) {
396:     if (!((PetscObject)newvec)->name) {
397:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since HDF5 format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
398:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
399:     }
400:     VecLoad_HDF5(newvec, viewer);
401:   } else
402: #endif
403:   {
404:     VecLoad_Binary(newvec, viewer);
405:   }
406:   return(0);
407: }

411: /*@
412:   VecChop - Set all values in the vector less than the tolerance to zero

414:   Input Parameters:
415: + v   - The vector
416: - tol - The zero tolerance

418:   Output Parameters:
419: . v - The chopped vector

421:   Level: intermediate

423: .seealso: VecCreate(), VecSet()
424: @*/
425: PetscErrorCode VecChop(Vec v, PetscReal tol)
426: {
427:   PetscScalar    *a;
428:   PetscInt       n, i;

432:   VecGetLocalSize(v, &n);
433:   VecGetArray(v, &a);
434:   for (i = 0; i < n; ++i) {
435:     if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
436:   }
437:   VecRestoreArray(v, &a);
438:   return(0);
439: }