Actual source code: vecio.c
petsc-3.11.3 2019-06-26
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>
10: #include <petsc/private/vecimpl.h>
11: #include <petsc/private/viewerimpl.h>
12: #include <petscviewerhdf5.h>
14: static PetscErrorCode PetscViewerBinaryReadVecHeader_Private(PetscViewer viewer,PetscInt *rows)
15: {
17: MPI_Comm comm;
18: PetscInt tr[2],type;
21: PetscObjectGetComm((PetscObject)viewer,&comm);
22: /* Read vector header */
23: PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
24: type = tr[0];
25: if (type != VEC_FILE_CLASSID) {
26: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
27: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a vector next in file");
28: }
29: *rows = tr[1];
30: return(0);
31: }
33: #if defined(PETSC_HAVE_MPIIO)
34: static PetscErrorCode VecLoad_Binary_MPIIO(Vec vec, PetscViewer viewer)
35: {
37: PetscMPIInt lsize;
38: PetscScalar *avec;
39: MPI_File mfdes;
40: MPI_Offset off;
43: VecGetArray(vec,&avec);
44: PetscMPIIntCast(vec->map->n,&lsize);
46: PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
47: PetscViewerBinaryGetMPIIOOffset(viewer,&off);
48: off += vec->map->rstart*sizeof(PetscScalar);
49: MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
50: MPIU_File_read_all(mfdes,avec,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
51: PetscViewerBinaryAddMPIIOOffset(viewer,vec->map->N*sizeof(PetscScalar));
53: VecRestoreArray(vec,&avec);
54: VecAssemblyBegin(vec);
55: VecAssemblyEnd(vec);
56: return(0);
57: }
58: #endif
60: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
61: {
62: PetscMPIInt size,rank,tag;
63: int fd;
64: PetscInt i,rows = 0,n,*range,N,bs;
66: PetscBool flag,skipheader;
67: PetscScalar *avec,*avecwork;
68: MPI_Comm comm;
69: MPI_Request request;
70: MPI_Status status;
71: #if defined(PETSC_HAVE_MPIIO)
72: PetscBool useMPIIO;
73: #endif
76: /* force binary viewer to load .info file if it has not yet done so */
77: PetscViewerSetUp(viewer);
78: PetscObjectGetComm((PetscObject)viewer,&comm);
79: MPI_Comm_rank(comm,&rank);
80: MPI_Comm_size(comm,&size);
82: PetscViewerBinaryGetDescriptor(viewer,&fd);
83: PetscViewerBinaryGetSkipHeader(viewer,&skipheader);
84: if (!skipheader) {
85: PetscViewerBinaryReadVecHeader_Private(viewer,&rows);
86: } else {
87: VecType vtype;
88: VecGetType(vec,&vtype);
89: 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");
90: VecGetSize(vec,&N);
91: 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");
92: rows = N;
93: }
94: /* Set Vec sizes,blocksize,and type if not already set. Block size first so that local sizes will be compatible. */
95: PetscOptionsGetInt(NULL,((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flag);
96: if (flag) {
97: VecSetBlockSize(vec, bs);
98: }
99: if (vec->map->n < 0 && vec->map->N < 0) {
100: VecSetSizes(vec,PETSC_DECIDE,rows);
101: }
103: /* If sizes and type already set,check if the vector global size is correct */
104: VecGetSize(vec, &N);
105: if (N != rows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", rows, N);
107: #if defined(PETSC_HAVE_MPIIO)
108: PetscViewerBinaryGetUseMPIIO(viewer,&useMPIIO);
109: if (useMPIIO) {
110: VecLoad_Binary_MPIIO(vec, viewer);
111: return(0);
112: }
113: #endif
115: VecGetLocalSize(vec,&n);
116: PetscObjectGetNewTag((PetscObject)viewer,&tag);
117: VecGetArray(vec,&avec);
118: if (!rank) {
119: PetscBinaryRead(fd,avec,n,PETSC_SCALAR);
121: if (size > 1) {
122: /* read in other chuncks and send to other processors */
123: /* determine maximum chunck owned by other */
124: range = vec->map->range;
125: n = 1;
126: for (i=1; i<size; i++) n = PetscMax(n,range[i+1] - range[i]);
128: PetscMalloc1(n,&avecwork);
129: for (i=1; i<size; i++) {
130: n = range[i+1] - range[i];
131: PetscBinaryRead(fd,avecwork,n,PETSC_SCALAR);
132: MPI_Isend(avecwork,n,MPIU_SCALAR,i,tag,comm,&request);
133: MPI_Wait(&request,&status);
134: }
135: PetscFree(avecwork);
136: }
137: } else {
138: MPI_Recv(avec,n,MPIU_SCALAR,0,tag,comm,&status);
139: }
141: VecRestoreArray(vec,&avec);
142: VecAssemblyBegin(vec);
143: VecAssemblyEnd(vec);
144: return(0);
145: }
147: #if defined(PETSC_HAVE_HDF5)
148: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
149: {
150: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
151: PetscScalar *x;
152: const char *vecname;
156: if (!((PetscObject)xin)->name) SETERRQ(PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
157: #if defined(PETSC_USE_REAL_SINGLE)
158: scalartype = H5T_NATIVE_FLOAT;
159: #elif defined(PETSC_USE_REAL___FLOAT128)
160: #error "HDF5 output with 128 bit floats not supported."
161: #elif defined(PETSC_USE_REAL___FP16)
162: #error "HDF5 output with 16 bit floats not supported."
163: #else
164: scalartype = H5T_NATIVE_DOUBLE;
165: #endif
166: PetscObjectGetName((PetscObject)xin, &vecname);
167: PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void**)&x);
168: VecSetUp(xin);
169: VecReplaceArray(xin, x);
170: return(0);
171: }
172: #endif
174: #if defined(PETSC_HAVE_ADIOS)
175: #include <adios.h>
176: #include <adios_read.h>
177: #include <petsc/private/vieweradiosimpl.h>
178: #include <petsc/private/viewerimpl.h>
180: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
181: {
182: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
183: PetscErrorCode ierr;
184: PetscScalar *x;
185: PetscInt Nfile,N,rstart,n;
186: uint64_t N_t,rstart_t;
187: const char *vecname;
188: ADIOS_VARINFO *v;
189: ADIOS_SELECTION *sel;
192: PetscObjectGetName((PetscObject) xin, &vecname);
194: v = adios_inq_var(adios->adios_fp, vecname);
195: if (v->ndim != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%D)", v->ndim);
196: Nfile = (PetscInt) v->dims[0];
198: /* Set Vec sizes,blocksize,and type if not already set */
199: if ((xin)->map->n < 0 && (xin)->map->N < 0) {
200: VecSetSizes(xin, PETSC_DECIDE, Nfile);
201: }
202: /* If sizes and type already set,check if the vector global size is correct */
203: VecGetSize(xin, &N);
204: VecGetLocalSize(xin, &n);
205: if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
206:
207: VecGetOwnershipRange(xin,&rstart,NULL);
208: rstart_t = rstart;
209: N_t = n;
210: sel = adios_selection_boundingbox (v->ndim, &rstart_t, &N_t);
211: VecGetArray(xin,&x);
212: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
213: adios_perform_reads (adios->adios_fp, 1);
214: VecRestoreArray(xin,&x);
215: adios_selection_delete(sel);
217: return(0);
218: }
219: #endif
221: #if defined(PETSC_HAVE_ADIOS2)
222: #include <adios2_c.h>
223: #include <petsc/private/vieweradios2impl.h>
224: #include <petsc/private/viewerimpl.h>
226: PetscErrorCode VecLoad_ADIOS2(Vec xin, PetscViewer viewer)
227: {
228: PetscViewer_ADIOS2 *adios2 = (PetscViewer_ADIOS2*)viewer->data;
229: PetscErrorCode ierr;
230: PetscScalar *x;
231: PetscInt Nfile,N,rstart,n;
232: const char *vecname;
235: PetscObjectGetName((PetscObject) xin, &vecname);
237: /* Set Vec sizes,blocksize,and type if not already set */
238: if ((xin)->map->n < 0 && (xin)->map->N < 0) {
239: VecSetSizes(xin, PETSC_DECIDE, Nfile);
240: }
241: /* If sizes and type already set,check if the vector global size is correct */
242: VecGetSize(xin, &N);
243: VecGetLocalSize(xin, &n);
244: if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
246: VecGetOwnershipRange(xin,&rstart,NULL);
247: VecGetArray(xin,&x);
248: VecRestoreArray(xin,&x);
249: return(0);
250: }
251: #endif
253: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
254: {
256: PetscBool isbinary;
257: #if defined(PETSC_HAVE_HDF5)
258: PetscBool ishdf5;
259: #endif
260: #if defined(PETSC_HAVE_ADIOS)
261: PetscBool isadios;
262: #endif
263: #if defined(PETSC_HAVE_ADIOS2)
264: PetscBool isadios2;
265: #endif
268: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
269: #if defined(PETSC_HAVE_HDF5)
270: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
271: #endif
272: #if defined(PETSC_HAVE_ADIOS)
273: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
274: #endif
275: #if defined(PETSC_HAVE_ADIOS2)
276: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS2,&isadios2);
277: #endif
279: #if defined(PETSC_HAVE_HDF5)
280: if (ishdf5) {
281: if (!((PetscObject)newvec)->name) {
282: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
283: 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()");
284: }
285: VecLoad_HDF5(newvec, viewer);
286: } else
287: #endif
288: #if defined(PETSC_HAVE_ADIOS)
289: if (isadios) {
290: if (!((PetscObject)newvec)->name) {
291: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
292: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since ADIOS format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
293: }
294: VecLoad_ADIOS(newvec, viewer);
295: } else
296: #endif
297: #if defined(PETSC_HAVE_ADIOS2)
298: if (isadios2) {
299: if (!((PetscObject)newvec)->name) {
300: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
301: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since ADIOS2 format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
302: }
303: VecLoad_ADIOS2(newvec, viewer);
304: } else
305: #endif
306: {
307: VecLoad_Binary(newvec, viewer);
308: }
309: return(0);
310: }
312: /*@
313: VecChop - Set all values in the vector with an absolute value less than the tolerance to zero
315: Input Parameters:
316: + v - The vector
317: - tol - The zero tolerance
319: Output Parameters:
320: . v - The chopped vector
322: Level: intermediate
324: .seealso: VecCreate(), VecSet()
325: @*/
326: PetscErrorCode VecChop(Vec v, PetscReal tol)
327: {
328: PetscScalar *a;
329: PetscInt n, i;
333: VecGetLocalSize(v, &n);
334: VecGetArray(v, &a);
335: for (i = 0; i < n; ++i) {
336: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
337: }
338: VecRestoreArray(v, &a);
339: return(0);
340: }