1 static char help[] = "Test sequential FFTW convolution\n\n";
2
3 /*
4 Compiling the code:
5 This code uses the complex numbers, so configure must be given --with-scalar-type=complex to enable this
6 */
7
8 #include <petscmat.h>
9
main(int argc,char ** args)10 int main(int argc, char **args)
11 {
12 typedef enum {
13 RANDOM,
14 CONSTANT,
15 TANH,
16 NUM_FUNCS
17 } FuncType;
18 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
19 Mat A;
20 PetscMPIInt size;
21 PetscInt n = 10, N, ndim = 4, dim[4], DIM, i, j;
22 Vec w, x, y1, y2, z1, z2;
23 PetscScalar *a, *a2, *a3;
24 PetscScalar s;
25 PetscRandom rdm;
26 PetscReal enorm;
27 PetscInt func = 0;
28 FuncType function = RANDOM;
29 PetscBool view = PETSC_FALSE;
30
31 PetscFunctionBeginUser;
32 PetscCall(PetscInitialize(&argc, &args, NULL, help));
33 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
34 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
35 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");
36 PetscCall(PetscOptionsEList("-function", "Function type", "ex121", funcNames, NUM_FUNCS, funcNames[function], &func, NULL));
37 PetscCall(PetscOptionsBool("-vec_view draw", "View the functions", "ex112", view, &view, NULL));
38 function = (FuncType)func;
39 PetscOptionsEnd();
40
41 for (DIM = 0; DIM < ndim; DIM++) dim[DIM] = n; /* size of transformation in DIM-dimension */
42 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
43 PetscCall(PetscRandomSetFromOptions(rdm));
44
45 for (DIM = 1; DIM < 5; DIM++) {
46 /* create vectors of length N=n^DIM */
47 for (i = 0, N = 1; i < DIM; i++) N *= dim[i];
48 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N));
49 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
50 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
51 PetscCall(VecDuplicate(x, &w));
52 PetscCall(PetscObjectSetName((PetscObject)w, "Window vector"));
53 PetscCall(VecDuplicate(x, &y1));
54 PetscCall(PetscObjectSetName((PetscObject)y1, "Frequency space vector"));
55 PetscCall(VecDuplicate(x, &y2));
56 PetscCall(PetscObjectSetName((PetscObject)y2, "Frequency space window vector"));
57 PetscCall(VecDuplicate(x, &z1));
58 PetscCall(PetscObjectSetName((PetscObject)z1, "Reconstructed convolution"));
59 PetscCall(VecDuplicate(x, &z2));
60 PetscCall(PetscObjectSetName((PetscObject)z2, "Real space convolution"));
61
62 if (function == RANDOM) {
63 PetscCall(VecSetRandom(x, rdm));
64 } else if (function == CONSTANT) {
65 PetscCall(VecSet(x, 1.0));
66 } else if (function == TANH) {
67 PetscCall(VecGetArray(x, &a));
68 for (i = 0; i < N; ++i) a[i] = tanh((i - N / 2.0) * (10.0 / N));
69 PetscCall(VecRestoreArray(x, &a));
70 }
71 if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
72
73 /* Create window function */
74 PetscCall(VecGetArray(w, &a));
75 for (i = 0; i < N; ++i) {
76 /* Step Function */
77 a[i] = (i > N / 4 && i < 3 * N / 4) ? 1.0 : 0.0;
78 /* Delta Function */
79 /*a[i] = (i == N/2)? 1.0: 0.0; */
80 }
81 PetscCall(VecRestoreArray(w, &a));
82 if (view) PetscCall(VecView(w, PETSC_VIEWER_DRAW_WORLD));
83
84 /* create FFTW object */
85 PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A));
86
87 /* Convolve x with w*/
88 PetscCall(MatMult(A, x, y1));
89 PetscCall(MatMult(A, w, y2));
90 PetscCall(VecPointwiseMult(y1, y1, y2));
91 if (view && i == 0) PetscCall(VecView(y1, PETSC_VIEWER_DRAW_WORLD));
92 PetscCall(MatMultTranspose(A, y1, z1));
93
94 /* Compute the real space convolution */
95 PetscCall(VecGetArray(x, &a));
96 PetscCall(VecGetArray(w, &a2));
97 PetscCall(VecGetArray(z2, &a3));
98 for (i = 0; i < N; ++i) {
99 /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/
100
101 a3[i] = 0.0;
102 for (j = -N / 2 + 1; j < N / 2; ++j) {
103 PetscInt xpInd = (j < 0) ? N + j : j;
104 PetscInt diffInd = (i - j < 0) ? N - (j - i) : (i - j > N - 1) ? i - j - N : i - j;
105
106 a3[i] += a[xpInd] * a2[diffInd];
107 }
108 }
109 PetscCall(VecRestoreArray(x, &a));
110 PetscCall(VecRestoreArray(w, &a2));
111 PetscCall(VecRestoreArray(z2, &a3));
112
113 /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */
114 s = 1.0 / (PetscReal)N;
115 PetscCall(VecScale(z1, s));
116 if (view) PetscCall(VecView(z1, PETSC_VIEWER_DRAW_WORLD));
117 if (view) PetscCall(VecView(z2, PETSC_VIEWER_DRAW_WORLD));
118 PetscCall(VecAXPY(z1, -1.0, z2));
119 PetscCall(VecNorm(z1, NORM_1, &enorm));
120 if (enorm > 1.e-11) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |z1 - z2| %g\n", (double)enorm));
121
122 /* free spaces */
123 PetscCall(VecDestroy(&x));
124 PetscCall(VecDestroy(&y1));
125 PetscCall(VecDestroy(&y2));
126 PetscCall(VecDestroy(&z1));
127 PetscCall(VecDestroy(&z2));
128 PetscCall(VecDestroy(&w));
129 PetscCall(MatDestroy(&A));
130 }
131 PetscCall(PetscRandomDestroy(&rdm));
132 PetscCall(PetscFinalize());
133 return 0;
134 }
135
136 /*TEST
137
138 build:
139 requires: fftw complex
140
141 test:
142 output_file: output/ex121.out
143 TODO: Example or FFTW interface is broken
144
145 TEST*/
146