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