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;

 25:   PetscInitialize(&argc, &args, (char *)0, help);

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

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

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

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

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

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

 69: /*TEST

 71:    build:
 72:      requires: hdf5 fftw

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