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