Actual source code: ex50.c

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

  3: /*
  4:   fftw vectors alloate 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:   PetscInitialize(&argc,&args,(char*)0,help);

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

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

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

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

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

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

 68: /*TEST

 70:    build:
 71:      requires: hdf5 fftw

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