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: }