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