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