Actual source code: ex143.c
1: static char help[] = "Illustrate how to use mpi FFTW and PETSc-FFTW interface \n\n";
3: /*
4: Compiling the code:
5: This code uses the complex numbers version of PETSc, so configure
6: must be run to enable this
8: Usage:
9: mpiexec -n <np> ./ex143 -use_FFTW_interface NO
10: mpiexec -n <np> ./ex143 -use_FFTW_interface YES
11: */
13: #include <petscmat.h>
14: #include <fftw3-mpi.h>
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;
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();
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));
37: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
38: PetscCall(PetscRandomSetFromOptions(rdm));
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;
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);
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);
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"));
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);
67: PetscCall(VecSetRandom(x, rdm));
68: if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
70: fftw_execute(fplan);
71: if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
73: fftw_execute(bplan);
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));
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));
93: } else {
94: /* Use PETSc-FFTW interface */
95: /*-------------------------------------------*/
96: PetscInt i, *dim, k;
97: Mat A;
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];
106: /* Create FFTW object */
107: if (rank == 0) printf("Use PETSc-FFTW interface...%d-DIM: %d\n", (int)DIM, (int)N);
109: PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));
111: /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */
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"));
118: /* Set values of space vector x */
119: PetscCall(VecSetRandom(x, rdm));
121: if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
123: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
124: PetscCall(MatMult(A, x, y));
125: if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
127: PetscCall(MatMultTranspose(A, y, z));
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));
137: PetscCall(VecDestroy(&x));
138: PetscCall(VecDestroy(&y));
139: PetscCall(VecDestroy(&z));
140: PetscCall(MatDestroy(&A));
142: PetscCall(PetscFree(dim));
143: }
144: }
146: PetscCall(PetscRandomDestroy(&rdm));
147: PetscCall(PetscFinalize());
148: return 0;
149: }
151: /*TEST
153: build:
154: requires: !mpiuni fftw complex
156: test:
157: output_file: output/ex143.out
159: test:
160: suffix: 2
161: nsize: 3
162: output_file: output/ex143.out
164: TEST*/