Actual source code: ex158.c
1: static char help[] = "Illustrate how to use mpi FFTW and PETSc-FFTW interface \n\n";
3: /*
4: Usage:
5: mpiexec -n <np> ./ex158 -use_FFTW_interface NO
6: mpiexec -n <np> ./ex158 -use_FFTW_interface YES
7: */
9: #include <petscmat.h>
10: #include <fftw3-mpi.h>
12: int main(int argc, char **args)
13: {
14: PetscMPIInt rank, size;
15: PetscInt N0 = 50, N1 = 20, N = N0 * N1;
16: PetscRandom rdm;
17: PetscScalar a;
18: PetscReal enorm;
19: Vec x, y, z;
20: PetscBool view = PETSC_FALSE, use_interface = PETSC_TRUE;
22: PetscFunctionBeginUser;
23: PetscCall(PetscInitialize(&argc, &args, NULL, help));
24: PetscCheck(!PetscDefined(USE_COMPLEX), PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex");
26: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex158");
27: PetscCall(PetscOptionsBool("-use_FFTW_interface", "Use PETSc-FFTW interface", "ex158", use_interface, &use_interface, NULL));
28: PetscOptionsEnd();
30: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
31: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
33: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
34: PetscCall(PetscRandomSetFromOptions(rdm));
36: if (!use_interface) {
37: /* Use mpi FFTW without PETSc-FFTW interface, 2D case only */
38: /*---------------------------------------------------------*/
39: fftw_plan fplan, bplan;
40: fftw_complex *data_in, *data_out, *data_out2;
41: ptrdiff_t alloc_local, local_n0, local_0_start;
43: if (rank == 0) printf("Use FFTW without PETSc-FFTW interface\n");
44: fftw_mpi_init();
45: N = N0 * N1;
46: alloc_local = fftw_mpi_local_size_2d(N0, N1, PETSC_COMM_WORLD, &local_n0, &local_0_start);
48: data_in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
49: data_out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
50: data_out2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
52: PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_in, &x));
53: PetscCall(PetscObjectSetName((PetscObject)x, "Real Space vector"));
54: PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out, &y));
55: PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
56: PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, (PetscInt)local_n0 * N1, (PetscInt)N, (const PetscScalar *)data_out2, &z));
57: PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
59: fplan = fftw_mpi_plan_dft_2d(N0, N1, data_in, data_out, PETSC_COMM_WORLD, FFTW_FORWARD, FFTW_ESTIMATE);
60: bplan = fftw_mpi_plan_dft_2d(N0, N1, data_out, data_out2, PETSC_COMM_WORLD, FFTW_BACKWARD, FFTW_ESTIMATE);
62: PetscCall(VecSetRandom(x, rdm));
63: if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
65: fftw_execute(fplan);
66: if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
68: fftw_execute(bplan);
70: /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
71: a = 1.0 / (PetscReal)N;
72: PetscCall(VecScale(z, a));
73: if (view) PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD));
74: PetscCall(VecAXPY(z, -1.0, x));
75: PetscCall(VecNorm(z, NORM_1, &enorm));
76: if (enorm > 1.e-11) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm));
78: /* Free spaces */
79: fftw_destroy_plan(fplan);
80: fftw_destroy_plan(bplan);
81: fftw_free(data_in);
82: PetscCall(VecDestroy(&x));
83: fftw_free(data_out);
84: PetscCall(VecDestroy(&y));
85: fftw_free(data_out2);
86: PetscCall(VecDestroy(&z));
88: } else {
89: /* Use PETSc-FFTW interface */
90: /*-------------------------------------------*/
91: PetscInt i, *dim, k, DIM;
92: Mat A;
93: Vec input, output;
95: N = 30;
96: for (i = 2; i < 3; i++) { /* (i=3,4: -- error in VecScatterPetscToFFTW(A,input,x); */
97: DIM = i;
98: PetscCall(PetscMalloc1(i, &dim));
99: for (k = 0; k < i; k++) dim[k] = 30;
100: N *= dim[i - 1];
102: /* Create FFTW object */
103: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Use PETSc-FFTW interface...%d-DIM:%d \n", DIM, N));
104: PetscCall(MatCreateFFT(PETSC_COMM_WORLD, DIM, dim, MATFFTW, &A));
106: /* Create FFTW vectors that are compatible with parallel layout of A */
107: PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
108: PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
109: PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
110: PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
112: /* Create and set PETSc vector */
113: PetscCall(VecCreate(PETSC_COMM_WORLD, &input));
114: PetscCall(VecSetSizes(input, PETSC_DECIDE, N));
115: PetscCall(VecSetFromOptions(input));
116: PetscCall(VecSetRandom(input, rdm));
117: PetscCall(VecDuplicate(input, &output));
118: if (view) PetscCall(VecView(input, PETSC_VIEWER_STDOUT_WORLD));
120: /* Vector input is copied to another vector x using VecScatterPetscToFFTW. This is because the user data
121: can have any parallel layout. But FFTW requires special parallel layout of the data. Hence the original
122: data which is in the vector "input" here, needs to be copied to a vector x, which has the correct parallel
123: layout for FFTW. Also, during parallel real transform, this pads extra zeros automatically
124: at the end of last dimension. This padding is required by FFTW to perform parallel real D.F.T. */
125: PetscCall(VecScatterPetscToFFTW(A, input, x)); /* buggy for dim = 3, 4... */
127: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
128: PetscCall(MatMult(A, x, y));
129: if (view) PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
130: PetscCall(MatMultTranspose(A, y, z));
132: /* Output from Backward DFT needs to be modified to obtain user readable data the routine VecScatterFFTWToPetsc
133: performs the job. In some sense this is the reverse operation of VecScatterPetscToFFTW. This routine gets rid of
134: the extra spaces that were artificially padded to perform real parallel transform. */
135: PetscCall(VecScatterFFTWToPetsc(A, z, output));
137: /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
138: a = 1.0 / (PetscReal)N;
139: PetscCall(VecScale(output, a));
140: if (view) PetscCall(VecView(output, PETSC_VIEWER_STDOUT_WORLD));
141: PetscCall(VecAXPY(output, -1.0, input));
142: PetscCall(VecNorm(output, NORM_1, &enorm));
143: if (enorm > 1.e-09 && rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %e\n", enorm));
145: /* Free spaces */
146: PetscCall(PetscFree(dim));
147: PetscCall(VecDestroy(&input));
148: PetscCall(VecDestroy(&output));
149: PetscCall(VecDestroy(&x));
150: PetscCall(VecDestroy(&y));
151: PetscCall(VecDestroy(&z));
152: PetscCall(MatDestroy(&A));
153: }
154: }
155: PetscCall(PetscRandomDestroy(&rdm));
156: PetscCall(PetscFinalize());
157: return 0;
158: }
160: /*TEST
162: build:
163: requires: !mpiuni fftw !complex
165: test:
166: output_file: output/ex158.out
168: test:
169: suffix: 2
170: nsize: 3
172: TEST*/