Actual source code: vecio.c
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: PetscBool skipHeader;
17: PetscLayout map;
18: PetscInt tr[2], n, s, N;
19: const PetscScalar *xarray;
21: PetscFunctionBegin;
22: PetscCheckSameComm(vec, 1, viewer, 2);
23: PetscCall(PetscViewerSetUp(viewer));
24: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
26: PetscCall(VecGetLayout(vec, &map));
27: PetscCall(PetscLayoutGetLocalSize(map, &n));
28: PetscCall(PetscLayoutGetRange(map, &s, NULL));
29: PetscCall(PetscLayoutGetSize(map, &N));
31: tr[0] = VEC_FILE_CLASSID;
32: tr[1] = N;
33: if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, tr, 2, PETSC_INT));
35: PetscCall(VecGetArrayRead(vec, &xarray));
36: PetscCall(PetscViewerBinaryWriteAll(viewer, xarray, n, s, N, PETSC_SCALAR));
37: PetscCall(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: PetscCall(PetscViewerBinaryGetInfoPointer(viewer, &info));
46: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)vec), &rank));
48: PetscCall(PetscViewerGetFormat(viewer, &format));
49: if (format == PETSC_VIEWER_BINARY_MATLAB) {
50: PetscCall(PetscObjectGetName((PetscObject)vec, &name));
51: if (rank == 0 && info) {
52: PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "#--- begin code written by PetscViewerBinary for MATLAB format ---#\n"));
53: PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "#$$ Set.%s = PetscBinaryRead(fd);\n", name));
54: PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n"));
55: }
56: }
58: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)vec, &pre));
59: if (rank == 0 && info) PetscCall(PetscFPrintf(PETSC_COMM_SELF, info, "-%svecload_block_size %" PetscInt_FMT "\n", pre ? pre : "", PetscAbs(vec->map->bs)));
60: }
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
65: {
66: PetscBool skipHeader, flg;
67: PetscInt tr[2], 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_INT));
81: PetscCheck(tr[0] == VEC_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a vector next in file");
82: PetscCheck(tr[1] >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Vector size (%" PetscInt_FMT ") in file is negative", tr[1]);
83: PetscCheck(N < 0 || N == tr[1], PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different size (%" PetscInt_FMT ") than input vector (%" PetscInt_FMT ")", tr[1], N);
84: rows = tr[1];
85: } else {
86: 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");
87: rows = N;
88: }
90: /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
91: PetscCall(PetscLayoutGetBlockSize(map, &bs));
92: PetscCall(PetscOptionsGetInt(((PetscObject)viewer)->options, ((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flg));
93: if (flg) PetscCall(VecSetBlockSize(vec, bs));
94: PetscCall(PetscLayoutGetLocalSize(map, &n));
95: if (N < 0) PetscCall(VecSetSizes(vec, n, rows));
96: PetscCall(VecSetUp(vec));
98: /* get vector sizes and check global size */
99: PetscCall(VecGetSize(vec, &N));
100: PetscCall(VecGetLocalSize(vec, &n));
101: PetscCall(VecGetOwnershipRange(vec, &s, NULL));
102: PetscCheck(N == rows, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different size (%" PetscInt_FMT ") than input vector (%" PetscInt_FMT ")", rows, N);
104: /* read vector values */
105: PetscCall(VecGetArray(vec, &array));
106: PetscCall(PetscViewerBinaryReadAll(viewer, array, n, s, N, PETSC_SCALAR));
107: PetscCall(VecRestoreArray(vec, &array));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: #if defined(PETSC_HAVE_HDF5)
112: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
113: {
114: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
115: PetscScalar *x, *arr;
116: const char *vecname;
118: PetscFunctionBegin;
119: PetscCheck(((PetscObject)xin)->name, PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
120: #if defined(PETSC_USE_REAL_SINGLE)
121: scalartype = H5T_NATIVE_FLOAT;
122: #elif defined(PETSC_USE_REAL___FLOAT128)
123: #error "HDF5 output with 128 bit floats not supported."
124: #elif defined(PETSC_USE_REAL___FP16)
125: #error "HDF5 output with 16 bit floats not supported."
126: #else
127: scalartype = H5T_NATIVE_DOUBLE;
128: #endif
129: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
130: PetscCall(PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void **)&x));
131: PetscCall(VecSetUp(xin)); /* VecSetSizes might have not been called so ensure ops->create has been called */
132: if (!xin->ops->replacearray) {
133: PetscCall(VecGetArray(xin, &arr));
134: PetscCall(PetscArraycpy(arr, x, xin->map->n));
135: PetscCall(PetscFree(x));
136: PetscCall(VecRestoreArray(xin, &arr));
137: } else {
138: PetscCall(VecReplaceArray(xin, x));
139: }
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
142: #endif
144: #if defined(PETSC_HAVE_ADIOS)
145: #include <adios.h>
146: #include <adios_read.h>
147: #include <petsc/private/vieweradiosimpl.h>
148: #include <petsc/private/viewerimpl.h>
150: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
151: {
152: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
153: PetscScalar *x;
154: PetscInt Nfile, N, rstart, n;
155: uint64_t N_t, rstart_t;
156: const char *vecname;
157: ADIOS_VARINFO *v;
158: ADIOS_SELECTION *sel;
160: PetscFunctionBegin;
161: PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
163: v = adios_inq_var(adios->adios_fp, vecname);
164: PetscCheck(v->ndim == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%" PetscInt_FMT ")", v->ndim);
165: Nfile = (PetscInt)v->dims[0];
167: /* Set Vec sizes,blocksize,and type if not already set */
168: if ((xin)->map->n < 0 && (xin)->map->N < 0) PetscCall(VecSetSizes(xin, PETSC_DECIDE, Nfile));
169: /* If sizes and type already set,check if the vector global size is correct */
170: PetscCall(VecGetSize(xin, &N));
171: PetscCall(VecGetLocalSize(xin, &n));
172: PetscCheck(N == Nfile, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%" PetscInt_FMT ") then input vector (%" PetscInt_FMT ")", Nfile, N);
174: PetscCall(VecGetOwnershipRange(xin, &rstart, NULL));
175: rstart_t = rstart;
176: N_t = n;
177: sel = adios_selection_boundingbox(v->ndim, &rstart_t, &N_t);
178: PetscCall(VecGetArray(xin, &x));
179: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
180: adios_perform_reads(adios->adios_fp, 1);
181: PetscCall(VecRestoreArray(xin, &x));
182: adios_selection_delete(sel);
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
186: #endif
188: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
189: {
190: PetscBool isbinary;
191: #if defined(PETSC_HAVE_HDF5)
192: PetscBool ishdf5;
193: #endif
194: #if defined(PETSC_HAVE_ADIOS)
195: PetscBool isadios;
196: #endif
198: PetscFunctionBegin;
199: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
200: #if defined(PETSC_HAVE_HDF5)
201: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
202: #endif
203: #if defined(PETSC_HAVE_ADIOS)
204: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
205: #endif
207: #if defined(PETSC_HAVE_HDF5)
208: if (ishdf5) {
209: if (!((PetscObject)newvec)->name) {
210: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
211: 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()");
212: }
213: PetscCall(VecLoad_HDF5(newvec, viewer));
214: } else
215: #endif
216: #if defined(PETSC_HAVE_ADIOS)
217: if (isadios) {
218: if (!((PetscObject)newvec)->name) {
219: PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
220: 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()");
221: }
222: PetscCall(VecLoad_ADIOS(newvec, viewer));
223: } else
224: #endif
225: {
226: PetscCall(VecLoad_Binary(newvec, viewer));
227: }
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: VecChop - Set all values in the vector with an absolute value less than the tolerance to zero
234: Input Parameters:
235: + v - The vector
236: - tol - The zero tolerance
238: Output Parameter:
239: . v - The chopped vector
241: Level: intermediate
243: .seealso: `VecCreate()`, `VecSet()`
244: @*/
245: PetscErrorCode VecChop(Vec v, PetscReal tol)
246: {
247: PetscScalar *a;
248: PetscInt n, i;
250: PetscFunctionBegin;
251: PetscCall(VecGetLocalSize(v, &n));
252: PetscCall(VecGetArray(v, &a));
253: for (i = 0; i < n; ++i) {
254: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
255: }
256: PetscCall(VecRestoreArray(v, &a));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }