1 static char help[] = "Test duplication/destruction of FFTW vecs \n\n";
2
3 /*
4 Compiling the code:
5 This code uses the FFTW interface.
6 Use one of the options below to configure:
7 --with-fftw-dir=/.... or --download-fftw
8 Usage:
9 mpiexec -np <np> ./ex228
10 */
11
12 #include <petscmat.h>
main(int argc,char ** args)13 int main(int argc, char **args)
14 {
15 Mat A; /* FFT Matrix */
16 Vec x, y, z; /* Work vectors */
17 Vec x1, y1, z1; /* Duplicate vectors */
18 PetscRandom rdm; /* for creating random input */
19 PetscScalar a; /* used to scale output */
20 PetscReal enorm; /* norm for sanity check */
21 PetscInt n = 10, N = 1; /* FFT dimension params */
22 PetscInt DIM, dim[5]; /* FFT params */
23
24 PetscFunctionBeginUser;
25 PetscCall(PetscInitialize(&argc, &args, NULL, help));
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
27
28 /* To create random input vector */
29 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
30 PetscCall(PetscRandomSetFromOptions(rdm));
31
32 /* Iterate over dimensions, use PETSc-FFTW interface */
33 for (PetscInt i = 1; i < 5; i++) {
34 DIM = i;
35 N = 1;
36 for (PetscInt k = 0; k < i; k++) {
37 dim[k] = n;
38 N *= n;
39 }
40
41 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N));
42
43 /* create FFTW object */
44 PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A));
45 /* create vectors of length N */
46 PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
47
48 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
49 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
50 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
51
52 /* Test vector duplication*/
53 PetscCall(VecDuplicate(x, &x1));
54 PetscCall(VecDuplicate(y, &y1));
55 PetscCall(VecDuplicate(z, &z1));
56
57 /* Set values of space vector x, copy to duplicate */
58 PetscCall(VecSetRandom(x, rdm));
59 PetscCall(VecCopy(x, x1));
60
61 /* Apply FFTW_FORWARD and FFTW_BACKWARD */
62 PetscCall(MatMult(A, x, y));
63 PetscCall(MatMultTranspose(A, y, z));
64
65 /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
66 PetscCall(MatMult(A, x1, y1));
67 PetscCall(MatMultTranspose(A, y1, z1));
68
69 /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
70 a = 1.0 / (PetscReal)N;
71 PetscCall(VecScale(z1, a));
72 PetscCall(VecAXPY(z1, -1.0, x));
73 PetscCall(VecNorm(z1, NORM_1, &enorm));
74 if (enorm > 1.e-9) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z1| %g dimension %" PetscInt_FMT "\n", enorm, i));
75
76 /* free spaces */
77 PetscCall(VecDestroy(&x1));
78 PetscCall(VecDestroy(&y1));
79 PetscCall(VecDestroy(&z1));
80 PetscCall(VecDestroy(&x));
81 PetscCall(VecDestroy(&y));
82 PetscCall(VecDestroy(&z));
83 PetscCall(MatDestroy(&A));
84 }
85
86 PetscCall(PetscRandomDestroy(&rdm));
87 PetscCall(PetscFinalize());
88 return 0;
89 }
90
91 /*TEST
92
93 build:
94 requires: fftw complex
95
96 test:
97 suffix: 2
98 nsize : 4
99 args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16
100
101 test:
102 suffix: 3
103 nsize : 2
104 args: -mat_fftw_plannerflags FFTW_MEASURE -n 12
105
106 test:
107 suffix: 4
108 nsize : 2
109 args: -mat_fftw_plannerflags FFTW_PATIENT -n 10
110
111 test:
112 suffix: 5
113 nsize : 1
114 args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5
115
116 TEST*/
117