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:   PetscInt     tr[2], rows, N, n, s, bs;
 67:   PetscScalar *array;
 68:   PetscLayout  map;

 70:   PetscFunctionBegin;
 71:   PetscCall(PetscViewerSetUp(viewer));
 72:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));

 74:   PetscCall(VecGetLayout(vec, &map));
 75:   PetscCall(PetscLayoutGetSize(map, &N));

 77:   /* read vector header */
 78:   if (!skipHeader) {
 79:     PetscCall(PetscViewerBinaryRead(viewer, tr, 2, NULL, PETSC_INT));
 80:     PetscCheck(tr[0] == VEC_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a vector next in file");
 81:     PetscCheck(tr[1] >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Vector size (%" PetscInt_FMT ") in file is negative", tr[1]);
 82:     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);
 83:     rows = tr[1];
 84:   } else {
 85:     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");
 86:     rows = N;
 87:   }

 89:   /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
 90:   PetscCall(PetscLayoutGetBlockSize(map, &bs));
 91:   PetscCall(PetscOptionsGetInt(((PetscObject)viewer)->options, ((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flg));
 92:   if (flg) PetscCall(VecSetBlockSize(vec, bs));
 93:   PetscCall(PetscLayoutGetLocalSize(map, &n));
 94:   if (N < 0) PetscCall(VecSetSizes(vec, n, rows));
 95:   PetscCall(VecSetUp(vec));

 97:   /* get vector sizes and check global size */
 98:   PetscCall(VecGetSize(vec, &N));
 99:   PetscCall(VecGetLocalSize(vec, &n));
100:   PetscCall(VecGetOwnershipRange(vec, &s, NULL));
101:   PetscCheck(N == rows, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different size (%" PetscInt_FMT ") than input vector (%" PetscInt_FMT ")", rows, N);

103:   /* read vector values */
104:   PetscCall(VecGetArray(vec, &array));
105:   PetscCall(PetscViewerBinaryReadAll(viewer, array, n, s, N, PETSC_SCALAR));
106:   PetscCall(VecRestoreArray(vec, &array));
107:   PetscFunctionReturn(PETSC_SUCCESS);
108: }

110: #if defined(PETSC_HAVE_HDF5)
111: static PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
112: {
113:   hid_t        scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
114:   PetscScalar *x, *arr;
115:   const char  *vecname;

117:   PetscFunctionBegin;
118:   PetscCheck(((PetscObject)xin)->name, PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
119:   #if defined(PETSC_USE_REAL_SINGLE)
120:   scalartype = H5T_NATIVE_FLOAT;
121:   #elif defined(PETSC_USE_REAL___FLOAT128)
122:     #error "HDF5 output with 128 bit floats not supported."
123:   #elif defined(PETSC_USE_REAL___FP16)
124:     #error "HDF5 output with 16 bit floats not supported."
125:   #else
126:   scalartype = H5T_NATIVE_DOUBLE;
127:   #endif
128:   PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));
129:   PetscCall(PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void **)&x));
130:   PetscCall(VecSetUp(xin)); /* VecSetSizes might have not been called so ensure ops->create has been called */
131:   if (!xin->ops->replacearray) {
132:     PetscCall(VecGetArray(xin, &arr));
133:     PetscCall(PetscArraycpy(arr, x, xin->map->n));
134:     PetscCall(PetscFree(x));
135:     PetscCall(VecRestoreArray(xin, &arr));
136:   } else {
137:     PetscCall(VecReplaceArray(xin, x));
138:   }
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }
141: #endif

143: #if defined(PETSC_HAVE_ADIOS)
144:   #include <adios.h>
145:   #include <adios_read.h>
146: #include <petsc/private/vieweradiosimpl.h>
147: #include <petsc/private/viewerimpl.h>

149: static PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
150: {
151:   PetscViewer_ADIOS *adios = (PetscViewer_ADIOS *)viewer->data;
152:   PetscScalar       *x;
153:   PetscInt           Nfile, N, rstart, n;
154:   uint64_t           N_t, rstart_t;
155:   const char        *vecname;
156:   ADIOS_VARINFO     *v;
157:   ADIOS_SELECTION   *sel;

159:   PetscFunctionBegin;
160:   PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));

162:   v = adios_inq_var(adios->adios_fp, vecname);
163:   PetscCheck(v->ndim == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%" PetscInt_FMT ")", v->ndim);
164:   Nfile = (PetscInt)v->dims[0];

166:   /* Set Vec sizes,blocksize,and type if not already set */
167:   if ((xin)->map->n < 0 && (xin)->map->N < 0) PetscCall(VecSetSizes(xin, PETSC_DECIDE, Nfile));
168:   /* If sizes and type already set,check if the vector global size is correct */
169:   PetscCall(VecGetSize(xin, &N));
170:   PetscCall(VecGetLocalSize(xin, &n));
171:   PetscCheck(N == Nfile, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%" PetscInt_FMT ") then input vector (%" PetscInt_FMT ")", Nfile, N);

173:   PetscCall(VecGetOwnershipRange(xin, &rstart, NULL));
174:   rstart_t = rstart;
175:   N_t      = n;
176:   sel      = adios_selection_boundingbox(v->ndim, &rstart_t, &N_t);
177:   PetscCall(VecGetArray(xin, &x));
178:   adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
179:   adios_perform_reads(adios->adios_fp, 1);
180:   PetscCall(VecRestoreArray(xin, &x));
181:   adios_selection_delete(sel);

183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }
185: #endif

187: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
188: {
189:   PetscBool isbinary;
190: #if defined(PETSC_HAVE_HDF5)
191:   PetscBool ishdf5;
192: #endif
193: #if defined(PETSC_HAVE_ADIOS)
194:   PetscBool isadios;
195: #endif

197:   PetscFunctionBegin;
198:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
199: #if defined(PETSC_HAVE_HDF5)
200:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
201: #endif
202: #if defined(PETSC_HAVE_ADIOS)
203:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERADIOS, &isadios));
204: #endif

206: #if defined(PETSC_HAVE_HDF5)
207:   if (ishdf5) {
208:     if (!((PetscObject)newvec)->name) {
209:       PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
210:       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()");
211:     }
212:     PetscCall(VecLoad_HDF5(newvec, viewer));
213:   } else
214: #endif
215: #if defined(PETSC_HAVE_ADIOS)
216:     if (isadios) {
217:     if (!((PetscObject)newvec)->name) {
218:       PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
219:       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()");
220:     }
221:     PetscCall(VecLoad_ADIOS(newvec, viewer));
222:   } else
223: #endif
224:   {
225:     PetscCall(VecLoad_Binary(newvec, viewer));
226:   }
227:   PetscFunctionReturn(PETSC_SUCCESS);
228: }

230: /*@
231:   VecFilter - Set all values in the vector with an absolute value less than or equal to the tolerance to zero

233:   Input Parameters:
234: + v   - The vector
235: - tol - The zero tolerance

237:   Output Parameter:
238: . v - The filtered vector

240:   Level: intermediate

242: .seealso: `VecCreate()`, `VecSet()`, `MatFilter()`
243: @*/
244: PetscErrorCode VecFilter(Vec v, PetscReal tol)
245: {
246:   PetscScalar *a;
247:   PetscInt     n;

249:   PetscFunctionBegin;
250:   PetscCall(VecGetLocalSize(v, &n));
251:   PetscCall(VecGetArray(v, &a));
252:   for (PetscInt i = 0; i < n; ++i) {
253:     if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
254:   }
255:   PetscCall(VecRestoreArray(v, &a));
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }