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:   uint32_t     tr[2];
     67:   PetscInt     token, 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_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:   }

     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));

    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);

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

    118: #if defined(PETSC_HAVE_HDF5)
    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;

    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

    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>

    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;

    167:   PetscFunctionBegin;
    168:   PetscCall(PetscObjectGetName((PetscObject)xin, &vecname));

    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];

    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);

    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

    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

    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

    213: #if defined(PETSC_HAVE_HDF5)
    214:   if (ishdf5) {
    215:     if (!((PetscObject)newvec)->name) {
    216:       PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
    217:       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()");
    218:     }
    219:     PetscCall(VecLoad_HDF5(newvec, viewer));
    220:   } else
    221: #endif
    222: #if defined(PETSC_HAVE_ADIOS)
    223:     if (isadios) {
    224:     if (!((PetscObject)newvec)->name) {
    225:       PetscCall(PetscLogEventEnd(VEC_Load, viewer, 0, 0, 0));
    226:       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()");
    227:     }
    228:     PetscCall(VecLoad_ADIOS(newvec, viewer));
    229:   } else
    230: #endif
    231:   {
    232:     PetscCall(VecLoad_Binary(newvec, viewer));
    233:   }
    234:   PetscFunctionReturn(PETSC_SUCCESS);
    235: }

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

    240:   Input Parameters:
    241: + v   - The vector
    242: - tol - The zero tolerance

    244:   Output Parameter:
    245: . v - The filtered vector

    247:   Level: intermediate

    249: .seealso: `VecCreate()`, `VecSet()`, `MatFilter()`
    250: @*/
    251: PetscErrorCode VecFilter(Vec v, PetscReal tol)
    252: {
    253:   PetscScalar *a;
    254:   PetscInt     n;

    256:   PetscFunctionBegin;
    257:   PetscCall(VecGetLocalSize(v, &n));
    258:   PetscCall(VecGetArray(v, &a));
    259:   for (PetscInt i = 0; i < n; ++i) {
    260:     if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
    261:   }
    262:   PetscCall(VecRestoreArray(v, &a));
    263:   PetscFunctionReturn(PETSC_SUCCESS);
    264: }