1 static char help[] = "Illustrate how to use mpi FFTW and PETSc-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 Usage:
9 mpiexec -n <np> ./ex143 -use_FFTW_interface NO
10 mpiexec -n <np> ./ex143 -use_FFTW_interface YES
11 */
12
13 #include <petscmat.h>
14 #include <fftw3-mpi.h>
15
main(int argc,char ** args)16 int main(int argc, char **args)
17 {
18 PetscMPIInt rank, size;
19 PetscInt N0 = 50, N1 = 20, N = N0 * N1, DIM;
20 PetscRandom rdm;
21 PetscScalar a;
22 PetscReal enorm;
23 Vec x, y, z;
24 PetscBool view = PETSC_FALSE, use_interface = PETSC_TRUE;
25
26 PetscFunctionBeginUser;
27 PetscCall(PetscInitialize(&argc, &args, NULL, help));
28 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex143");
29 PetscCall(PetscOptionsBool("-vec_view draw", "View the vectors", "ex143", view, &view, NULL));
30 PetscCall(PetscOptionsBool("-use_FFTW_interface", "Use PETSc-FFTW interface", "ex143", use_interface, &use_interface, NULL));
31 PetscOptionsEnd();
32
33 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_FFTW_interface", &use_interface, NULL));
34 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
35 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
36
37 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
38 PetscCall(PetscRandomSetFromOptions(rdm));
39
40 if (!use_interface) {
41 /* Use mpi FFTW without PETSc-FFTW interface, 2D case only */
42 /*---------------------------------------------------------*/
43 fftw_plan fplan, bplan;
44 fftw_complex *data_in, *data_out, *data_out2;
45 ptrdiff_t alloc_local, local_n0, local_0_start;
46
47 DIM = 2;
48 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Use FFTW without PETSc-FFTW interface, DIM %" PetscInt_FMT "\n", DIM));
49 fftw_mpi_init();
50 N = N0 * N1;
51 alloc_local = fftw_mpi_local_size_2d(N0, N1, PETSC_COMM_WORLD, &local_n0, &local_0_start);
52
53 data_in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
54 data_out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
55 data_out2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
56
57 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_in, &x));
58 PetscCall(PetscObjectSetName((PetscObject)x, "Real Space vector"));
59 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out, &y));
60 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
61 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out2, &z));
62 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
63
64 fplan = fftw_mpi_plan_dft_2d(N0, N1, data_in, data_out, PETSC_COMM_WORLD, FFTW_FORWARD, FFTW_ESTIMATE);
65 bplan = fftw_mpi_plan_dft_2d(N0, N1, data_out, data_out2, PETSC_COMM_WORLD, FFTW_BACKWARD, FFTW_ESTIMATE);
66
67 PetscCall(VecSetRandom(x, rdm));
68 if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
69
70 fftw_execute(fplan);
71 if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
72
73 fftw_execute(bplan);
74
75 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
76 a = 1.0 / (PetscReal)N;
77 PetscCall(VecScale(z, a));
78 if (view) PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD));
79 PetscCall(VecAXPY(z, -1.0, x));
80 PetscCall(VecNorm(z, NORM_1, &enorm));
81 if (enorm > 1.e-11 && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm));
82
83 /* Free spaces */
84 fftw_destroy_plan(fplan);
85 fftw_destroy_plan(bplan);
86 fftw_free(data_in);
87 PetscCall(VecDestroy(&x));
88 fftw_free(data_out);
89 PetscCall(VecDestroy(&y));
90 fftw_free(data_out2);
91 PetscCall(VecDestroy(&z));
92
93 } else {
94 /* Use PETSc-FFTW interface */
95 /*-------------------------------------------*/
96 PetscInt i, *dim, k;
97 Mat A;
98
99 N = 1;
100 for (i = 1; i < 5; i++) {
101 DIM = i;
102 PetscCall(PetscMalloc1(i, &dim));
103 for (k = 0; k < i; k++) dim[k] = 30;
104 N *= dim[i - 1];
105
106 /* Create FFTW object */
107 if (rank == 0) printf("Use PETSc-FFTW interface...%d-DIM: %d\n", (int)DIM, (int)N);
108
109 PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));
110
111 /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */
112
113 PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
114 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
115 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
116 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
117
118 /* Set values of space vector x */
119 PetscCall(VecSetRandom(x, rdm));
120
121 if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
122
123 /* Apply FFTW_FORWARD and FFTW_BACKWARD */
124 PetscCall(MatMult(A, x, y));
125 if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
126
127 PetscCall(MatMultTranspose(A, y, z));
128
129 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
130 a = 1.0 / (PetscReal)N;
131 PetscCall(VecScale(z, a));
132 if (view) PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD));
133 PetscCall(VecAXPY(z, -1.0, x));
134 PetscCall(VecNorm(z, NORM_1, &enorm));
135 if (enorm > 1.e-9 && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %e\n", enorm));
136
137 PetscCall(VecDestroy(&x));
138 PetscCall(VecDestroy(&y));
139 PetscCall(VecDestroy(&z));
140 PetscCall(MatDestroy(&A));
141
142 PetscCall(PetscFree(dim));
143 }
144 }
145
146 PetscCall(PetscRandomDestroy(&rdm));
147 PetscCall(PetscFinalize());
148 return 0;
149 }
150
151 /*TEST
152
153 build:
154 requires: !mpiuni fftw complex
155
156 test:
157 output_file: output/ex143.out
158
159 test:
160 suffix: 2
161 nsize: 3
162 output_file: output/ex143.out
163
164 TEST*/
165