1c4762a1bSJed Brown static char help[] = "Test duplication/destruction of FFTW vecs \n\n";
2c4762a1bSJed Brown
3c4762a1bSJed Brown /*
4c4762a1bSJed Brown Compiling the code:
5c4762a1bSJed Brown This code uses the FFTW interface.
6c4762a1bSJed Brown Use one of the options below to configure:
7c4762a1bSJed Brown --with-fftw-dir=/.... or --download-fftw
8c4762a1bSJed Brown Usage:
9c4762a1bSJed Brown mpiexec -np <np> ./ex228
10c4762a1bSJed Brown */
11c4762a1bSJed Brown
12c4762a1bSJed Brown #include <petscmat.h>
main(int argc,char ** args)13d71ae5a4SJacob Faibussowitsch int main(int argc, char **args)
14d71ae5a4SJacob Faibussowitsch {
15c4762a1bSJed Brown Mat A; /* FFT Matrix */
16c4762a1bSJed Brown Vec x, y, z; /* Work vectors */
17c4762a1bSJed Brown Vec x1, y1, z1; /* Duplicate vectors */
18c4762a1bSJed Brown PetscRandom rdm; /* for creating random input */
19c4762a1bSJed Brown PetscScalar a; /* used to scale output */
20c4762a1bSJed Brown PetscReal enorm; /* norm for sanity check */
21c4762a1bSJed Brown PetscInt n = 10, N = 1; /* FFT dimension params */
22c4762a1bSJed Brown PetscInt DIM, dim[5]; /* FFT params */
23c4762a1bSJed Brown
24327415f7SBarry Smith PetscFunctionBeginUser;
25c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &args, NULL, help));
269566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
27c4762a1bSJed Brown
28c4762a1bSJed Brown /* To create random input vector */
299566063dSJacob Faibussowitsch PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
309566063dSJacob Faibussowitsch PetscCall(PetscRandomSetFromOptions(rdm));
31c4762a1bSJed Brown
32c4762a1bSJed Brown /* Iterate over dimensions, use PETSc-FFTW interface */
33*e8bdd179SBarry Smith for (PetscInt i = 1; i < 5; i++) {
34c4762a1bSJed Brown DIM = i;
35c4762a1bSJed Brown N = 1;
36*e8bdd179SBarry Smith for (PetscInt k = 0; k < i; k++) {
379371c9d4SSatish Balay dim[k] = n;
389371c9d4SSatish Balay N *= n;
399371c9d4SSatish Balay }
40c4762a1bSJed Brown
419566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N));
42c4762a1bSJed Brown
43c4762a1bSJed Brown /* create FFTW object */
449566063dSJacob Faibussowitsch PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A));
45c4762a1bSJed Brown /* create vectors of length N */
469566063dSJacob Faibussowitsch PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
47c4762a1bSJed Brown
489566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
499566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
509566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
51c4762a1bSJed Brown
52c4762a1bSJed Brown /* Test vector duplication*/
539566063dSJacob Faibussowitsch PetscCall(VecDuplicate(x, &x1));
549566063dSJacob Faibussowitsch PetscCall(VecDuplicate(y, &y1));
559566063dSJacob Faibussowitsch PetscCall(VecDuplicate(z, &z1));
56c4762a1bSJed Brown
57c4762a1bSJed Brown /* Set values of space vector x, copy to duplicate */
589566063dSJacob Faibussowitsch PetscCall(VecSetRandom(x, rdm));
599566063dSJacob Faibussowitsch PetscCall(VecCopy(x, x1));
60c4762a1bSJed Brown
61c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD */
629566063dSJacob Faibussowitsch PetscCall(MatMult(A, x, y));
639566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(A, y, z));
64c4762a1bSJed Brown
65c4762a1bSJed Brown /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
669566063dSJacob Faibussowitsch PetscCall(MatMult(A, x1, y1));
679566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(A, y1, z1));
68c4762a1bSJed Brown
69c4762a1bSJed Brown /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
70c4762a1bSJed Brown a = 1.0 / (PetscReal)N;
719566063dSJacob Faibussowitsch PetscCall(VecScale(z1, a));
729566063dSJacob Faibussowitsch PetscCall(VecAXPY(z1, -1.0, x));
739566063dSJacob Faibussowitsch PetscCall(VecNorm(z1, NORM_1, &enorm));
74*e8bdd179SBarry Smith if (enorm > 1.e-9) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z1| %g dimension %" PetscInt_FMT "\n", enorm, i));
75c4762a1bSJed Brown
76c4762a1bSJed Brown /* free spaces */
779566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x1));
789566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y1));
799566063dSJacob Faibussowitsch PetscCall(VecDestroy(&z1));
809566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x));
819566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y));
829566063dSJacob Faibussowitsch PetscCall(VecDestroy(&z));
839566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A));
84c4762a1bSJed Brown }
85c4762a1bSJed Brown
869566063dSJacob Faibussowitsch PetscCall(PetscRandomDestroy(&rdm));
879566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
88b122ec5aSJacob Faibussowitsch return 0;
89c4762a1bSJed Brown }
90c4762a1bSJed Brown
91c4762a1bSJed Brown /*TEST
92c4762a1bSJed Brown
93c4762a1bSJed Brown build:
94c4762a1bSJed Brown requires: fftw complex
95c4762a1bSJed Brown
96c4762a1bSJed Brown test:
97c4762a1bSJed Brown suffix: 2
98c4762a1bSJed Brown nsize : 4
99c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16
100c4762a1bSJed Brown
101c4762a1bSJed Brown test:
102c4762a1bSJed Brown suffix: 3
103c4762a1bSJed Brown nsize : 2
104c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_MEASURE -n 12
105c4762a1bSJed Brown
106c4762a1bSJed Brown test:
107c4762a1bSJed Brown suffix: 4
108c4762a1bSJed Brown nsize : 2
109c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_PATIENT -n 10
110c4762a1bSJed Brown
111c4762a1bSJed Brown test:
112c4762a1bSJed Brown suffix: 5
113c4762a1bSJed Brown nsize : 1
114c4762a1bSJed Brown args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5
115c4762a1bSJed Brown
116c4762a1bSJed Brown TEST*/
117