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;
28: PetscInitialize(&argc,&args,(char*)0,help);
29: MPI_Comm_size(PETSC_COMM_WORLD, &size);
30: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
32: PetscRandomCreate(PETSC_COMM_WORLD,&rnd);
34: 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);
35: #if defined(DEBUGGING)
36: printf("The value alloc_local is %ld from process %d\n",alloc_local,rank);
37: printf("The value local_n0 is %ld from process %d\n",local_n0,rank);
38: printf("The value local_0_start is %ld from process %d\n",local_0_start,rank);
39: /* printf("The value local_n1 is %ld from process %d\n",local_n1,rank); */
40: /* printf("The value local_1_start is %ld from process %d\n",local_1_start,rank); */
41: /* printf("The value local_n0 is %ld from process %d\n",local_n0,rank); */
42: #endif
44: /* Allocate space for input and output arrays */
45: in1=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
46: in2=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
47: out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
49: N = 2*N0*(N1/2+1);
50: N_factor = N0*N1;
51: n = 2*local_n0*(N1/2+1);
53: /* printf("The value N is %d from process %d\n",N,rank); */
54: /* printf("The value n is %d from process %d\n",n,rank); */
55: /* printf("The value n1 is %d from process %d\n",n1,rank);*/
56: /* Creating data vector and accompanying array with VeccreateMPIWithArray */
57: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in1,&fin);
58: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)out,&fout);
59: VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in2,&fout1);
61: /* Set the vector with random data */
62: VecSet(fin,zero);
63: /* for (i=0;i<N0*N1;i++) */
64: /* { */
65: /* VecSetValues(fin,1,&i,&one,INSERT_VALUES); */
66: /* } */
68: /* VecSet(fin,one); */
69: i =0;
70: VecSetValues(fin,1,&i,&one,INSERT_VALUES);
71: i =1;
72: VecSetValues(fin,1,&i,&two,INSERT_VALUES);
73: i =4;
74: VecSetValues(fin,1,&i,&three,INSERT_VALUES);
75: i =5;
76: VecSetValues(fin,1,&i,&four,INSERT_VALUES);
77: VecAssemblyBegin(fin);
78: VecAssemblyEnd(fin);
80: VecSet(fout,zero);
81: VecSet(fout1,zero);
83: /* Get the meaningful portion of array */
84: VecGetArray(fin,&x_arr);
85: VecGetArray(fout1,&z_arr);
86: VecGetArray(fout,&y_arr);
88: fplan=fftw_mpi_plan_dft_r2c_2d(N0,N1,(double*)x_arr,(fftw_complex*)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
89: bplan=fftw_mpi_plan_dft_c2r_2d(N0,N1,(fftw_complex*)y_arr,(double*)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
91: fftw_execute(fplan);
92: fftw_execute(bplan);
94: VecRestoreArray(fin,&x_arr);
95: VecRestoreArray(fout1,&z_arr);
96: VecRestoreArray(fout,&y_arr);
98: /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */
99: VecCreate(PETSC_COMM_WORLD,&ini);
100: VecCreate(PETSC_COMM_WORLD,&final);
101: VecSetSizes(ini,local_n0*N1,N0*N1);
102: VecSetSizes(final,local_n0*N1,N0*N1);
103: VecSetFromOptions(ini);
104: VecSetFromOptions(final);
106: if (N1%2==0) {
107: NM = N1+2;
108: } else {
109: NM = N1+1;
110: }
111: /*printf("The Value of NM is %d",NM); */
112: VecGetOwnershipRange(fin,&low,NULL);
113: /*printf("The local index is %d from %d\n",low,rank); */
114: PetscMalloc1(local_n0*N1,&indx3);
115: PetscMalloc1(local_n0*N1,&indx4);
116: for (i=0;i<local_n0;i++) {
117: for (j=0;j<N1;j++) {
118: tempindx = i*N1 + j;
119: tempindx1 = i*NM + j;
121: indx3[tempindx]=local_0_start*N1+tempindx;
122: indx4[tempindx]=low+tempindx1;
123: /* printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
124: /* printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
125: }
126: }
128: PetscMalloc2(local_n0*N1,&x_arr,local_n0*N1,&y_arr); /* arr must be allocated for VecGetValues() */
129: VecGetValues(fin,local_n0*N1,indx4,(PetscScalar*)x_arr);
130: VecSetValues(ini,local_n0*N1,indx3,x_arr,INSERT_VALUES);
132: VecAssemblyBegin(ini);
133: VecAssemblyEnd(ini);
135: VecGetValues(fout1,local_n0*N1,indx4,y_arr);
136: VecSetValues(final,local_n0*N1,indx3,y_arr,INSERT_VALUES);
137: VecAssemblyBegin(final);
138: VecAssemblyEnd(final);
139: PetscFree2(x_arr,y_arr);
141: /*
142: VecScatter vecscat;
143: IS indx1,indx2;
144: for (i=0;i<N0;i++) {
145: indx = i*NM;
146: ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx1);
147: indx = i*N1;
148: ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx2);
149: VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
150: VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
151: VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
152: VecScatterBegin(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
153: VecScatterEnd(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
154: }
155: */
157: a = 1.0/(PetscReal)N_factor;
158: VecScale(fout1,a);
159: VecScale(final,a);
161: /* VecView(ini,PETSC_VIEWER_STDOUT_WORLD); */
162: /* VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
163: VecAXPY(final,-1.0,ini);
165: VecNorm(final,NORM_1,&enorm);
166: if (enorm > 1.e-10) {
167: PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z| = %e\n",enorm);
168: }
170: /* Execute fftw with function fftw_execute and destroy it after execution */
171: fftw_destroy_plan(fplan);
172: fftw_destroy_plan(bplan);
173: fftw_free(in1); VecDestroy(&fin);
174: fftw_free(out); VecDestroy(&fout);
175: fftw_free(in2); VecDestroy(&fout1);
177: VecDestroy(&ini);
178: VecDestroy(&final);
180: PetscRandomDestroy(&rnd);
181: PetscFree(indx3);
182: PetscFree(indx4);
183: PetscFinalize();
184: return 0;
185: }
187: /*TEST
189: build:
190: requires: !mpiuni fftw !complex
192: test:
193: output_file: output/ex144.out
195: test:
196: suffix: 2
197: nsize: 3
198: output_file: output/ex144.out
200: TEST*/