Actual source code: ex50.c

  1: static char help[] = "Test if VecLoad_HDF5 can correctly handle FFTW vectors\n\n";

  3: /*
  4:   fftw vectors allocate their data array through fftw_malloc() and have their specialized VecDestroy().
  5:   When doing VecLoad on these vectors, we must take care of the v->array, v->array_allocated properly
  6:   to avoid memory leaks and double-free.

  8:   Contributed-by: Sajid Ali <sajidsyed2021@u.northwestern.edu>
  9: */

 11: #include <petscmat.h>
 12: #include <petscviewerhdf5.h>

 14: int main(int argc, char **args)
 15: {
 16:   PetscInt    i, low, high, ldim, iglobal;
 17:   PetscInt    m = 64, dim[2] = {8, 8}, DIM = 2; /* FFT parameters */
 18:   Vec         u, u_, H;                         /* wave, work and transfer function vectors */
 19:   Vec         slice_rid;                        /* vector to hold the refractive index */
 20:   Mat         A;                                /* FFT-matrix to call FFTW via interface */
 21:   PetscViewer viewer;                           /* Load refractive index */
 22:   PetscScalar v;

 24:   PetscFunctionBeginUser;
 25:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

 27:   /* Generate vector */
 28:   PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
 29:   PetscCall(PetscObjectSetName((PetscObject)u, "ref_index"));
 30:   PetscCall(VecSetSizes(u, PETSC_DECIDE, m));
 31:   PetscCall(VecSetFromOptions(u));
 32:   PetscCall(VecGetOwnershipRange(u, &low, &high));
 33:   PetscCall(VecGetLocalSize(u, &ldim));

 35:   for (i = 0; i < ldim; i++) {
 36:     iglobal = i + low;
 37:     v       = (PetscScalar)(i + low);
 38:     PetscCall(VecSetValues(u, 1, &iglobal, &v, INSERT_VALUES));
 39:   }
 40:   PetscCall(VecAssemblyBegin(u));
 41:   PetscCall(VecAssemblyEnd(u));
 42:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "ex50tmp.h5", FILE_MODE_WRITE, &viewer));
 43:   PetscCall(VecView(u, viewer));
 44:   PetscCall(VecDestroy(&u));
 45:   PetscCall(PetscViewerDestroy(&viewer));

 47:   /* Make FFT matrix (via interface) and create vecs aligned to it. */
 48:   PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));

 50:   /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */
 51:   PetscCall(MatCreateVecsFFTW(A, &u, &u_, &H));
 52:   PetscCall(VecDuplicate(u, &slice_rid));

 54:   /* Load refractive index from file */
 55:   PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "ex50tmp.h5", FILE_MODE_READ, &viewer));
 56:   PetscCall(PetscObjectSetName((PetscObject)slice_rid, "ref_index"));
 57:   PetscCall(VecLoad(slice_rid, viewer)); /* Test if VecLoad_HDF5 can correctly handle FFTW vectors */
 58:   PetscCall(PetscViewerDestroy(&viewer));

 60:   PetscCall(MatDestroy(&A));
 61:   PetscCall(VecDestroy(&slice_rid));
 62:   PetscCall(VecDestroy(&u));
 63:   PetscCall(VecDestroy(&u_));
 64:   PetscCall(VecDestroy(&H));
 65:   PetscCall(PetscFinalize());
 66:   return 0;
 67: }

 69: /*TEST

 71:    build:
 72:      requires: hdf5 fftw

 74:    test:
 75:      nsize: 2
 76:      requires: complex
 77: TEST*/