Actual source code: vecio.c
1: /*
2: This file contains simple binary input routines for vectors. The
3: analogous output routines are within each vector implementation's
4: VecView (with viewer types PETSCVIEWERBINARY)
5: */
7: #include <petscsys.h>
8: #include <petscvec.h>
9: #include <petsc/private/vecimpl.h>
10: #include <petsc/private/viewerimpl.h>
11: #include <petsclayouthdf5.h>
13: PetscErrorCode VecView_Binary(Vec vec, PetscViewer viewer)
14: {
15: PetscBool skipHeader;
16: PetscLayout map;
17: PetscInt tr[2], n, s, N;
18: const PetscScalar *xarray;
20: PetscFunctionBegin;
21: PetscCheckSameComm(vec, 1, viewer, 2);
22: PetscCall(PetscViewerSetUp(viewer));
23: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
25: PetscCall(VecGetLayout(vec, &map));
26: PetscCall(PetscLayoutGetLocalSize(map, &n));
27: PetscCall(PetscLayoutGetRange(map, &s, NULL));
28: PetscCall(PetscLayoutGetSize(map, &N));
30: tr[0] = VEC_FILE_CLASSID;
31: tr[1] = N;
32: if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, tr, 2, PETSC_INT));
34: PetscCall(VecGetArrayRead(vec, &xarray));
35: PetscCall(PetscViewerBinaryWriteAll(viewer, xarray, n, s, N, PETSC_SCALAR));
36: PetscCall(VecRestoreArrayRead(vec, &xarray));
38: { /* write to the viewer's .info file */
39: FILE *info;
40: PetscMPIInt rank;
41: PetscViewerFormat format;
42: const char *name, *pre;
44: PetscCall(PetscViewerBinaryGetInfoPointer(viewer, &info));
45: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)vec), &rank));
47: PetscCall(PetscViewerGetFormat(viewer, &format));
48: if (format == PETSC_VIEWER_BINARY_MATLAB) {
49: PetscCall(PetscObjectGetName((PetscObject)vec, &name));
50: if (rank == 0 && info) {
51: PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "#--- begin code written by PetscViewerBinary for MATLAB format ---#\n"));
52: PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "#$$ Set.%s = PetscBinaryRead(fd);\n", name));
53: PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n"));
54: }
55: }
57: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)vec, &pre));
58: if (rank == 0 && info) PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "-%svecload_block_size %" PetscInt_FMT "\n", pre ? pre : "", PetscAbs(vec->map->bs)));
59: }
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
64: {
65: PetscBool skipHeader, flg;
66: uint32_t tr[2];
67: PetscInt token, rows, N, n, s, bs;
68: PetscScalar *array;
69: PetscLayout map;
71: PetscFunctionBegin;
72: PetscCall(PetscViewerSetUp(viewer));
73: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
75: PetscCall(VecGetLayout(vec, &map));
76: PetscCall(PetscLayoutGetSize(map, &N));
78: /* read vector header */
79: if (!skipHeader) {
80: PetscCall(PetscViewerBinaryRead(viewer, tr, 2, NULL, PETSC_INT32));
81: if (tr[0] == VEC_FILE_CLASSID) { // File was written with 32-bit ints
82: token = tr[0];
83: rows = tr[1];
84: } else { // Assume file was written with 64-bit ints so reconstruct token and read number of rows
85: PetscInt64 rows64;
86: token = ((uint64_t)tr[0] << 32) + tr[1];
87: PetscCheck(token == VEC_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a vector next in file");
88: PetscCall(PetscViewerBinaryRead(viewer, &rows64, 1, NULL, PETSC_INT64));
89: PetscCall(PetscIntCast(rows64, &rows));
90: }
91: PetscCheck(rows >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Vector size (%" PetscInt_FMT ") in file is negative", rows);
92: PetscCheck(N < 0 || N == rows, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different size (%" PetscInt_FMT ") than input vector (%" PetscInt_FMT ")", rows, N);
93: } else {
94: PetscCheck(N >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Vector binary file header was skipped, thus the user must specify the global size of input vector");
95: rows = N;
96: }
98: /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
99: PetscCall(PetscLayoutGetBlockSize(map, &bs));
100: PetscCall(PetscOptionsGetInt(((PetscObject)viewer)->options, ((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flg));
101: if (flg) PetscCall(VecSetBlockSize(vec, bs));
102: PetscCall(PetscLayoutGetLocalSize(map, &n));
103: if (N < 0) PetscCall(VecSetSizes(vec, n, rows));
104: PetscCall(VecSetUp(vec));
106: /* get vector sizes and check global size */
107: PetscCall(VecGetSize(vec, &N));
108: PetscCall(VecGetLocalSize(vec, &n));
109: PetscCall(VecGetOwnershipRange(vec, &s, NULL));
110: PetscCheck(N == rows, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different size (%" PetscInt_FMT ") than input vector (%" PetscInt_FMT ")", rows, N);
112: /* read vector values */
113: PetscCall(VecGetArray(vec, &array));
114: PetscCall(PetscViewerBinaryReadAll(viewer, array, n, s, N, PETSC_SCALAR));
115: PetscCall(VecRestoreArray(vec, &array));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: #if defined(PETSC_HAVE_HDF5)
120: static PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
121: {
122: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
123: PetscScalar *x, *arr;
124: const char *vecname;
126: PetscFunctionBegin;
127: PetscCheck(((PetscObject)xin)->name, PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
128: #if defined(PETSC_USE_REAL_SINGLE)
129: scalartype = H5T_NATIVE_FLOAT;
130: #elif defined(PETSC_USE_REAL___FLOAT128)
131: #error "HDF5 output with 128 bit floats not supported."
132: #elif defined(PETSC_USE_REAL___FP16)
133: #error "HDF5 output with 16 bit floats not supported."
134: #else
135: scalartype = H5T_NATIVE_DOUBLE;
136: #endif
137: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
138: PetscCall(PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void **)&x));
139: PetscCall(VecSetUp(xin)); /* VecSetSizes might have not been called so ensure ops->create has been called */
140: if (!xin->ops->replacearray) {
141: PetscCall(VecGetArray(xin, &arr));
142: PetscCall(PetscArraycpy(arr, x, xin->map->n));
143: PetscCall(PetscFree(x));
144: PetscCall(VecRestoreArray(xin, &arr));
145: } else {
146: PetscCall(VecReplaceArray(xin, x));
147: }
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
150: #endif
152: #if defined(PETSC_HAVE_ADIOS)
153: #include <adios.h>
154: #include <adios_read.h>
155: #include <petsc/private/vieweradiosimpl.h>
156: #include <petsc/private/viewerimpl.h>
158: static PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
159: {
160: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
161: PetscScalar *x;
162: PetscInt Nfile, N, rstart, n;
163: uint64_t N_t, rstart_t;
164: const char *vecname;
165: ADIOS_VARINFO *v;
166: ADIOS_SELECTION *sel;
168: PetscFunctionBegin;
169: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
171: v = adios_inq_var(adios->adios_fp, vecname);
172: PetscCheck(v->ndim == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%" PetscInt_FMT ")", v->ndim);
173: Nfile = (PetscInt)v->dims[0];
175: /* Set Vec sizes,blocksize,and type if not already set */
176: if ((xin)->map->n < 0 && (xin)->map->N < 0) PetscCall(VecSetSizes(xin, PETSC_DECIDE, Nfile));
177: /* If sizes and type already set,check if the vector global size is correct */
178: PetscCall(VecGetSize(xin, &N));
179: PetscCall(VecGetLocalSize(xin, &n));
180: PetscCheck(N == Nfile, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%" PetscInt_FMT ") then input vector (%" PetscInt_FMT ")", Nfile, N);
182: PetscCall(VecGetOwnershipRange(xin, &rstart, NULL));
183: rstart_t = rstart;
184: N_t = n;
185: sel = adios_selection_boundingbox(v->ndim, &rstart_t, &N_t);
186: PetscCall(VecGetArray(xin, &x));
187: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
188: adios_perform_reads(adios->adios_fp, 1);
189: PetscCall(VecRestoreArray(xin, &x));
190: adios_selection_delete(sel);
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
193: #endif
195: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
196: {
197: PetscBool isbinary;
198: #if defined(PETSC_HAVE_HDF5)
199: PetscBool ishdf5;
200: #endif
201: #if defined(PETSC_HAVE_ADIOS)
202: PetscBool isadios;
203: #endif
205: PetscFunctionBegin;
206: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
207: #if defined(PETSC_HAVE_HDF5)
208: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
209: #endif
210: #if defined(PETSC_HAVE_ADIOS)
211: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
212: #endif
214: #if defined(PETSC_HAVE_HDF5)
215: if (ishdf5) {
216: if (!((PetscObject)newvec)->name) {
217: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
218: 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()");
219: }
220: PetscCall(VecLoad_HDF5(newvec, viewer));
221: } else
222: #endif
223: #if defined(PETSC_HAVE_ADIOS)
224: if (isadios) {
225: if (!((PetscObject)newvec)->name) {
226: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
227: 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()");
228: }
229: PetscCall(VecLoad_ADIOS(newvec, viewer));
230: } else
231: #endif
232: {
233: PetscCall(VecLoad_Binary(newvec, viewer));
234: }
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: VecFilter - Set all values in the vector with an absolute value less than or equal to the tolerance to zero
241: Input Parameters:
242: + v - The vector
243: - tol - The zero tolerance
245: Output Parameter:
246: . v - The filtered vector
248: Level: intermediate
250: .seealso: `VecCreate()`, `VecSet()`, `MatFilter()`
251: @*/
252: PetscErrorCode VecFilter(Vec v, PetscReal tol)
253: {
254: PetscScalar *a;
255: PetscInt n;
257: PetscFunctionBegin;
258: PetscCall(VecGetLocalSize(v, &n));
259: PetscCall(VecGetArray(v, &a));
260: for (PetscInt i = 0; i < n; ++i) {
261: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
262: }
263: PetscCall(VecRestoreArray(v, &a));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }