Actual source code: ex144.c
1: /* This program illustrates use of parallel real FFT */
2: static char help[] = "This program illustrates the use of parallel real 2D fft using fftw without PETSc interface";
4: #include <petscmat.h>
5: #include <fftw3.h>
6: #include <fftw3-mpi.h>
8: int main(int argc, char **args)
9: {
10: const ptrdiff_t N0 = 2056, N1 = 2056;
11: fftw_plan bplan, fplan;
12: fftw_complex *out;
13: double *in1, *in2;
14: ptrdiff_t alloc_local, local_n0, local_0_start;
15: ptrdiff_t local_n1, local_1_start;
16: PetscInt i, j;
17: PetscMPIInt size, rank;
18: int n, N, N_factor, NM;
19: PetscScalar one = 2.0, zero = 0.5;
20: PetscScalar two = 4.0, three = 8.0, four = 16.0;
21: PetscScalar a, *x_arr, *y_arr, *z_arr;
22: PetscReal enorm;
23: Vec fin, fout, fout1;
24: Vec ini, final;
25: PetscRandom rnd;
26: PetscInt *indx3, tempindx, low, *indx4, tempindx1;
29: PetscInitialize(&argc, &args, (char *)0, help);
30: MPI_Comm_size(PETSC_COMM_WORLD, &size);
31: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
33: PetscRandomCreate(PETSC_COMM_WORLD, &rnd);
35: alloc_local = fftw_mpi_local_size_2d_transposed(N0, N1 / 2 + 1, PETSC_COMM_WORLD, &local_n0, &local_0_start, &local_n1, &local_1_start);
36: #if defined(DEBUGGING)
37: printf("The value alloc_local is %ld from process %d\n", alloc_local, rank);
38: printf("The value local_n0 is %ld from process %d\n", local_n0, rank);
39: printf("The value local_0_start is %ld from process %d\n", local_0_start, rank);
40: /* printf("The value local_n1 is %ld from process %d\n",local_n1,rank); */
41: /* printf("The value local_1_start is %ld from process %d\n",local_1_start,rank); */
42: /* printf("The value local_n0 is %ld from process %d\n",local_n0,rank); */
43: #endif
45: /* Allocate space for input and output arrays */
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);
50: N = 2 * N0 * (N1 / 2 + 1);
51: N_factor = N0 * N1;
52: n = 2 * local_n0 * (N1 / 2 + 1);
54: /* printf("The value N is %d from process %d\n",N,rank); */
55: /* printf("The value n is %d from process %d\n",n,rank); */
56: /* printf("The value n1 is %d from process %d\n",n1,rank);*/
57: /* Creating data vector and accompanying array with VeccreateMPIWithArray */
58: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in1, &fin);
59: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)out, &fout);
60: VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in2, &fout1);
62: /* Set the vector with random data */
63: VecSet(fin, zero);
64: /* for (i=0;i<N0*N1;i++) */
65: /* { */
66: /* VecSetValues(fin,1,&i,&one,INSERT_VALUES); */
67: /* } */
69: /* VecSet(fin,one); */
70: i = 0;
71: VecSetValues(fin, 1, &i, &one, INSERT_VALUES);
72: i = 1;
73: VecSetValues(fin, 1, &i, &two, INSERT_VALUES);
74: i = 4;
75: VecSetValues(fin, 1, &i, &three, INSERT_VALUES);
76: i = 5;
77: VecSetValues(fin, 1, &i, &four, INSERT_VALUES);
78: VecAssemblyBegin(fin);
79: VecAssemblyEnd(fin);
81: VecSet(fout, zero);
82: VecSet(fout1, zero);
84: /* Get the meaningful portion of array */
85: VecGetArray(fin, &x_arr);
86: VecGetArray(fout1, &z_arr);
87: VecGetArray(fout, &y_arr);
89: fplan = fftw_mpi_plan_dft_r2c_2d(N0, N1, (double *)x_arr, (fftw_complex *)y_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);
90: bplan = fftw_mpi_plan_dft_c2r_2d(N0, N1, (fftw_complex *)y_arr, (double *)z_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);
92: fftw_execute(fplan);
93: fftw_execute(bplan);
95: VecRestoreArray(fin, &x_arr);
96: VecRestoreArray(fout1, &z_arr);
97: VecRestoreArray(fout, &y_arr);
99: /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */
100: VecCreate(PETSC_COMM_WORLD, &ini);
101: VecCreate(PETSC_COMM_WORLD, &final);
102: VecSetSizes(ini, local_n0 * N1, N0 * N1);
103: VecSetSizes(final, local_n0 * N1, N0 * N1);
104: VecSetFromOptions(ini);
105: VecSetFromOptions(final);
107: if (N1 % 2 == 0) {
108: NM = N1 + 2;
109: } else {
110: NM = N1 + 1;
111: }
112: /*printf("The Value of NM is %d",NM); */
113: VecGetOwnershipRange(fin, &low, NULL);
114: /*printf("The local index is %d from %d\n",low,rank); */
115: PetscMalloc1(local_n0 * N1, &indx3);
116: PetscMalloc1(local_n0 * N1, &indx4);
117: for (i = 0; i < local_n0; i++) {
118: for (j = 0; j < N1; j++) {
119: tempindx = i * N1 + j;
120: tempindx1 = i * NM + j;
122: indx3[tempindx] = local_0_start * N1 + tempindx;
123: indx4[tempindx] = low + tempindx1;
124: /* printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
125: /* printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
126: }
127: }
129: PetscMalloc2(local_n0 * N1, &x_arr, local_n0 * N1, &y_arr); /* arr must be allocated for VecGetValues() */
130: VecGetValues(fin, local_n0 * N1, indx4, (PetscScalar *)x_arr);
131: VecSetValues(ini, local_n0 * N1, indx3, x_arr, INSERT_VALUES);
133: VecAssemblyBegin(ini);
134: VecAssemblyEnd(ini);
136: VecGetValues(fout1, local_n0 * N1, indx4, y_arr);
137: VecSetValues(final, local_n0 * N1, indx3, y_arr, INSERT_VALUES);
138: VecAssemblyBegin(final);
139: VecAssemblyEnd(final);
140: PetscFree2(x_arr, y_arr);
142: /*
143: VecScatter vecscat;
144: IS indx1,indx2;
145: for (i=0;i<N0;i++) {
146: indx = i*NM;
147: ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx1);
148: indx = i*N1;
149: ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx2);
150: VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
151: VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
152: VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
153: VecScatterBegin(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
154: VecScatterEnd(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
155: }
156: */
158: a = 1.0 / (PetscReal)N_factor;
159: VecScale(fout1, a);
160: VecScale(final, a);
162: /* VecView(ini,PETSC_VIEWER_STDOUT_WORLD); */
163: /* VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
164: VecAXPY(final, -1.0, ini);
166: VecNorm(final, NORM_1, &enorm);
167: if (enorm > 1.e-10) PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z| = %e\n", enorm);
169: /* Execute fftw with function fftw_execute and destroy it after execution */
170: fftw_destroy_plan(fplan);
171: fftw_destroy_plan(bplan);
172: fftw_free(in1);
173: VecDestroy(&fin);
174: fftw_free(out);
175: VecDestroy(&fout);
176: fftw_free(in2);
177: VecDestroy(&fout1);
179: VecDestroy(&ini);
180: VecDestroy(&final);
182: PetscRandomDestroy(&rnd);
183: PetscFree(indx3);
184: PetscFree(indx4);
185: PetscFinalize();
186: return 0;
187: }
189: /*TEST
191: build:
192: requires: !mpiuni fftw !complex
194: test:
195: output_file: output/ex144.out
197: test:
198: suffix: 2
199: nsize: 3
200: output_file: output/ex144.out
202: TEST*/