1 static char help[] = "Test sequential FFTW interface \n\n";
2
3 /*
4 Compiling the code:
5 This code uses the complex numbers version of PETSc, so configure
6 must be run to enable this
7
8 */
9
10 #include <petscmat.h>
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 Mat A;
21 PetscMPIInt size;
22 PetscInt n = 10, N, ndim = 4, dim[4], DIM, i;
23 Vec x, y, z;
24 PetscScalar s;
25 PetscRandom rdm;
26 PetscReal enorm, tol = PETSC_SMALL;
27 PetscInt func;
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", "ex112", 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 for (i = 0, N = 1; i < DIM; i++) N *= dim[i];
47 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N));
48
49 /* create FFTW object */
50 PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A));
51
52 /* create vectors of length N=n^DIM */
53 PetscCall(MatCreateVecs(A, &x, &y));
54 PetscCall(MatCreateVecs(A, &z, NULL));
55 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
56 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
57 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
58
59 /* set values of space vector x */
60 if (function == RANDOM) {
61 PetscCall(VecSetRandom(x, rdm));
62 } else if (function == CONSTANT) {
63 PetscCall(VecSet(x, 1.0));
64 } else if (function == TANH) {
65 PetscScalar *a;
66 PetscCall(VecGetArray(x, &a));
67 for (i = 0; i < N; ++i) a[i] = tanh((i - N / 2.0) * (10.0 / N));
68 PetscCall(VecRestoreArray(x, &a));
69 }
70 if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
71
72 /* apply FFTW_FORWARD and FFTW_BACKWARD several times on same x, y, and z */
73 for (i = 0; i < 3; i++) {
74 PetscCall(MatMult(A, x, y));
75 if (view && i == 0) PetscCall(VecView(y, PETSC_VIEWER_DRAW_WORLD));
76 PetscCall(MatMultTranspose(A, y, z));
77
78 /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
79 s = 1.0 / (PetscReal)N;
80 PetscCall(VecScale(z, s));
81 if (view && i == 0) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD));
82 PetscCall(VecAXPY(z, -1.0, x));
83 PetscCall(VecNorm(z, NORM_1, &enorm));
84 if (enorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm));
85 }
86
87 /* apply FFTW_FORWARD and FFTW_BACKWARD several times on different x */
88 for (i = 0; i < 3; i++) {
89 PetscCall(VecDestroy(&x));
90 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
91 PetscCall(VecSetRandom(x, rdm));
92
93 PetscCall(MatMult(A, x, y));
94 PetscCall(MatMultTranspose(A, y, z));
95
96 /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
97 s = 1.0 / (PetscReal)N;
98 PetscCall(VecScale(z, s));
99 if (view && i == 0) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD));
100 PetscCall(VecAXPY(z, -1.0, x));
101 PetscCall(VecNorm(z, NORM_1, &enorm));
102 if (enorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of new |x - z| %g\n", (double)enorm));
103 }
104
105 /* free spaces */
106 PetscCall(VecDestroy(&x));
107 PetscCall(VecDestroy(&y));
108 PetscCall(VecDestroy(&z));
109 PetscCall(MatDestroy(&A));
110 }
111 PetscCall(PetscRandomDestroy(&rdm));
112 PetscCall(PetscFinalize());
113 return 0;
114 }
115
116 /*TEST
117
118 build:
119 requires: fftw complex
120
121 test:
122 args: -mat_fftw_plannerflags FFTW_ESTIMATE
123 output_file: output/ex112.out
124
125 test:
126 suffix: 2
127 args: -mat_fftw_plannerflags FFTW_MEASURE
128 output_file: output/ex112.out
129 requires: !defined(PETSC_USE_CXXCOMPLEX)
130
131 test:
132 suffix: 3
133 args: -mat_fftw_plannerflags FFTW_PATIENT
134 output_file: output/ex112.out
135
136 test:
137 suffix: 4
138 args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE
139 output_file: output/ex112.out
140
141 TEST*/
142