Actual source code: ex146.c

petsc-3.3-p6 2013-02-11
  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:     MPI_Comm_size(PETSC_COMM_WORLD, &size);
 32:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

 34:     PetscRandomCreate(PETSC_COMM_WORLD,&rnd);


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

 39: //    printf("The value alloc_local is %ld from process %d\n",alloc_local,rank);
 40:     printf("The value local_n0 is %ld from process %d\n",local_n0,rank);
 41: //    printf("The value local_0_start is  %ld from process %d\n",local_0_start,rank);
 42: //    printf("The value local_n1 is  %ld from process %d\n",local_n1,rank);
 43: //    printf("The value local_1_start is  %ld from process %d\n",local_1_start,rank);

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

 47:     in1=(double *)fftw_malloc(sizeof(double)*alloc_local*2);
 48:     in2=(double *)fftw_malloc(sizeof(double)*alloc_local*2);
 49:     out=(fftw_complex *)fftw_malloc(sizeof(fftw_complex)*alloc_local);


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

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

 63: //    VecGetSize(fin,&size);
 64: //    printf("The size is %d\n",size);

 66:     VecSet(fin,one);
 67:     VecSet(fout,zero);
 68:     VecSet(fout1,zero);

 70:     VecAssemblyBegin(fin);
 71:     VecAssemblyEnd(fin);
 72: //    VecView(fin,PETSC_VIEWER_STDOUT_WORLD);


 75:     VecGetArray(fin,&x_arr);
 76:     VecGetArray(fout1,&z_arr);
 77:     VecGetArray(fout,&y_arr);
 78: 
 79:     fplan=fftw_mpi_plan_dft_r2c_3d(N0,N1,N2,(double *)x_arr,(fftw_complex *)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
 80:     bplan=fftw_mpi_plan_dft_c2r_3d(N0,N1,N2,(fftw_complex *)y_arr,(double *)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);

 82:     fftw_execute(fplan);
 83:     fftw_execute(bplan);

 85:     VecRestoreArray(fin,&x_arr);
 86:     VecRestoreArray(fout1,&z_arr);
 87:     VecRestoreArray(fout,&y_arr);


 90: //    a = 1.0/(PetscReal)N_factor;
 91: //    VecScale(fout1,a);
 92:     VecCreate(PETSC_COMM_WORLD,&ini);
 93:     VecCreate(PETSC_COMM_WORLD,&final);
 94:     VecSetSizes(ini,local_n0*N1*N2,N_factor);
 95:     VecSetSizes(final,local_n0*N1*N2,N_factor);
 96: //    VecSetSizes(ini,PETSC_DECIDE,N_factor);
 97: //    VecSetSizes(final,PETSC_DECIDE,N_factor);
 98:     VecSetFromOptions(ini);
 99:     VecSetFromOptions(final);
100: 
101:     if (N2%2==0)
102:       NM=N2+2;
103:     else
104:       NM=N2+1;

106:     VecGetOwnershipRange(fin,&low,PETSC_NULL);
107:     printf("The local index is %d from %d\n",low,rank);
108:     PetscMalloc(sizeof(PetscInt)*local_n0*N1*N2,&indx3);
109:     PetscMalloc(sizeof(PetscInt)*local_n0*N1*N2,&indx4);
110:     for (i=0;i<local_n0;i++){
111:        for (j=0;j<N1;j++){
112:           for (k=0;k<N2;k++){
113:             tempindx = i*N1*N2 + j*N2 + k;
114:             tempindx1 = i*N1*NM + j*NM + k;
115:             indx3[tempindx]=local_0_start*N1*N2+tempindx;
116:             indx4[tempindx]=low+tempindx1;
117:           }
118:   //          printf("index3 %d from proc %d is \n",indx3[tempindx],rank);
119:   //          printf("index4 %d from proc %d is \n",indx4[tempindx],rank);
120:        }
121:     }
122:     VecGetValues(fin,local_n0*N1*N2,indx4,x_arr);
123:     VecSetValues(ini,local_n0*N1*N2,indx3,x_arr,INSERT_VALUES);
124:     VecAssemblyBegin(ini);
125:     VecAssemblyEnd(ini);

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

153:     VecAssemblyBegin(ini);
154:     VecAssemblyEnd(ini);
155: 
156:     VecAssemblyBegin(final);
157:     VecAssemblyEnd(final);

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

169:     PetscFinalize();
170:     return 0;
171: }