Actual source code: ex50.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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: {
 17:   PetscInt       i,low,high,ldim,iglobal;
 18:   PetscInt       m=64,dim[2]={8,8},DIM=2; /* FFT parameters */
 19:   Vec            u,u_,H;    /* wave, work and transfer function vectors */
 20:   Vec            slice_rid; /* vector to hold the refractive index */
 21:   Mat            A;         /* FFT-matrix to call FFTW via interface */
 22:   PetscViewer    viewer;    /* Load refractive index */
 23:   PetscScalar    v;

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

 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*/