1 /* This program illustrates use of parallel real FFT */
2 static char help[] = "This program illustrates the use of parallel real multi-dimensional fftw (without PETSc interface)";
3 #include <petscmat.h>
4 #include <fftw3.h>
5 #include <fftw3-mpi.h>
6
main(int argc,char ** args)7 int main(int argc, char **args)
8 {
9 ptrdiff_t N0 = 2, N1 = 2, N2 = 2, N3 = 2, dim[4], N, D;
10 fftw_plan bplan, fplan;
11 fftw_complex *out;
12 double *in1, *in2;
13 ptrdiff_t alloc_local, local_n0, local_0_start;
14 ptrdiff_t local_n1, local_1_start;
15 PetscInt i, j, indx[100], n1;
16 PetscInt size, rank, n, *in, N_factor;
17 PetscScalar *data_fin, value1, one = 1.0, zero = 0.0;
18 PetscScalar a, *x_arr, *y_arr, *z_arr, enorm;
19 Vec fin, fout, fout1, x, y;
20 PetscRandom rnd;
21
22 PetscFunctionBeginUser;
23 PetscCall(PetscInitialize(&argc, &args, NULL, help));
24 #if defined(PETSC_USE_COMPLEX)
25 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex");
26 #endif
27 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
28 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
29
30 PetscRandomCreate(PETSC_COMM_WORLD, &rnd);
31 D = 4;
32 dim[0] = N0;
33 dim[1] = N1;
34 dim[2] = N2;
35 dim[3] = N3 / 2 + 1;
36
37 alloc_local = fftw_mpi_local_size_transposed(D, dim, PETSC_COMM_WORLD, &local_n0, &local_0_start, &local_n1, &local_1_start);
38
39 printf("The value alloc_local is %ld from process %d\n", alloc_local, rank);
40 printf("The value local_n0 is %ld from process %d\n", local_n0, rank);
41 printf("The value local_0_start is %ld from process %d\n", local_0_start, rank);
42 printf("The value local_n1 is %ld from process %d\n", local_n1, rank);
43 printf("The value local_1_start is %ld from process %d\n", local_1_start, rank);
44
45 /* Allocate space for input and output arrays */
46
47 in1 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2);
48 in2 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2);
49 out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
50
51 N = 2 * N0 * N1 * N2 * (N3 / 2 + 1);
52 N_factor = N0 * N1 * N2 * N3;
53 n = 2 * local_n0 * N1 * N2 * (N3 / 2 + 1);
54 n1 = local_n1 * N0 * 2 * N1 * N2;
55
56 /* printf("The value N is %d from process %d\n",N,rank); */
57 /* printf("The value n is %d from process %d\n",n,rank); */
58 /* printf("The value n1 is %d from process %d\n",n1,rank); */
59 /* Creating data vector and accompanying array with VeccreateMPIWithArray */
60 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in1, &fin));
61 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)out, &fout));
62 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in2, &fout1));
63
64 /* VecGetSize(fin,&size); */
65 /* printf("The size is %d\n",size); */
66
67 VecSet(fin, one);
68 /* VecAssemblyBegin(fin); */
69 /* VecAssemblyEnd(fin); */
70 /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */
71
72 VecGetArray(fin, &x_arr);
73 VecGetArray(fout1, &z_arr);
74 VecGetArray(fout, &y_arr);
75
76 dim[3] = N3;
77
78 fplan = fftw_mpi_plan_dft_r2c(D, dim, (double *)x_arr, (fftw_complex *)y_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);
79 bplan = fftw_mpi_plan_dft_c2r(D, dim, (fftw_complex *)y_arr, (double *)z_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);
80
81 fftw_execute(fplan);
82 fftw_execute(bplan);
83
84 VecRestoreArray(fin, &x_arr);
85 VecRestoreArray(fout1, &z_arr);
86 VecRestoreArray(fout, &y_arr);
87
88 /* a = 1.0/(PetscReal)N_factor; */
89 /* PetscCall(VecScale(fout1,a)); */
90
91 VecAssemblyBegin(fout1);
92 VecAssemblyEnd(fout1);
93
94 VecView(fout1, PETSC_VIEWER_STDOUT_WORLD);
95
96 fftw_destroy_plan(fplan);
97 fftw_destroy_plan(bplan);
98 fftw_free(in1);
99 PetscCall(VecDestroy(&fin));
100 fftw_free(out);
101 PetscCall(VecDestroy(&fout));
102 fftw_free(in2);
103 PetscCall(VecDestroy(&fout1));
104
105 PetscCall(PetscFinalize());
106 return 0;
107 }
108