Actual source code: ex144.c

petsc-3.4.2 2013-07-02
  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;

 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);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,1,n,N,(PetscScalar*)in1,&fin);
 57:   VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)out,&fout);
 58:   VecCreateMPIWithArray(PETSC_COMM_WORLD,1,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: /*     } */

 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);

 79:   VecSet(fout,zero);
 80:   VecSet(fout1,zero);

 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);

 90:   fftw_execute(fplan);
 91:   fftw_execute(bplan);

 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);

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,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;

120:       indx3[tempindx]=local_0_start*N1+tempindx;
121:       indx4[tempindx]=low+tempindx1;
122:       /*          printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
123:       /*          printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
124:     }
125:   }

127:   VecGetValues(fin,local_n0*N1,indx4,x_arr);
128:   VecSetValues(ini,local_n0*N1,indx3,x_arr,INSERT_VALUES);
129:   VecAssemblyBegin(ini);
130:   VecAssemblyEnd(ini);

132:   VecGetValues(fout1,local_n0*N1,indx4,y_arr);
133:   VecSetValues(final,local_n0*N1,indx3,y_arr,INSERT_VALUES);
134:   VecAssemblyBegin(final);
135:   VecAssemblyEnd(final);

137: /*
138:     VecScatter      vecscat;
139:     IS              indx1,indx2;
140:     for (i=0;i<N0;i++) {
141:        indx = i*NM;
142:        ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx1);
143:        indx = i*N1;
144:        ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx2);
145:        VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
146:        VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
147:        VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
148:        VecScatterBegin(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
149:        VecScatterEnd(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
150:     }
151: */

153:   a    = 1.0/(PetscReal)N_factor;
154:   VecScale(fout1,a);
155:   VecScale(final,a);


158: /*    VecView(ini,PETSC_VIEWER_STDOUT_WORLD);   */
159: /*    VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
160:   VecAXPY(final,-1.0,ini);

162:   VecNorm(final,NORM_1,&enorm);
163:   if (enorm > 1.e-10) {
164:     PetscPrintf(PETSC_COMM_WORLD,"  Error norm of |x - z|  = %e\n",enorm);
165:   }

167:   /* Execute fftw with function fftw_execute and destory it after execution */
168:   fftw_destroy_plan(fplan);
169:   fftw_destroy_plan(bplan);
170:   fftw_free(in1);  VecDestroy(&fin);
171:   fftw_free(out);  VecDestroy(&fout);
172:   fftw_free(in2);  VecDestroy(&fout1);

174:   VecDestroy(&ini);
175:   VecDestroy(&final);

177:   PetscRandomDestroy(&rnd);
178:   PetscFree(indx3);
179:   PetscFree(indx4);
180:   PetscFinalize();
181:   return 0;
182: }