1 /* This program illustrates use of parallel real FFT*/
2 static char help[] = "This program illustrates the use of parallel real 3D 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 = 256, N1 = 256, N2 = 256, N3 = 2, dim[4];
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, n1;
16 PetscInt size, rank, n, N, *in, N_factor, NM;
17 PetscScalar *data_fin, value1, one = 1.57, zero = 0.0;
18 PetscScalar a, *x_arr, *y_arr, *z_arr, enorm;
19 Vec fin, fout, fout1, ini, final;
20 PetscRandom rnd;
21 VecScatter vecscat, vecscat1;
22 IS indx1, indx2;
23 PetscInt *indx3, k, l, *indx4;
24 PetscInt low, tempindx, tempindx1;
25
26 PetscFunctionBeginUser;
27 PetscCall(PetscInitialize(&argc, &args, NULL, help));
28 #if defined(PETSC_USE_COMPLEX)
29 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex");
30 #endif
31 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
32 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
33
34 PetscRandomCreate(PETSC_COMM_WORLD, &rnd);
35
36 alloc_local = fftw_mpi_local_size_3d_transposed(N0, N1, N2 / 2 + 1, PETSC_COMM_WORLD, &local_n0, &local_0_start, &local_n1, &local_1_start);
37
38 /* printf("The value alloc_local is %ld from process %d\n",alloc_local,rank); */
39 printf("The value local_n0 is %ld from process %d\n", local_n0, rank);
40 /* printf("The value local_0_start is %ld from process %d\n",local_0_start,rank);*/
41 /* printf("The value local_n1 is %ld from process %d\n",local_n1,rank); */
42 /* printf("The value local_1_start is %ld from process %d\n",local_1_start,rank);*/
43
44 /* Allocate space for input and output arrays */
45
46 in1 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2);
47 in2 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2);
48 out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local);
49
50 N = 2 * N0 * N1 * (N2 / 2 + 1);
51 N_factor = N0 * N1 * N2;
52 n = 2 * local_n0 * N1 * (N2 / 2 + 1);
53 n1 = local_n1 * N0 * 2 * N1;
54
55 /* printf("The value N is %d from process %d\n",N,rank); */
56 /* printf("The value n is %d from process %d\n",n,rank); */
57 /* printf("The value n1 is %d from process %d\n",n1,rank); */
58 /* Creating data vector and accompanying array with VeccreateMPIWithArray */
59 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in1, &fin));
60 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)out, &fout));
61 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in2, &fout1));
62
63 /* VecGetSize(fin,&size); */
64 /* printf("The size is %d\n",size); */
65
66 VecSet(fin, one);
67 VecSet(fout, zero);
68 VecSet(fout1, zero);
69
70 VecAssemblyBegin(fin);
71 VecAssemblyEnd(fin);
72 /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */
73
74 VecGetArray(fin, &x_arr);
75 VecGetArray(fout1, &z_arr);
76 VecGetArray(fout, &y_arr);
77
78 fplan = fftw_mpi_plan_dft_r2c_3d(N0, N1, N2, (double *)x_arr, (fftw_complex *)y_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);
79 bplan = fftw_mpi_plan_dft_c2r_3d(N0, N1, N2, (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 VecCreate(PETSC_COMM_WORLD, &ini);
91 VecCreate(PETSC_COMM_WORLD, &final);
92 VecSetSizes(ini, local_n0 * N1 * N2, N_factor);
93 VecSetSizes(final, local_n0 * N1 * N2, N_factor);
94 /* VecSetSizes(ini,PETSC_DECIDE,N_factor); */
95 /* VecSetSizes(final,PETSC_DECIDE,N_factor); */
96 VecSetFromOptions(ini);
97 VecSetFromOptions(final);
98
99 if (N2 % 2 == 0) NM = N2 + 2;
100 else NM = N2 + 1;
101
102 PetscCall(VecGetOwnershipRange(fin, &low, NULL));
103 printf("The local index is %d from %d\n", low, rank);
104 PetscCall(PetscMalloc1(local_n0 * N1 * N2, &indx3));
105 PetscCall(PetscMalloc1(local_n0 * N1 * N2, &indx4));
106 for (i = 0; i < local_n0; i++) {
107 for (j = 0; j < N1; j++) {
108 for (k = 0; k < N2; k++) {
109 tempindx = i * N1 * N2 + j * N2 + k;
110 tempindx1 = i * N1 * NM + j * NM + k;
111
112 indx3[tempindx] = local_0_start * N1 * N2 + tempindx;
113 indx4[tempindx] = low + tempindx1;
114 }
115 /* printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
116 /* printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
117 }
118 }
119 VecGetValues(fin, local_n0 * N1 * N2, indx4, x_arr);
120 VecSetValues(ini, local_n0 * N1 * N2, indx3, x_arr, INSERT_VALUES);
121 VecAssemblyBegin(ini);
122 VecAssemblyEnd(ini);
123
124 VecGetValues(fout1, local_n0 * N1 * N2, indx4, y_arr);
125 VecSetValues(final, local_n0 * N1 * N2, indx3, y_arr, INSERT_VALUES);
126 VecAssemblyBegin(final);
127 VecAssemblyEnd(final);
128
129 printf("The local index value is %ld from %d", local_n0 * N1 * N2, rank);
130 /*
131 for (i=0;i<N0;i++) {
132 for (j=0;j<N1;j++) {
133 indx=i*N1*NM+j*NM;
134 ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx1);
135 indx=i*N1*N2+j*N2;
136 ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx2);
137 VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
138 VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
139 VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
140 VecScatterCreate(fout1,indx1,final,indx2,&vecscat1);
141 VecScatterBegin(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
142 VecScatterEnd(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
143 }
144 }
145 */
146 a = 1.0 / (PetscReal)N_factor;
147 PetscCall(VecScale(fout1, a));
148 PetscCall(VecScale(final, a));
149
150 VecAssemblyBegin(ini);
151 VecAssemblyEnd(ini);
152
153 VecAssemblyBegin(final);
154 VecAssemblyEnd(final);
155
156 /* VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
157 PetscCall(VecAXPY(final, -1.0, ini));
158 PetscCall(VecNorm(final, NORM_1, &enorm));
159 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z| = %e\n", enorm));
160 fftw_destroy_plan(fplan);
161 fftw_destroy_plan(bplan);
162 fftw_free(in1);
163 PetscCall(VecDestroy(&fin));
164 fftw_free(out);
165 PetscCall(VecDestroy(&fout));
166 fftw_free(in2);
167 PetscCall(VecDestroy(&fout1));
168
169 PetscCall(PetscFinalize());
170 return 0;
171 }
172