1 static char help[] = "Test sequential r2c/c2r FFTW without PETSc interface \n\n";
2
3 /*
4 Compiling the code:
5 This code uses the real numbers version of PETSc
6 */
7
8 #include <petscmat.h>
9 #include <fftw3.h>
10
main(int argc,char ** args)11 int main(int argc, char **args)
12 {
13 typedef enum {
14 RANDOM,
15 CONSTANT,
16 TANH,
17 NUM_FUNCS
18 } FuncType;
19 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
20 PetscMPIInt size;
21 int n = 10, N, Ny, ndim = 4, i, dim[4], DIM;
22 Vec x, y, z;
23 PetscScalar s;
24 PetscRandom rdm;
25 PetscReal enorm;
26 PetscInt func = RANDOM;
27 FuncType function = RANDOM;
28 PetscBool view = PETSC_FALSE;
29 PetscScalar *x_array, *y_array, *z_array;
30 fftw_plan fplan, bplan;
31
32 PetscFunctionBeginUser;
33 PetscCall(PetscInitialize(&argc, &args, NULL, help));
34 #if defined(PETSC_USE_COMPLEX)
35 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers");
36 #endif
37
38 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
39 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
40 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex142");
41 PetscCall(PetscOptionsEList("-function", "Function type", "ex142", funcNames, NUM_FUNCS, funcNames[function], &func, NULL));
42 PetscCall(PetscOptionsBool("-vec_view draw", "View the functions", "ex142", view, &view, NULL));
43 function = (FuncType)func;
44 PetscOptionsEnd();
45
46 for (DIM = 0; DIM < ndim; DIM++) dim[DIM] = n; /* size of real space vector in DIM-dimension */
47 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
48 PetscCall(PetscRandomSetFromOptions(rdm));
49
50 for (DIM = 1; DIM < 5; DIM++) {
51 /* create vectors of length N=dim[0]*dim[1]* ...*dim[DIM-1] */
52 /*----------------------------------------------------------*/
53 N = Ny = 1;
54 for (i = 0; i < DIM - 1; i++) N *= dim[i];
55 Ny = N;
56 Ny *= 2 * (dim[DIM - 1] / 2 + 1); /* add padding elements to output vector y */
57 N *= dim[DIM - 1];
58
59 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N));
60 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
61 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
62
63 PetscCall(VecCreateSeq(PETSC_COMM_SELF, Ny, &y));
64 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
65
66 PetscCall(VecDuplicate(x, &z));
67 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
68
69 /* Set fftw plan */
70 /*----------------------------------*/
71 PetscCall(VecGetArray(x, &x_array));
72 PetscCall(VecGetArray(y, &y_array));
73 PetscCall(VecGetArray(z, &z_array));
74
75 unsigned int flags = FFTW_ESTIMATE; /*or FFTW_MEASURE */
76 /* The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so such planning
77 should be done before the input is initialized by the user. */
78 PetscCall(PetscPrintf(PETSC_COMM_SELF, "DIM: %d, N %d, Ny %d\n", DIM, N, Ny));
79
80 switch (DIM) {
81 case 1:
82 fplan = fftw_plan_dft_r2c_1d(dim[0], (double *)x_array, (fftw_complex *)y_array, flags);
83 bplan = fftw_plan_dft_c2r_1d(dim[0], (fftw_complex *)y_array, (double *)z_array, flags);
84 break;
85 case 2:
86 fplan = fftw_plan_dft_r2c_2d(dim[0], dim[1], (double *)x_array, (fftw_complex *)y_array, flags);
87 bplan = fftw_plan_dft_c2r_2d(dim[0], dim[1], (fftw_complex *)y_array, (double *)z_array, flags);
88 break;
89 case 3:
90 fplan = fftw_plan_dft_r2c_3d(dim[0], dim[1], dim[2], (double *)x_array, (fftw_complex *)y_array, flags);
91 bplan = fftw_plan_dft_c2r_3d(dim[0], dim[1], dim[2], (fftw_complex *)y_array, (double *)z_array, flags);
92 break;
93 default:
94 fplan = fftw_plan_dft_r2c(DIM, (int *)dim, (double *)x_array, (fftw_complex *)y_array, flags);
95 bplan = fftw_plan_dft_c2r(DIM, (int *)dim, (fftw_complex *)y_array, (double *)z_array, flags);
96 break;
97 }
98
99 PetscCall(VecRestoreArray(x, &x_array));
100 PetscCall(VecRestoreArray(y, &y_array));
101 PetscCall(VecRestoreArray(z, &z_array));
102
103 /* Initialize Real space vector x:
104 The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so planning
105 should be done before the input is initialized by the user.
106 --------------------------------------------------------*/
107 if (function == RANDOM) {
108 PetscCall(VecSetRandom(x, rdm));
109 } else if (function == CONSTANT) {
110 PetscCall(VecSet(x, 1.0));
111 } else if (function == TANH) {
112 PetscCall(VecGetArray(x, &x_array));
113 for (i = 0; i < N; ++i) x_array[i] = tanh((i - N / 2.0) * (10.0 / N));
114 PetscCall(VecRestoreArray(x, &x_array));
115 }
116 if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
117
118 /* FFT - also test repeated transformation */
119 /*-------------------------------------------*/
120 PetscCall(VecGetArray(x, &x_array));
121 PetscCall(VecGetArray(y, &y_array));
122 PetscCall(VecGetArray(z, &z_array));
123 for (i = 0; i < 4; i++) {
124 /* FFTW_FORWARD */
125 fftw_execute(fplan);
126
127 /* FFTW_BACKWARD: destroys its input array 'y_array' even for out-of-place transforms! */
128 fftw_execute(bplan);
129 }
130 PetscCall(VecRestoreArray(x, &x_array));
131 PetscCall(VecRestoreArray(y, &y_array));
132 PetscCall(VecRestoreArray(z, &z_array));
133
134 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
135 /*------------------------------------------------------------------*/
136 s = 1.0 / (PetscReal)N;
137 PetscCall(VecScale(z, s));
138 if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
139 if (view) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD));
140 PetscCall(VecAXPY(z, -1.0, x));
141 PetscCall(VecNorm(z, NORM_1, &enorm));
142 if (enorm > 1.e-11) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm));
143
144 /* free spaces */
145 fftw_destroy_plan(fplan);
146 fftw_destroy_plan(bplan);
147 PetscCall(VecDestroy(&x));
148 PetscCall(VecDestroy(&y));
149 PetscCall(VecDestroy(&z));
150 }
151 PetscCall(PetscRandomDestroy(&rdm));
152 PetscCall(PetscFinalize());
153 return 0;
154 }
155
156 /*TEST
157
158 build:
159 requires: fftw !complex
160
161 test:
162 output_file: output/ex142.out
163
164 TEST*/
165