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: }