Actual source code: ex121.c
1: static char help[] = "Test sequential FFTW convolution\n\n";
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: */
8: #include <petscmat.h>
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;
32: PetscInitialize(&argc, &args, (char *)0, help);
33: MPI_Comm_size(PETSC_COMM_WORLD, &size);
35: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");
36: PetscOptionsEList("-function", "Function type", "ex121", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
37: PetscOptionsBool("-vec_view draw", "View the functions", "ex112", view, &view, NULL);
38: function = (FuncType)func;
39: PetscOptionsEnd();
41: for (DIM = 0; DIM < ndim; DIM++) { dim[DIM] = n; /* size of transformation in DIM-dimension */ }
42: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
43: PetscRandomSetFromOptions(rdm);
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: PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N);
49: VecCreateSeq(PETSC_COMM_SELF, N, &x);
50: PetscObjectSetName((PetscObject)x, "Real space vector");
51: VecDuplicate(x, &w);
52: PetscObjectSetName((PetscObject)w, "Window vector");
53: VecDuplicate(x, &y1);
54: PetscObjectSetName((PetscObject)y1, "Frequency space vector");
55: VecDuplicate(x, &y2);
56: PetscObjectSetName((PetscObject)y2, "Frequency space window vector");
57: VecDuplicate(x, &z1);
58: PetscObjectSetName((PetscObject)z1, "Reconstructed convolution");
59: VecDuplicate(x, &z2);
60: PetscObjectSetName((PetscObject)z2, "Real space convolution");
62: if (function == RANDOM) {
63: VecSetRandom(x, rdm);
64: } else if (function == CONSTANT) {
65: VecSet(x, 1.0);
66: } else if (function == TANH) {
67: VecGetArray(x, &a);
68: for (i = 0; i < N; ++i) a[i] = tanh((i - N / 2.0) * (10.0 / N));
69: VecRestoreArray(x, &a);
70: }
71: if (view) VecView(x, PETSC_VIEWER_DRAW_WORLD);
73: /* Create window function */
74: 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: VecRestoreArray(w, &a);
82: if (view) VecView(w, PETSC_VIEWER_DRAW_WORLD);
84: /* create FFTW object */
85: MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A);
87: /* Convolve x with w*/
88: MatMult(A, x, y1);
89: MatMult(A, w, y2);
90: VecPointwiseMult(y1, y1, y2);
91: if (view && i == 0) VecView(y1, PETSC_VIEWER_DRAW_WORLD);
92: MatMultTranspose(A, y1, z1);
94: /* Compute the real space convolution */
95: VecGetArray(x, &a);
96: VecGetArray(w, &a2);
97: VecGetArray(z2, &a3);
98: for (i = 0; i < N; ++i) {
99: /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/
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;
106: a3[i] += a[xpInd] * a2[diffInd];
107: }
108: }
109: VecRestoreArray(x, &a);
110: VecRestoreArray(w, &a2);
111: VecRestoreArray(z2, &a3);
113: /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */
114: s = 1.0 / (PetscReal)N;
115: VecScale(z1, s);
116: if (view) VecView(z1, PETSC_VIEWER_DRAW_WORLD);
117: if (view) VecView(z2, PETSC_VIEWER_DRAW_WORLD);
118: VecAXPY(z1, -1.0, z2);
119: VecNorm(z1, NORM_1, &enorm);
120: if (enorm > 1.e-11) PetscPrintf(PETSC_COMM_SELF, " Error norm of |z1 - z2| %g\n", (double)enorm);
122: /* free spaces */
123: VecDestroy(&x);
124: VecDestroy(&y1);
125: VecDestroy(&y2);
126: VecDestroy(&z1);
127: VecDestroy(&z2);
128: VecDestroy(&w);
129: MatDestroy(&A);
130: }
131: PetscRandomDestroy(&rdm);
132: PetscFinalize();
133: return 0;
134: }
136: /*TEST
138: build:
139: requires: fftw complex
141: test:
142: output_file: output/ex121.out
143: TODO: Example or FFTW interface is broken
145: TEST*/