Actual source code: ex146.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: /* This program illustrates use of paralllel real FFT*/
  2: static char help[]="This program illustrates the use of parallel real 3D 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:   ptrdiff_t      N0=256,N1=256,N2=256,N3=2,dim[4];
 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,indx,n1;
 18:   PetscInt       size,rank,n,N,*in,N_factor,NM;
 19:   PetscScalar    *data_fin,value1,one=1.57,zero=0.0;
 20:   PetscScalar    a,*x_arr,*y_arr,*z_arr,enorm;
 21:   Vec            fin,fout,fout1,ini,final;
 22:   PetscRandom    rnd;
 24:   VecScatter     vecscat,vecscat1;
 25:   IS             indx1,indx2;
 26:   PetscInt       *indx3,k,l,*indx4;
 27:   PetscInt       low,tempindx,tempindx1;


 30:   PetscInitialize(&argc,&args,(char*)0,help);
 31: #if defined(PETSC_USE_COMPLEX)
 32:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex");
 33: #endif
 34:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 35:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

 37:   PetscRandomCreate(PETSC_COMM_WORLD,&rnd);


 40:   alloc_local = fftw_mpi_local_size_3d_transposed(N0,N1,N2/2+1,PETSC_COMM_WORLD,&local_n0,&local_0_start,&local_n1,&local_1_start);

 42: /*    printf("The value alloc_local is %ld from process %d\n",alloc_local,rank);     */
 43:   printf("The value local_n0 is %ld from process %d\n",local_n0,rank);
 44: /*    printf("The value local_0_start is  %ld from process %d\n",local_0_start,rank);*/
 45: /*    printf("The value local_n1 is  %ld from process %d\n",local_n1,rank);          */
 46: /*    printf("The value local_1_start is  %ld from process %d\n",local_1_start,rank);*/

 48:   /* Allocate space for input and output arrays  */

 50:   in1=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
 51:   in2=(double*)fftw_malloc(sizeof(double)*alloc_local*2);
 52:   out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);


 55:   N=2*N0*N1*(N2/2+1);N_factor=N0*N1*N2;
 56:   n=2*local_n0*N1*(N2/2+1);n1=local_n1*N0*2*N1;

 58: /*    printf("The value N is  %d from process %d\n",N,rank);   */
 59: /*    printf("The value n is  %d from process %d\n",n,rank);   */
 60: /*    printf("The value n1 is  %d from process %d\n",n1,rank); */
 61:   /* Creating data vector and accompanying array with VeccreateMPIWithArray */
 62:   VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in1,&fin);
 63:   VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)out,&fout);
 64:   VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in2,&fout1);

 66: /*    VecGetSize(fin,&size); */
 67: /*    printf("The size is %d\n",size); */

 69:   VecSet(fin,one);
 70:   VecSet(fout,zero);
 71:   VecSet(fout1,zero);

 73:   VecAssemblyBegin(fin);
 74:   VecAssemblyEnd(fin);
 75: /*    VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */


 78:   VecGetArray(fin,&x_arr);
 79:   VecGetArray(fout1,&z_arr);
 80:   VecGetArray(fout,&y_arr);

 82:   fplan=fftw_mpi_plan_dft_r2c_3d(N0,N1,N2,(double*)x_arr,(fftw_complex*)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
 83:   bplan=fftw_mpi_plan_dft_c2r_3d(N0,N1,N2,(fftw_complex*)y_arr,(double*)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);

 85:   fftw_execute(fplan);
 86:   fftw_execute(bplan);

 88:   VecRestoreArray(fin,&x_arr);
 89:   VecRestoreArray(fout1,&z_arr);
 90:   VecRestoreArray(fout,&y_arr);


 93: /*    a = 1.0/(PetscReal)N_factor; */
 94: /*    VecScale(fout1,a); */
 95:   VecCreate(PETSC_COMM_WORLD,&ini);
 96:   VecCreate(PETSC_COMM_WORLD,&final);
 97:   VecSetSizes(ini,local_n0*N1*N2,N_factor);
 98:   VecSetSizes(final,local_n0*N1*N2,N_factor);
 99: /*    VecSetSizes(ini,PETSC_DECIDE,N_factor); */
100: /*    VecSetSizes(final,PETSC_DECIDE,N_factor); */
101:   VecSetFromOptions(ini);
102:   VecSetFromOptions(final);

104:   if (N2%2==0) NM=N2+2;
105:   else NM=N2+1;

107:   VecGetOwnershipRange(fin,&low,NULL);
108:   printf("The local index is %d from %d\n",low,rank);
109:   PetscMalloc1(local_n0*N1*N2,&indx3);
110:   PetscMalloc1(local_n0*N1*N2,&indx4);
111:   for (i=0; i<local_n0; i++) {
112:     for (j=0;j<N1;j++) {
113:       for (k=0;k<N2;k++) {
114:         tempindx  = i*N1*N2 + j*N2 + k;
115:         tempindx1 = i*N1*NM + j*NM + k;

117:         indx3[tempindx]=local_0_start*N1*N2+tempindx;
118:         indx4[tempindx]=low+tempindx1;
119:       }
120:       /*          printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
121:       /*          printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
122:     }
123:   }
124:   VecGetValues(fin,local_n0*N1*N2,indx4,x_arr);
125:   VecSetValues(ini,local_n0*N1*N2,indx3,x_arr,INSERT_VALUES);
126:   VecAssemblyBegin(ini);
127:   VecAssemblyEnd(ini);

129:   VecGetValues(fout1,local_n0*N1*N2,indx4,y_arr);
130:   VecSetValues(final,local_n0*N1*N2,indx3,y_arr,INSERT_VALUES);
131:   VecAssemblyBegin(final);
132:   VecAssemblyEnd(final);

134:   printf("The local index value is %ld from %d",local_n0*N1*N2,rank);
135: /*
136:   for (i=0;i<N0;i++) {
137:      for (j=0;j<N1;j++) {
138:         indx=i*N1*NM+j*NM;
139:         ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx1);
140:         indx=i*N1*N2+j*N2;
141:         ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx2);
142:         VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
143:         VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
144:         VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
145:         VecScatterCreate(fout1,indx1,final,indx2,&vecscat1);
146:         VecScatterBegin(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
147:         VecScatterEnd(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
148:      }
149:   }
150: */
151:   a    = 1.0/(PetscReal)N_factor;
152:   VecScale(fout1,a);
153:   VecScale(final,a);

155:   VecAssemblyBegin(ini);
156:   VecAssemblyEnd(ini);

158:   VecAssemblyBegin(final);
159:   VecAssemblyEnd(final);

161: /*    VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
162:   VecAXPY(final,-1.0,ini);
163:   VecNorm(final,NORM_1,&enorm);
164:   PetscPrintf(PETSC_COMM_WORLD,"  Error norm of |x - z|  = %e\n",enorm);
165:   fftw_destroy_plan(fplan);
166:   fftw_destroy_plan(bplan);
167:   fftw_free(in1); VecDestroy(&fin);
168:   fftw_free(out); VecDestroy(&fout);
169:   fftw_free(in2); VecDestroy(&fout1);

171:   PetscFinalize();
172:   return 0;
173: }