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;
27: PetscInitialize(&argc, &args, (char *)0, help);
28: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex143");
29: PetscOptionsBool("-vec_view draw", "View the vectors", "ex143", view, &view, NULL);
30: PetscOptionsBool("-use_FFTW_interface", "Use PETSc-FFTW interface", "ex143", use_interface, &use_interface, NULL);
31: PetscOptionsEnd();
33: PetscOptionsGetBool(NULL, NULL, "-use_FFTW_interface", &use_interface, NULL);
34: MPI_Comm_size(PETSC_COMM_WORLD, &size);
35: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
37: PetscRandomCreate(PETSC_COMM_WORLD, &rdm);
38: 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) 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: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_in, &x);
58: PetscObjectSetName((PetscObject)x, "Real Space vector");
59: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out, &y);
60: PetscObjectSetName((PetscObject)y, "Frequency space vector");
61: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out2, &z);
62: 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: VecSetRandom(x, rdm);
68: if (view) VecView(x, PETSC_VIEWER_STDOUT_WORLD);
70: fftw_execute(fplan);
71: if (view) 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: VecScale(z, a);
78: if (view) VecView(z, PETSC_VIEWER_STDOUT_WORLD);
79: VecAXPY(z, -1.0, x);
80: VecNorm(z, NORM_1, &enorm);
81: if (enorm > 1.e-11 && rank == 0) 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: VecDestroy(&x);
88: fftw_free(data_out);
89: VecDestroy(&y);
90: fftw_free(data_out2);
91: 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: 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: MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A);
111: /* Create vectors that are compatible with parallel layout of A - must call MatCreateVecs()! */
113: MatCreateVecsFFTW(A, &x, &y, &z);
114: PetscObjectSetName((PetscObject)x, "Real space vector");
115: PetscObjectSetName((PetscObject)y, "Frequency space vector");
116: PetscObjectSetName((PetscObject)z, "Reconstructed vector");
118: /* Set values of space vector x */
119: VecSetRandom(x, rdm);
121: if (view) VecView(x, PETSC_VIEWER_STDOUT_WORLD);
123: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
124: MatMult(A, x, y);
125: if (view) VecView(y, PETSC_VIEWER_STDOUT_WORLD);
127: 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: VecScale(z, a);
132: if (view) VecView(z, PETSC_VIEWER_STDOUT_WORLD);
133: VecAXPY(z, -1.0, x);
134: VecNorm(z, NORM_1, &enorm);
135: if (enorm > 1.e-9 && rank == 0) PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %e\n", enorm);
137: VecDestroy(&x);
138: VecDestroy(&y);
139: VecDestroy(&z);
140: MatDestroy(&A);
142: PetscFree(dim);
143: }
144: }
146: PetscRandomDestroy(&rdm);
147: 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*/