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)";
3: #include <petscmat.h>
4: #include <fftw3.h>
5: #include <fftw3-mpi.h>
9: PetscInt main(PetscInt argc,char **args)
10: {
11: const ptrdiff_t N0=2056,N1=2056;
12: fftw_plan bplan,fplan;
13: fftw_complex *out;
14: double *in1,*in2;
15: ptrdiff_t alloc_local,local_n0,local_0_start;
16: ptrdiff_t local_n1,local_1_start;
17: PetscInt i,j,n1;
18: PetscInt size,rank,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,enorm;
22: Vec fin,fout,fout1;
23: Vec ini,final;
24: PetscRandom rnd;
25: PetscErrorCode ierr;
26: PetscInt *indx3,tempindx,low,*indx4,tempindx1;
27:
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);
48:
49: N=2*N0*(N1/2+1);N_factor=N0*N1;
50: n=2*local_n0*(N1/2+1);n1=local_n1*N0*2;
52: // printf("The value N is %d from process %d\n",N,rank);
53: // printf("The value n is %d from process %d\n",n,rank);
54: // printf("The value n1 is %d from process %d\n",n1,rank);
55: /* Creating data vector and accompanying array with VeccreateMPIWithArray */
56: VecCreateMPIWithArray(PETSC_COMM_WORLD,n,N,( PetscScalar*)in1,&fin);
57: VecCreateMPIWithArray(PETSC_COMM_WORLD,n,N,(PetscScalar*)out,&fout);
58: VecCreateMPIWithArray(PETSC_COMM_WORLD,n,N,(PetscScalar*)in2,&fout1);
60: /* Set the vector with random data */
61: VecSet(fin,zero);
62: // for(i=0;i<N0*N1;i++)
63: // {
64: // VecSetValues(fin,1,&i,&one,INSERT_VALUES);
65: // }
66:
67: // VecSet(fin,one);
68: i=0;
69: VecSetValues(fin,1,&i,&one,INSERT_VALUES);
70: i=1;
71: VecSetValues(fin,1,&i,&two,INSERT_VALUES);
72: i=4;
73: VecSetValues(fin,1,&i,&three,INSERT_VALUES);
74: i=5;
75: VecSetValues(fin,1,&i,&four,INSERT_VALUES);
76: VecAssemblyBegin(fin);
77: VecAssemblyEnd(fin);
78:
79: VecSet(fout,zero);
80: VecSet(fout1,zero);
81:
82: // Get the meaningful portion of array
83: VecGetArray(fin,&x_arr);
84: VecGetArray(fout1,&z_arr);
85: VecGetArray(fout,&y_arr);
87: fplan=fftw_mpi_plan_dft_r2c_2d(N0,N1,(double *)x_arr,(fftw_complex *)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
88: bplan=fftw_mpi_plan_dft_c2r_2d(N0,N1,(fftw_complex *)y_arr,(double *)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
89:
90: fftw_execute(fplan);
91: fftw_execute(bplan);
92:
93: VecRestoreArray(fin,&x_arr);
94: VecRestoreArray(fout1,&z_arr);
95: VecRestoreArray(fout,&y_arr);
97: // VecView(fin,PETSC_VIEWER_STDOUT_WORLD);
98: VecCreate(PETSC_COMM_WORLD,&ini);
99: VecCreate(PETSC_COMM_WORLD,&final);
100: VecSetSizes(ini,local_n0*N1,N0*N1);
101: VecSetSizes(final,local_n0*N1,N0*N1);
102: VecSetFromOptions(ini);
103: VecSetFromOptions(final);
104:
105: if (N1%2==0){
106: NM = N1+2;
107: } else {
108: NM = N1+1;
109: }
110: //printf("The Value of NM is %d",NM);
111: VecGetOwnershipRange(fin,&low,PETSC_NULL);
112: //printf("The local index is %d from %d\n",low,rank);
113: PetscMalloc(sizeof(PetscInt)*local_n0*N1,&indx3);
114: PetscMalloc(sizeof(PetscInt)*local_n0*N1,&indx4);
115: for (i=0;i<local_n0;i++){
116: for (j=0;j<N1;j++){
117: tempindx = i*N1 + j;
118: tempindx1 = i*NM + j;
119: indx3[tempindx]=local_0_start*N1+tempindx;
120: indx4[tempindx]=low+tempindx1;
121: // printf("index3 %d from proc %d is \n",indx3[tempindx],rank);
122: // printf("index4 %d from proc %d is \n",indx4[tempindx],rank);
123: }
124: }
126: VecGetValues(fin,local_n0*N1,indx4,x_arr);
127: VecSetValues(ini,local_n0*N1,indx3,x_arr,INSERT_VALUES);
128: VecAssemblyBegin(ini);
129: VecAssemblyEnd(ini);
131: VecGetValues(fout1,local_n0*N1,indx4,y_arr);
132: VecSetValues(final,local_n0*N1,indx3,y_arr,INSERT_VALUES);
133: VecAssemblyBegin(final);
134: VecAssemblyEnd(final);
136: /*
137: VecScatter vecscat;
138: IS indx1,indx2;
139: for (i=0;i<N0;i++){
140: indx = i*NM;
141: ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx1);
142: indx = i*N1;
143: ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx2);
144: VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
145: VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
146: VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
147: VecScatterBegin(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
148: VecScatterEnd(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
149: }
150: */
152: a = 1.0/(PetscReal)N_factor;
153: VecScale(fout1,a);
154: VecScale(final,a);
155:
157: // VecView(ini,PETSC_VIEWER_STDOUT_WORLD);
158: // VecView(final,PETSC_VIEWER_STDOUT_WORLD);
159: VecAXPY(final,-1.0,ini);
160:
161: VecNorm(final,NORM_1,&enorm);
162: if (enorm > 1.e-10){
163: PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z| = %e\n",enorm);
164: }
165:
166: // Execute fftw with function fftw_execute and destory it after execution
167: fftw_destroy_plan(fplan);
168: fftw_destroy_plan(bplan);
169: fftw_free(in1); VecDestroy(&fin);
170: fftw_free(out); VecDestroy(&fout);
171: fftw_free(in2); VecDestroy(&fout1);
173: VecDestroy(&ini);
174: VecDestroy(&final);
176: PetscRandomDestroy(&rnd);
177: PetscFree(indx3);
178: PetscFree(indx4);
179: PetscFinalize();
180: return 0;
181: }