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;
22: PetscViewerSetUp(viewer);
23: PetscViewerBinaryGetSkipHeader(viewer, &skipHeader);
25: VecGetLayout(vec, &map);
26: PetscLayoutGetLocalSize(map, &n);
27: PetscLayoutGetRange(map, &s, NULL);
28: PetscLayoutGetSize(map, &N);
30: tr[0] = VEC_FILE_CLASSID;
31: tr[1] = N;
32: if (!skipHeader) PetscViewerBinaryWrite(viewer, tr, 2, PETSC_INT);
34: VecGetArrayRead(vec, &xarray);
35: PetscViewerBinaryWriteAll(viewer, xarray, n, s, N, PETSC_SCALAR);
36: 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: PetscViewerBinaryGetInfoPointer(viewer, &info);
45: MPI_Comm_rank(PetscObjectComm((PetscObject)vec), &rank);
47: PetscViewerGetFormat(viewer, &format);
48: if (format == PETSC_VIEWER_BINARY_MATLAB) {
49: PetscObjectGetName((PetscObject)vec, &name);
50: if (rank == 0 && info) {
51: PetscFPrintf(PETSC_COMM_SELF, info, "#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
52: PetscFPrintf(PETSC_COMM_SELF, info, "#$$ Set.%s = PetscBinaryRead(fd);\n", name);
53: PetscFPrintf(PETSC_COMM_SELF, info, "#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
54: }
55: }
57: PetscObjectGetOptionsPrefix((PetscObject)vec, &pre);
58: if (rank == 0 && info) PetscFPrintf(PETSC_COMM_SELF, info, "-%svecload_block_size %" PetscInt_FMT "\n", pre ? pre : "", PetscAbs(vec->map->bs));
59: }
60: return 0;
61: }
63: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
64: {
65: PetscBool skipHeader, flg;
66: PetscInt tr[2], rows, N, n, s, bs;
67: PetscScalar *array;
68: PetscLayout map;
70: PetscViewerSetUp(viewer);
71: PetscViewerBinaryGetSkipHeader(viewer, &skipHeader);
73: VecGetLayout(vec, &map);
74: PetscLayoutGetSize(map, &N);
76: /* read vector header */
77: if (!skipHeader) {
78: PetscViewerBinaryRead(viewer, tr, 2, NULL, PETSC_INT);
82: rows = tr[1];
83: } else {
85: rows = N;
86: }
88: /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
89: PetscLayoutGetBlockSize(map, &bs);
90: PetscOptionsGetInt(((PetscObject)viewer)->options, ((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flg);
91: if (flg) VecSetBlockSize(vec, bs);
92: PetscLayoutGetLocalSize(map, &n);
93: if (N < 0) VecSetSizes(vec, n, rows);
94: VecSetUp(vec);
96: /* get vector sizes and check global size */
97: VecGetSize(vec, &N);
98: VecGetLocalSize(vec, &n);
99: VecGetOwnershipRange(vec, &s, NULL);
102: /* read vector values */
103: VecGetArray(vec, &array);
104: PetscViewerBinaryReadAll(viewer, array, n, s, N, PETSC_SCALAR);
105: VecRestoreArray(vec, &array);
106: return 0;
107: }
109: #if defined(PETSC_HAVE_HDF5)
110: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
111: {
112: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
113: PetscScalar *x, *arr;
114: const char *vecname;
117: #if defined(PETSC_USE_REAL_SINGLE)
118: scalartype = H5T_NATIVE_FLOAT;
119: #elif defined(PETSC_USE_REAL___FLOAT128)
120: #error "HDF5 output with 128 bit floats not supported."
121: #elif defined(PETSC_USE_REAL___FP16)
122: #error "HDF5 output with 16 bit floats not supported."
123: #else
124: scalartype = H5T_NATIVE_DOUBLE;
125: #endif
126: PetscObjectGetName((PetscObject)xin, &vecname);
127: PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void **)&x);
128: VecSetUp(xin); /* VecSetSizes might have not been called so ensure ops->create has been called */
129: if (!xin->ops->replacearray) {
130: VecGetArray(xin, &arr);
131: PetscArraycpy(arr, x, xin->map->n);
132: PetscFree(x);
133: VecRestoreArray(xin, &arr);
134: } else {
135: VecReplaceArray(xin, x);
136: }
137: return 0;
138: }
139: #endif
141: #if defined(PETSC_HAVE_ADIOS)
142: #include <adios.h>
143: #include <adios_read.h>
144: #include <petsc/private/vieweradiosimpl.h>
145: #include <petsc/private/viewerimpl.h>
147: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
148: {
149: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
150: PetscScalar *x;
151: PetscInt Nfile, N, rstart, n;
152: uint64_t N_t, rstart_t;
153: const char *vecname;
154: ADIOS_VARINFO *v;
155: ADIOS_SELECTION *sel;
157: PetscObjectGetName((PetscObject)xin, &vecname);
159: v = adios_inq_var(adios->adios_fp, vecname);
161: Nfile = (PetscInt)v->dims[0];
163: /* Set Vec sizes,blocksize,and type if not already set */
164: if ((xin)->map->n < 0 && (xin)->map->N < 0) VecSetSizes(xin, PETSC_DECIDE, Nfile);
165: /* If sizes and type already set,check if the vector global size is correct */
166: VecGetSize(xin, &N);
167: VecGetLocalSize(xin, &n);
170: VecGetOwnershipRange(xin, &rstart, NULL);
171: rstart_t = rstart;
172: N_t = n;
173: sel = adios_selection_boundingbox(v->ndim, &rstart_t, &N_t);
174: VecGetArray(xin, &x);
175: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
176: adios_perform_reads(adios->adios_fp, 1);
177: VecRestoreArray(xin, &x);
178: adios_selection_delete(sel);
180: return 0;
181: }
182: #endif
184: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
185: {
186: PetscBool isbinary;
187: #if defined(PETSC_HAVE_HDF5)
188: PetscBool ishdf5;
189: #endif
190: #if defined(PETSC_HAVE_ADIOS)
191: PetscBool isadios;
192: #endif
194: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
195: #if defined(PETSC_HAVE_HDF5)
196: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
197: #endif
198: #if defined(PETSC_HAVE_ADIOS)
199: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios);
200: #endif
202: #if defined(PETSC_HAVE_HDF5)
203: if (ishdf5) {
204: if (!((PetscObject)newvec)->name) {
205: PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0);
206: 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()");
207: }
208: VecLoad_HDF5(newvec, viewer);
209: } else
210: #endif
211: #if defined(PETSC_HAVE_ADIOS)
212: if (isadios) {
213: if (!((PetscObject)newvec)->name) {
214: PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0);
215: 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()");
216: }
217: VecLoad_ADIOS(newvec, viewer);
218: } else
219: #endif
220: {
221: VecLoad_Binary(newvec, viewer);
222: }
223: return 0;
224: }
226: /*@
227: VecChop - Set all values in the vector with an absolute value less than the tolerance to zero
229: Input Parameters:
230: + v - The vector
231: - tol - The zero tolerance
233: Output Parameters:
234: . v - The chopped vector
236: Level: intermediate
238: .seealso: `VecCreate()`, `VecSet()`
239: @*/
240: PetscErrorCode VecChop(Vec v, PetscReal tol)
241: {
242: PetscScalar *a;
243: PetscInt n, i;
245: VecGetLocalSize(v, &n);
246: VecGetArray(v, &a);
247: for (i = 0; i < n; ++i) {
248: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
249: }
250: VecRestoreArray(v, &a);
251: return 0;
252: }