Actual source code: vecio.c
petsc-3.12.2 2019-11-22
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 <petsclayouthdf5.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,NULL,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,NULL,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,*arr;
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); /* VecSetSizes might have not been called so ensure ops->create has been called */
169: if (!xin->ops->replacearray) {
170: VecGetArray(xin,&arr);
171: PetscArraycpy(arr,x,xin->map->n);
172: PetscFree(x);
173: VecRestoreArray(xin,&arr);
174: } else {
175: VecReplaceArray(xin,x);
176: }
177: return(0);
178: }
179: #endif
181: #if defined(PETSC_HAVE_ADIOS)
182: #include <adios.h>
183: #include <adios_read.h>
184: #include <petsc/private/vieweradiosimpl.h>
185: #include <petsc/private/viewerimpl.h>
187: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
188: {
189: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
190: PetscErrorCode ierr;
191: PetscScalar *x;
192: PetscInt Nfile,N,rstart,n;
193: uint64_t N_t,rstart_t;
194: const char *vecname;
195: ADIOS_VARINFO *v;
196: ADIOS_SELECTION *sel;
199: PetscObjectGetName((PetscObject) xin, &vecname);
201: v = adios_inq_var(adios->adios_fp, vecname);
202: if (v->ndim != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%D)", v->ndim);
203: Nfile = (PetscInt) v->dims[0];
205: /* Set Vec sizes,blocksize,and type if not already set */
206: if ((xin)->map->n < 0 && (xin)->map->N < 0) {
207: VecSetSizes(xin, PETSC_DECIDE, Nfile);
208: }
209: /* If sizes and type already set,check if the vector global size is correct */
210: VecGetSize(xin, &N);
211: VecGetLocalSize(xin, &n);
212: if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
213:
214: VecGetOwnershipRange(xin,&rstart,NULL);
215: rstart_t = rstart;
216: N_t = n;
217: sel = adios_selection_boundingbox (v->ndim, &rstart_t, &N_t);
218: VecGetArray(xin,&x);
219: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
220: adios_perform_reads (adios->adios_fp, 1);
221: VecRestoreArray(xin,&x);
222: adios_selection_delete(sel);
224: return(0);
225: }
226: #endif
228: #if defined(PETSC_HAVE_ADIOS2)
229: #include <adios2_c.h>
230: #include <petsc/private/vieweradios2impl.h>
231: #include <petsc/private/viewerimpl.h>
233: PetscErrorCode VecLoad_ADIOS2(Vec xin, PetscViewer viewer)
234: {
235: PetscViewer_ADIOS2 *adios2 = (PetscViewer_ADIOS2*)viewer->data;
236: PetscErrorCode ierr;
237: PetscScalar *x;
238: PetscInt Nfile,N,rstart,n;
239: const char *vecname;
242: PetscObjectGetName((PetscObject) xin, &vecname);
244: /* Set Vec sizes,blocksize,and type if not already set */
245: if ((xin)->map->n < 0 && (xin)->map->N < 0) {
246: VecSetSizes(xin, PETSC_DECIDE, Nfile);
247: }
248: /* If sizes and type already set,check if the vector global size is correct */
249: VecGetSize(xin, &N);
250: VecGetLocalSize(xin, &n);
251: if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
253: VecGetOwnershipRange(xin,&rstart,NULL);
254: VecGetArray(xin,&x);
255: VecRestoreArray(xin,&x);
256: return(0);
257: }
258: #endif
260: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
261: {
263: PetscBool isbinary;
264: #if defined(PETSC_HAVE_HDF5)
265: PetscBool ishdf5;
266: #endif
267: #if defined(PETSC_HAVE_ADIOS)
268: PetscBool isadios;
269: #endif
270: #if defined(PETSC_HAVE_ADIOS2)
271: PetscBool isadios2;
272: #endif
275: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
276: #if defined(PETSC_HAVE_HDF5)
277: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
278: #endif
279: #if defined(PETSC_HAVE_ADIOS)
280: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
281: #endif
282: #if defined(PETSC_HAVE_ADIOS2)
283: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS2,&isadios2);
284: #endif
286: #if defined(PETSC_HAVE_HDF5)
287: if (ishdf5) {
288: if (!((PetscObject)newvec)->name) {
289: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
290: 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()");
291: }
292: VecLoad_HDF5(newvec, viewer);
293: } else
294: #endif
295: #if defined(PETSC_HAVE_ADIOS)
296: if (isadios) {
297: if (!((PetscObject)newvec)->name) {
298: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
299: 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()");
300: }
301: VecLoad_ADIOS(newvec, viewer);
302: } else
303: #endif
304: #if defined(PETSC_HAVE_ADIOS2)
305: if (isadios2) {
306: if (!((PetscObject)newvec)->name) {
307: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
308: 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()");
309: }
310: VecLoad_ADIOS2(newvec, viewer);
311: } else
312: #endif
313: {
314: VecLoad_Binary(newvec, viewer);
315: }
316: return(0);
317: }
319: /*@
320: VecChop - Set all values in the vector with an absolute value less than the tolerance to zero
322: Input Parameters:
323: + v - The vector
324: - tol - The zero tolerance
326: Output Parameters:
327: . v - The chopped vector
329: Level: intermediate
331: .seealso: VecCreate(), VecSet()
332: @*/
333: PetscErrorCode VecChop(Vec v, PetscReal tol)
334: {
335: PetscScalar *a;
336: PetscInt n, i;
340: VecGetLocalSize(v, &n);
341: VecGetArray(v, &a);
342: for (i = 0; i < n; ++i) {
343: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
344: }
345: VecRestoreArray(v, &a);
346: return(0);
347: }