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;
24: PetscFunctionBeginUser;
25: PetscCall(PetscInitialize(&argc, &args, NULL, help));
27: /* Generate vector */
28: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
29: PetscCall(PetscObjectSetName((PetscObject)u, "ref_index"));
30: PetscCall(VecSetSizes(u, PETSC_DECIDE, m));
31: PetscCall(VecSetFromOptions(u));
32: PetscCall(VecGetOwnershipRange(u, &low, &high));
33: PetscCall(VecGetLocalSize(u, &ldim));
35: for (i = 0; i < ldim; i++) {
36: iglobal = i + low;
37: v = (PetscScalar)(i + low);
38: PetscCall(VecSetValues(u, 1, &iglobal, &v, INSERT_VALUES));
39: }
40: PetscCall(VecAssemblyBegin(u));
41: PetscCall(VecAssemblyEnd(u));
42: PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "ex50tmp.h5", FILE_MODE_WRITE, &viewer));
43: PetscCall(VecView(u, viewer));
44: PetscCall(VecDestroy(&u));
45: PetscCall(PetscViewerDestroy(&viewer));
47: /* Make FFT matrix (via interface) and create vecs aligned to it. */
48: PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));
50: /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */
51: PetscCall(MatCreateVecsFFTW(A, &u, &u_, &H));
52: PetscCall(VecDuplicate(u, &slice_rid));
54: /* Load refractive index from file */
55: PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "ex50tmp.h5", FILE_MODE_READ, &viewer));
56: PetscCall(PetscObjectSetName((PetscObject)slice_rid, "ref_index"));
57: PetscCall(VecLoad(slice_rid, viewer)); /* Test if VecLoad_HDF5 can correctly handle FFTW vectors */
58: PetscCall(PetscViewerDestroy(&viewer));
60: PetscCall(MatDestroy(&A));
61: PetscCall(VecDestroy(&slice_rid));
62: PetscCall(VecDestroy(&u));
63: PetscCall(VecDestroy(&u_));
64: PetscCall(VecDestroy(&H));
65: PetscCall(PetscFinalize());
66: return 0;
67: }
69: /*TEST
71: build:
72: requires: hdf5 fftw
74: test:
75: nsize: 2
76: requires: complex
77: TEST*/