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 */
6
7 #include <petscsys.h>
8 #include <petscvec.h> /*I "petscvec.h" I*/
9 #include <petsc/private/vecimpl.h>
10 #include <petsc/private/viewerimpl.h>
11 #include <petsclayouthdf5.h>
12
VecView_Binary(Vec vec,PetscViewer viewer)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;
19
20 PetscFunctionBegin;
21 PetscCheckSameComm(vec, 1, viewer, 2);
22 PetscCall(PetscViewerSetUp(viewer));
23 PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
24
25 PetscCall(VecGetLayout(vec, &map));
26 PetscCall(PetscLayoutGetLocalSize(map, &n));
27 PetscCall(PetscLayoutGetRange(map, &s, NULL));
28 PetscCall(PetscLayoutGetSize(map, &N));
29
30 tr[0] = VEC_FILE_CLASSID;
31 tr[1] = N;
32 if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, tr, 2, PETSC_INT));
33
34 PetscCall(VecGetArrayRead(vec, &xarray));
35 PetscCall(PetscViewerBinaryWriteAll(viewer, xarray, n, s, N, PETSC_SCALAR));
36 PetscCall(VecRestoreArrayRead(vec, &xarray));
37
38 { /* write to the viewer's .info file */
39 FILE *info;
40 PetscMPIInt rank;
41 PetscViewerFormat format;
42 const char *name, *pre;
43
44 PetscCall(PetscViewerBinaryGetInfoPointer(viewer, &info));
45 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)vec), &rank));
46
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 }
56
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 : "", vec->map->bs));
59 }
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
VecLoad_Binary(Vec vec,PetscViewer viewer)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;
70
71 PetscFunctionBegin;
72 PetscCall(PetscViewerSetUp(viewer));
73 PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
74
75 PetscCall(VecGetLayout(vec, &map));
76 PetscCall(PetscLayoutGetSize(map, &N));
77
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 rows = tr[1];
83 } else { // Assume file was written with 64-bit ints so reconstruct token and read number of rows
84 PetscInt64 rows64;
85 token = (PetscInt)((uint64_t)tr[0] << 32) + tr[1];
86 PetscCheck(token == VEC_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a vector next in file");
87 PetscCall(PetscViewerBinaryRead(viewer, &rows64, 1, NULL, PETSC_INT64));
88 PetscCall(PetscIntCast(rows64, &rows));
89 }
90 PetscCheck(rows >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Vector size (%" PetscInt_FMT ") in file is negative", rows);
91 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);
92 } else {
93 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");
94 rows = N;
95 }
96
97 /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
98 PetscCall(PetscLayoutGetBlockSize(map, &bs));
99 PetscCall(PetscOptionsGetInt(((PetscObject)viewer)->options, ((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flg));
100 if (flg) PetscCall(VecSetBlockSize(vec, bs));
101 PetscCall(PetscLayoutGetLocalSize(map, &n));
102 if (N < 0) PetscCall(VecSetSizes(vec, n, rows));
103 PetscCall(VecSetUp(vec));
104
105 /* get vector sizes and check global size */
106 PetscCall(VecGetSize(vec, &N));
107 PetscCall(VecGetLocalSize(vec, &n));
108 PetscCall(VecGetOwnershipRange(vec, &s, NULL));
109 PetscCheck(N == rows, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different size (%" PetscInt_FMT ") than input vector (%" PetscInt_FMT ")", rows, N);
110
111 /* read vector values */
112 PetscCall(VecGetArray(vec, &array));
113 PetscCall(PetscViewerBinaryReadAll(viewer, array, n, s, N, PETSC_SCALAR));
114 PetscCall(VecRestoreArray(vec, &array));
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 #if defined(PETSC_HAVE_HDF5)
VecLoad_HDF5(Vec xin,PetscViewer viewer)119 static PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
120 {
121 hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
122 PetscScalar *x, *arr;
123 const char *vecname;
124
125 PetscFunctionBegin;
126 PetscCheck(((PetscObject)xin)->name, PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
127 #if defined(PETSC_USE_REAL_SINGLE)
128 scalartype = H5T_NATIVE_FLOAT;
129 #elif defined(PETSC_USE_REAL___FLOAT128)
130 #error "HDF5 output with 128 bit floats not supported."
131 #elif defined(PETSC_USE_REAL___FP16)
132 #error "HDF5 output with 16 bit floats not supported."
133 #else
134 scalartype = H5T_NATIVE_DOUBLE;
135 #endif
136 PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
137 PetscCall(PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void **)&x));
138 PetscCall(VecSetUp(xin)); /* VecSetSizes might have not been called so ensure ops->create has been called */
139 if (!xin->ops->replacearray) {
140 PetscCall(VecGetArray(xin, &arr));
141 PetscCall(PetscArraycpy(arr, x, xin->map->n));
142 PetscCall(PetscFree(x));
143 PetscCall(VecRestoreArray(xin, &arr));
144 } else {
145 PetscCall(VecReplaceArray(xin, x));
146 }
147 PetscFunctionReturn(PETSC_SUCCESS);
148 }
149 #endif
150
151 #if defined(PETSC_HAVE_ADIOS)
152 #include <adios.h>
153 #include <adios_read.h>
154 #include <petsc/private/vieweradiosimpl.h>
155 #include <petsc/private/viewerimpl.h>
156
VecLoad_ADIOS(Vec xin,PetscViewer viewer)157 static PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
158 {
159 PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
160 PetscScalar *x;
161 PetscInt Nfile, N, rstart, n;
162 uint64_t N_t, rstart_t;
163 const char *vecname;
164 ADIOS_VARINFO *v;
165 ADIOS_SELECTION *sel;
166
167 PetscFunctionBegin;
168 PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
169
170 v = adios_inq_var(adios->adios_fp, vecname);
171 PetscCheck(v->ndim == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%" PetscInt_FMT ")", v->ndim);
172 Nfile = (PetscInt)v->dims[0];
173
174 /* Set Vec sizes,blocksize,and type if not already set */
175 if ((xin)->map->n < 0 && (xin)->map->N < 0) PetscCall(VecSetSizes(xin, PETSC_DECIDE, Nfile));
176 /* If sizes and type already set,check if the vector global size is correct */
177 PetscCall(VecGetSize(xin, &N));
178 PetscCall(VecGetLocalSize(xin, &n));
179 PetscCheck(N == Nfile, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%" PetscInt_FMT ") then input vector (%" PetscInt_FMT ")", Nfile, N);
180
181 PetscCall(VecGetOwnershipRange(xin, &rstart, NULL));
182 rstart_t = rstart;
183 N_t = n;
184 sel = adios_selection_boundingbox(v->ndim, &rstart_t, &N_t);
185 PetscCall(VecGetArray(xin, &x));
186 adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
187 adios_perform_reads(adios->adios_fp, 1);
188 PetscCall(VecRestoreArray(xin, &x));
189 adios_selection_delete(sel);
190 PetscFunctionReturn(PETSC_SUCCESS);
191 }
192 #endif
193
VecLoad_Default(Vec newvec,PetscViewer viewer)194 PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
195 {
196 PetscBool isbinary;
197 #if defined(PETSC_HAVE_HDF5)
198 PetscBool ishdf5;
199 #endif
200 #if defined(PETSC_HAVE_ADIOS)
201 PetscBool isadios;
202 #endif
203
204 PetscFunctionBegin;
205 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
206 #if defined(PETSC_HAVE_HDF5)
207 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
208 #endif
209 #if defined(PETSC_HAVE_ADIOS)
210 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
211 #endif
212
213 if (isbinary) {
214 PetscCall(VecLoad_Binary(newvec, viewer));
215 #if defined(PETSC_HAVE_HDF5)
216 } else if (ishdf5) {
217 if (!((PetscObject)newvec)->name) {
218 PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
219 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()");
220 }
221 PetscCall(VecLoad_HDF5(newvec, viewer));
222 #endif
223 #if defined(PETSC_HAVE_ADIOS)
224 } else 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 #endif
231 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
232 PetscFunctionReturn(PETSC_SUCCESS);
233 }
234
235 /*@
236 VecFilter - Set all values in the vector with an absolute value less than or equal to the tolerance to zero
237
238 Input Parameters:
239 + v - The vector
240 - tol - The zero tolerance
241
242 Output Parameter:
243 . v - The filtered vector
244
245 Level: intermediate
246
247 .seealso: `VecCreate()`, `VecSet()`, `MatFilter()`
248 @*/
VecFilter(Vec v,PetscReal tol)249 PetscErrorCode VecFilter(Vec v, PetscReal tol)
250 {
251 PetscScalar *a;
252 PetscInt n;
253
254 PetscFunctionBegin;
255 PetscCall(VecGetLocalSize(v, &n));
256 PetscCall(VecGetArray(v, &a));
257 for (PetscInt i = 0; i < n; ++i) {
258 if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
259 }
260 PetscCall(VecRestoreArray(v, &a));
261 PetscFunctionReturn(PETSC_SUCCESS);
262 }
263