Actual source code: vecio.c

petsc-master 2015-01-26
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;
 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:   PetscViewerBinaryReadVecHeader_Private(viewer,&rows);
 89:   /* Set Vec sizes,blocksize,and type if not already set. Block size first so that local sizes will be compatible. */
 90:   PetscOptionsGetInt(((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flag);
 91:   if (flag) {
 92:     VecSetBlockSize(vec, bs);
 93:   }
 94:   if (vec->map->n < 0 && vec->map->N < 0) {
 95:     VecSetSizes(vec,PETSC_DECIDE,rows);
 96:   }

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

102: #if defined(PETSC_HAVE_MPIIO)
103:   PetscViewerBinaryGetMPIIO(viewer,&useMPIIO);
104:   if (useMPIIO) {
105:     VecLoad_Binary_MPIIO(vec, viewer);
106:     return(0);
107:   }
108: #endif

110:   VecGetLocalSize(vec,&n);
111:   PetscObjectGetNewTag((PetscObject)viewer,&tag);
112:   VecGetArray(vec,&avec);
113:   if (!rank) {
114:     PetscBinaryRead(fd,avec,n,PETSC_SCALAR);

116:     if (size > 1) {
117:       /* read in other chuncks and send to other processors */
118:       /* determine maximum chunck owned by other */
119:       range = vec->map->range;
120:       n = 1;
121:       for (i=1; i<size; i++) n = PetscMax(n,range[i+1] - range[i]);

123:       PetscMalloc1(n,&avecwork);
124:       for (i=1; i<size; i++) {
125:         n    = range[i+1] - range[i];
126:         PetscBinaryRead(fd,avecwork,n,PETSC_SCALAR);
127:         MPI_Isend(avecwork,n,MPIU_SCALAR,i,tag,comm,&request);
128:         MPI_Wait(&request,&status);
129:       }
130:       PetscFree(avecwork);
131:     }
132:   } else {
133:     MPI_Recv(avec,n,MPIU_SCALAR,0,tag,comm,&status);
134:   }

136:   VecRestoreArray(vec,&avec);
137:   VecAssemblyBegin(vec);
138:   VecAssemblyEnd(vec);
139:   return(0);
140: }

142: #if defined(PETSC_HAVE_HDF5)
145: PetscErrorCode PetscViewerHDF5OpenGroup(PetscViewer viewer, hid_t *fileId, hid_t *groupId)
146: {
147:   hid_t          file_id, group;
148:   htri_t         found;
149:   const char     *groupName = NULL;

153:   PetscViewerHDF5GetFileId(viewer, &file_id);
154:   PetscViewerHDF5GetGroup(viewer, &groupName);
155:   /* Open group */
156:   if (groupName) {
157:     PetscBool root;

159:     PetscStrcmp(groupName, "/", &root);
160:     PetscStackCall("H5Lexists",found = H5Lexists(file_id, groupName, H5P_DEFAULT));
161:     if (!root && (found <= 0)) {
162: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
163:       PetscStackCallHDF5Return(group,H5Gcreate2,(file_id, groupName, 0, H5P_DEFAULT, H5P_DEFAULT));
164: #else /* deprecated HDF5 1.6 API */
165:       PetscStackCallHDF5Return(group,H5Gcreate,(file_id, groupName, 0));
166: #endif
167:       PetscStackCallHDF5(H5Gclose,(group));
168:     }
169: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
170:     PetscStackCallHDF5Return(group,H5Gopen2,(file_id, groupName, H5P_DEFAULT));
171: #else
172:     PetscStackCallHDF5Return(group,H5Gopen,file_id, groupName));
173: #endif
174:   } else group = file_id;

176:   *fileId  = file_id;
177:   *groupId = group;
178:   return(0);
179: }

183: PetscErrorCode PetscViewerHDF5ReadSizes(PetscViewer viewer, const char name[], PetscInt *bs, PetscInt *N)
184: {
185:   hid_t          file_id, group, dset_id, filespace;
186:   hsize_t        rdim, dim;
187:   hsize_t        dims[4];
188:   PetscInt       bsInd, lenInd, timestep;

192:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
193:   PetscViewerHDF5GetTimestep(viewer, &timestep);
194: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
195:   PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, name, H5P_DEFAULT));
196: #else
197:   PetscStackCallHDF5Return(dset_id,H5Dopen,(group, name));
198: #endif
199:   PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
200:   dim = 0;
201:   if (timestep >= 0) ++dim;
202:   ++dim; /* length in blocks */
203:   ++dim; /* block size */
204: #if defined(PETSC_USE_COMPLEX)
205:   ++dim;
206: #endif
207:   PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL));
208: #if defined(PETSC_USE_COMPLEX)
209:   bsInd = rdim-2;
210: #else
211:   bsInd = rdim-1;
212: #endif
213:   lenInd = timestep >= 0 ? 1 : 0;
214:   if (rdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected", rdim, dim);
215:   /* Close/release resources */
216:   PetscStackCallHDF5(H5Sclose,(filespace));
217:   PetscStackCallHDF5(H5Dclose,(dset_id));
218:   if (group != file_id) PetscStackCallHDF5(H5Gclose,(group));
219:   if (bs) *bs = (PetscInt) dims[bsInd];
220:   if (N)  *N  = (PetscInt) dims[lenInd]*dims[bsInd];
221:   return(0);
222: }

226: /*
227:      This should handle properly the cases where PetscInt is 32 or 64 and hsize_t is 32 or 64. These means properly casting with
228:    checks back and forth between the two types of variables.
229: */
230: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
231: {
232:   hid_t          file_id, group, dset_id, filespace, memspace, plist_id;
233:   hsize_t        rdim, dim;
234:   hsize_t        dims[4], count[4], offset[4];
235:   PetscInt       n, N, bs = 1, bsInd, lenInd, low, timestep;
236:   PetscScalar    *x;
237:   const char     *vecname;

241:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
242:   PetscViewerHDF5GetTimestep(viewer, &timestep);
243:   VecGetBlockSize(xin,&bs);
244:   /* Create the dataset with default properties and close filespace */
245:   PetscObjectGetName((PetscObject)xin,&vecname);
246: #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
247:   PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT));
248: #else
249:   PetscStackCallHDF5Return(dset_id,H5Dopen,(group, vecname));
250: #endif
251:   /* Retrieve the dataspace for the dataset */
252:   PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
253:   dim = 0;
254:   if (timestep >= 0) ++dim;
255:   ++dim;
256:   if (bs >= 1) ++dim;
257: #if defined(PETSC_USE_COMPLEX)
258:   ++dim;
259: #endif
260:   PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL));
261: #if defined(PETSC_USE_COMPLEX)
262:   bsInd = rdim-2;
263: #else
264:   bsInd = rdim-1;
265: #endif
266:   lenInd = timestep >= 0 ? 1 : 0;
267:   if (rdim != dim) {
268:     if (rdim == dim+1 && bs == -1) bs = dims[bsInd];
269:     else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file %d not %d as expected",rdim,dim);
270:   } else if (bs >= 1 && bs != (PetscInt) dims[bsInd]) {
271:     VecSetBlockSize(xin, dims[bsInd]);
272:     bs = dims[bsInd];
273:   }

275:   /* Set Vec sizes,blocksize,and type if not already set */
276:   if ((xin)->map->n < 0 && (xin)->map->N < 0) {
277:     VecSetSizes(xin, PETSC_DECIDE, dims[lenInd]*bs);
278:   }
279:   /* If sizes and type already set,check if the vector global size is correct */
280:   VecGetSize(xin, &N);
281:   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);

283:   /* Each process defines a dataset and reads it from the hyperslab in the file */
284:   VecGetLocalSize(xin, &n);
285:   dim  = 0;
286:   if (timestep >= 0) {
287:     count[dim] = 1;
288:     ++dim;
289:   }
290:   PetscHDF5IntCast(n/bs,count + dim);
291:   ++dim;
292:   if (bs >= 1) {
293:     count[dim] = bs;
294:     ++dim;
295:   }
296: #if defined(PETSC_USE_COMPLEX)
297:   count[dim] = 2;
298:   ++dim;
299: #endif
300:   PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));

302:   /* Select hyperslab in the file */
303:   VecGetOwnershipRange(xin, &low, NULL);
304:   dim  = 0;
305:   if (timestep >= 0) {
306:     offset[dim] = timestep;
307:     ++dim;
308:   }
309:   PetscHDF5IntCast(low/bs,offset + dim);
310:   ++dim;
311:   if (bs >= 1) {
312:     offset[dim] = 0;
313:     ++dim;
314:   }
315: #if defined(PETSC_USE_COMPLEX)
316:   offset[dim] = 0;
317:   ++dim;
318: #endif
319:   PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));

321:   /* Create property list for collective dataset read */
322:   PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER));
323: #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
324:   PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE));
325: #endif
326:   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */

328:   VecGetArray(xin, &x);
329:   PetscStackCallHDF5(H5Dread,(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x));
330:   VecRestoreArray(xin, &x);

332:   /* Close/release resources */
333:   if (group != file_id) {
334:     PetscStackCallHDF5(H5Gclose,(group));
335:   }
336:   PetscStackCallHDF5(H5Pclose,(plist_id));
337:   PetscStackCallHDF5(H5Sclose,(filespace));
338:   PetscStackCallHDF5(H5Sclose,(memspace));
339:   PetscStackCallHDF5(H5Dclose,(dset_id));

341:   VecAssemblyBegin(xin);
342:   VecAssemblyEnd(xin);
343:   return(0);
344: }
345: #endif


350: PetscErrorCode  VecLoad_Default(Vec newvec, PetscViewer viewer)
351: {
353:   PetscBool      isbinary;
354: #if defined(PETSC_HAVE_HDF5)
355:   PetscBool      ishdf5;
356: #endif

359:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
360: #if defined(PETSC_HAVE_HDF5)
361:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
362: #endif

364: #if defined(PETSC_HAVE_HDF5)
365:   if (ishdf5) {
366:     if (!((PetscObject)newvec)->name) {
367:       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()");
368:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
369:     }
370:     VecLoad_HDF5(newvec, viewer);
371:   } else
372: #endif
373:   {
374:     VecLoad_Binary(newvec, viewer);
375:   }
376:   return(0);
377: }

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

384:   Input Parameters:
385: + v   - The vector
386: - tol - The zero tolerance

388:   Output Parameters:
389: . v - The chopped vector

391:   Level: intermediate

393: .seealso: VecCreate(), VecSet()
394: @*/
395: PetscErrorCode VecChop(Vec v, PetscReal tol)
396: {
397:   PetscScalar    *a;
398:   PetscInt       n, i;

402:   VecGetLocalSize(v, &n);
403:   VecGetArray(v, &a);
404:   for (i = 0; i < n; ++i) {
405:     if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
406:   }
407:   VecRestoreArray(v, &a);
408:   return(0);
409: }