Actual source code: vecio.c
petsc-3.5.4 2015-05-23
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, ×tep);
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, ×tep);
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: }