Actual source code: ex50.c
petsc-3.13.4 2020-08-01
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*/