Actual source code: vecio.c
petsc-3.13.6 2020-09-29
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: PetscErrorCode VecView_Binary(Vec vec,PetscViewer viewer)
15: {
16: PetscErrorCode ierr;
17: PetscBool skipHeader;
18: PetscLayout map;
19: PetscInt tr[2],n,s,N;
20: const PetscScalar *xarray;
24: PetscViewerSetUp(viewer);
25: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
27: VecGetLayout(vec,&map);
28: PetscLayoutGetLocalSize(map,&n);
29: PetscLayoutGetRange(map,&s,NULL);
30: PetscLayoutGetSize(map,&N);
32: tr[0] = VEC_FILE_CLASSID; tr[1] = N;
33: if (!skipHeader) {PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT);}
35: VecGetArrayRead(vec,&xarray);
36: PetscViewerBinaryWriteAll(viewer,xarray,n,s,N,PETSC_SCALAR);
37: VecRestoreArrayRead(vec,&xarray);
39: { /* write to the viewer's .info file */
40: FILE *info;
41: PetscMPIInt rank;
42: PetscViewerFormat format;
43: const char *name,*pre;
45: PetscViewerBinaryGetInfoPointer(viewer,&info);
46: MPI_Comm_rank(PetscObjectComm((PetscObject)vec),&rank);
48: PetscViewerGetFormat(viewer,&format);
49: if (format == PETSC_VIEWER_BINARY_MATLAB) {
50: PetscObjectGetName((PetscObject)vec,&name);
51: if (!rank && info) {
52: PetscFPrintf(PETSC_COMM_SELF,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
53: PetscFPrintf(PETSC_COMM_SELF,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
54: PetscFPrintf(PETSC_COMM_SELF,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
55: }
56: }
58: PetscObjectGetOptionsPrefix((PetscObject)vec,&pre);
59: if (!rank && info) {PetscFPrintf(PETSC_COMM_SELF,info,"-%svecload_block_size %D\n",pre?pre:"",PetscAbs(vec->map->bs));}
60: }
61: return(0);
62: }
64: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
65: {
67: PetscBool skipHeader,flg;
68: PetscInt tr[2],rows,N,n,s,bs;
69: PetscScalar *array;
70: PetscLayout map;
73: PetscViewerSetUp(viewer);
74: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
76: VecGetLayout(vec,&map);
77: PetscLayoutGetSize(map,&N);
79: /* read vector header */
80: if (!skipHeader) {
81: PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
82: if (tr[0] != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a vector next in file");
83: if (tr[1] < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Vector size (%D) in file is negative",tr[1]);
84: if (N >= 0 && N != tr[1]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Vector in file different size (%D) than input vector (%D)",tr[1],N);
85: rows = tr[1];
86: } else {
87: if (N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Vector binary file header was skipped, thus the user must specify the global size of input vector");
88: rows = N;
89: }
91: /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
92: PetscLayoutGetBlockSize(map,&bs);
93: PetscOptionsGetInt(((PetscObject)viewer)->options,((PetscObject)vec)->prefix,"-vecload_block_size",&bs,&flg);
94: if (flg) {VecSetBlockSize(vec,bs);}
95: PetscLayoutGetLocalSize(map,&n);
96: if (N < 0) {VecSetSizes(vec,n,rows);}
97: VecSetUp(vec);
99: /* get vector sizes and check global size */
100: VecGetSize(vec,&N);
101: VecGetLocalSize(vec,&n);
102: VecGetOwnershipRange(vec,&s,NULL);
103: if (N != rows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Vector in file different size (%D) than input vector (%D)",rows,N);
105: /* read vector values */
106: VecGetArray(vec,&array);
107: PetscViewerBinaryReadAll(viewer,array,n,s,N,PETSC_SCALAR);
108: VecRestoreArray(vec,&array);
109: return(0);
110: }
112: #if defined(PETSC_HAVE_HDF5)
113: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
114: {
115: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
116: PetscScalar *x,*arr;
117: const char *vecname;
121: if (!((PetscObject)xin)->name) SETERRQ(PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
122: #if defined(PETSC_USE_REAL_SINGLE)
123: scalartype = H5T_NATIVE_FLOAT;
124: #elif defined(PETSC_USE_REAL___FLOAT128)
125: #error "HDF5 output with 128 bit floats not supported."
126: #elif defined(PETSC_USE_REAL___FP16)
127: #error "HDF5 output with 16 bit floats not supported."
128: #else
129: scalartype = H5T_NATIVE_DOUBLE;
130: #endif
131: PetscObjectGetName((PetscObject)xin, &vecname);
132: PetscViewerHDF5Load(viewer,vecname,xin->map,scalartype,(void**)&x);
133: VecSetUp(xin); /* VecSetSizes might have not been called so ensure ops->create has been called */
134: if (!xin->ops->replacearray) {
135: VecGetArray(xin,&arr);
136: PetscArraycpy(arr,x,xin->map->n);
137: PetscFree(x);
138: VecRestoreArray(xin,&arr);
139: } else {
140: VecReplaceArray(xin,x);
141: }
142: return(0);
143: }
144: #endif
146: #if defined(PETSC_HAVE_ADIOS)
147: #include <adios.h>
148: #include <adios_read.h>
149: #include <petsc/private/vieweradiosimpl.h>
150: #include <petsc/private/viewerimpl.h>
152: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
153: {
154: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
155: PetscErrorCode ierr;
156: PetscScalar *x;
157: PetscInt Nfile,N,rstart,n;
158: uint64_t N_t,rstart_t;
159: const char *vecname;
160: ADIOS_VARINFO *v;
161: ADIOS_SELECTION *sel;
164: PetscObjectGetName((PetscObject) xin, &vecname);
166: v = adios_inq_var(adios->adios_fp, vecname);
167: if (v->ndim != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%D)", v->ndim);
168: Nfile = (PetscInt) v->dims[0];
170: /* Set Vec sizes,blocksize,and type if not already set */
171: if ((xin)->map->n < 0 && (xin)->map->N < 0) {
172: VecSetSizes(xin, PETSC_DECIDE, Nfile);
173: }
174: /* If sizes and type already set,check if the vector global size is correct */
175: VecGetSize(xin, &N);
176: VecGetLocalSize(xin, &n);
177: if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
178:
179: VecGetOwnershipRange(xin,&rstart,NULL);
180: rstart_t = rstart;
181: N_t = n;
182: sel = adios_selection_boundingbox (v->ndim, &rstart_t, &N_t);
183: VecGetArray(xin,&x);
184: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
185: adios_perform_reads (adios->adios_fp, 1);
186: VecRestoreArray(xin,&x);
187: adios_selection_delete(sel);
189: return(0);
190: }
191: #endif
193: #if defined(PETSC_HAVE_ADIOS2)
194: #include <adios2_c.h>
195: #include <petsc/private/vieweradios2impl.h>
196: #include <petsc/private/viewerimpl.h>
198: PetscErrorCode VecLoad_ADIOS2(Vec xin, PetscViewer viewer)
199: {
200: PetscViewer_ADIOS2 *adios2 = (PetscViewer_ADIOS2*)viewer->data;
201: PetscErrorCode ierr;
202: PetscScalar *x;
203: PetscInt Nfile,N,rstart,n;
204: const char *vecname;
207: PetscObjectGetName((PetscObject) xin, &vecname);
209: /* Set Vec sizes,blocksize,and type if not already set */
210: if ((xin)->map->n < 0 && (xin)->map->N < 0) {
211: VecSetSizes(xin, PETSC_DECIDE, Nfile);
212: }
213: /* If sizes and type already set,check if the vector global size is correct */
214: VecGetSize(xin, &N);
215: VecGetLocalSize(xin, &n);
216: if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
218: VecGetOwnershipRange(xin,&rstart,NULL);
219: VecGetArray(xin,&x);
220: VecRestoreArray(xin,&x);
221: return(0);
222: }
223: #endif
225: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
226: {
228: PetscBool isbinary;
229: #if defined(PETSC_HAVE_HDF5)
230: PetscBool ishdf5;
231: #endif
232: #if defined(PETSC_HAVE_ADIOS)
233: PetscBool isadios;
234: #endif
235: #if defined(PETSC_HAVE_ADIOS2)
236: PetscBool isadios2;
237: #endif
240: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
241: #if defined(PETSC_HAVE_HDF5)
242: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
243: #endif
244: #if defined(PETSC_HAVE_ADIOS)
245: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
246: #endif
247: #if defined(PETSC_HAVE_ADIOS2)
248: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS2,&isadios2);
249: #endif
251: #if defined(PETSC_HAVE_HDF5)
252: if (ishdf5) {
253: if (!((PetscObject)newvec)->name) {
254: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
255: 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()");
256: }
257: VecLoad_HDF5(newvec, viewer);
258: } else
259: #endif
260: #if defined(PETSC_HAVE_ADIOS)
261: if (isadios) {
262: if (!((PetscObject)newvec)->name) {
263: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
264: 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()");
265: }
266: VecLoad_ADIOS(newvec, viewer);
267: } else
268: #endif
269: #if defined(PETSC_HAVE_ADIOS2)
270: if (isadios2) {
271: if (!((PetscObject)newvec)->name) {
272: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
273: 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()");
274: }
275: VecLoad_ADIOS2(newvec, viewer);
276: } else
277: #endif
278: {
279: VecLoad_Binary(newvec, viewer);
280: }
281: return(0);
282: }
284: /*@
285: VecChop - Set all values in the vector with an absolute value less than the tolerance to zero
287: Input Parameters:
288: + v - The vector
289: - tol - The zero tolerance
291: Output Parameters:
292: . v - The chopped vector
294: Level: intermediate
296: .seealso: VecCreate(), VecSet()
297: @*/
298: PetscErrorCode VecChop(Vec v, PetscReal tol)
299: {
300: PetscScalar *a;
301: PetscInt n, i;
305: VecGetLocalSize(v, &n);
306: VecGetArray(v, &a);
307: for (i = 0; i < n; ++i) {
308: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
309: }
310: VecRestoreArray(v, &a);
311: return(0);
312: }