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> 6 7 int main(int argc,char **args) 8 { 9 ptrdiff_t N0=256,N1=256,N2=256,N3=2,dim[4]; 10 fftw_plan bplan,fplan; 11 fftw_complex *out; 12 double *in1,*in2; 13 ptrdiff_t alloc_local,local_n0,local_0_start; 14 ptrdiff_t local_n1,local_1_start; 15 PetscInt i,j,indx,n1; 16 PetscInt size,rank,n,N,*in,N_factor,NM; 17 PetscScalar *data_fin,value1,one=1.57,zero=0.0; 18 PetscScalar a,*x_arr,*y_arr,*z_arr,enorm; 19 Vec fin,fout,fout1,ini,final; 20 PetscRandom rnd; 21 PetscErrorCode ierr; 22 VecScatter vecscat,vecscat1; 23 IS indx1,indx2; 24 PetscInt *indx3,k,l,*indx4; 25 PetscInt low,tempindx,tempindx1; 26 27 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 28 #if defined(PETSC_USE_COMPLEX) 29 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex"); 30 #endif 31 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); 32 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 33 34 PetscRandomCreate(PETSC_COMM_WORLD,&rnd); 35 36 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); 37 38 /* printf("The value alloc_local is %ld from process %d\n",alloc_local,rank); */ 39 printf("The value local_n0 is %ld from process %d\n",local_n0,rank); 40 /* printf("The value local_0_start is %ld from process %d\n",local_0_start,rank);*/ 41 /* printf("The value local_n1 is %ld from process %d\n",local_n1,rank); */ 42 /* printf("The value local_1_start is %ld from process %d\n",local_1_start,rank);*/ 43 44 /* Allocate space for input and output arrays */ 45 46 in1=(double*)fftw_malloc(sizeof(double)*alloc_local*2); 47 in2=(double*)fftw_malloc(sizeof(double)*alloc_local*2); 48 out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 49 50 N=2*N0*N1*(N2/2+1);N_factor=N0*N1*N2; 51 n=2*local_n0*N1*(N2/2+1);n1=local_n1*N0*2*N1; 52 53 /* printf("The value N is %d from process %d\n",N,rank); */ 54 /* printf("The value n is %d from process %d\n",n,rank); */ 55 /* printf("The value n1 is %d from process %d\n",n1,rank); */ 56 /* Creating data vector and accompanying array with VeccreateMPIWithArray */ 57 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in1,&fin);CHKERRQ(ierr); 58 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)out,&fout);CHKERRQ(ierr); 59 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in2,&fout1);CHKERRQ(ierr); 60 61 /* VecGetSize(fin,&size); */ 62 /* printf("The size is %d\n",size); */ 63 64 VecSet(fin,one); 65 VecSet(fout,zero); 66 VecSet(fout1,zero); 67 68 VecAssemblyBegin(fin); 69 VecAssemblyEnd(fin); 70 /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */ 71 72 VecGetArray(fin,&x_arr); 73 VecGetArray(fout1,&z_arr); 74 VecGetArray(fout,&y_arr); 75 76 fplan=fftw_mpi_plan_dft_r2c_3d(N0,N1,N2,(double*)x_arr,(fftw_complex*)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE); 77 bplan=fftw_mpi_plan_dft_c2r_3d(N0,N1,N2,(fftw_complex*)y_arr,(double*)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE); 78 79 fftw_execute(fplan); 80 fftw_execute(bplan); 81 82 VecRestoreArray(fin,&x_arr); 83 VecRestoreArray(fout1,&z_arr); 84 VecRestoreArray(fout,&y_arr); 85 86 /* a = 1.0/(PetscReal)N_factor; */ 87 /* ierr = VecScale(fout1,a);CHKERRQ(ierr); */ 88 VecCreate(PETSC_COMM_WORLD,&ini); 89 VecCreate(PETSC_COMM_WORLD,&final); 90 VecSetSizes(ini,local_n0*N1*N2,N_factor); 91 VecSetSizes(final,local_n0*N1*N2,N_factor); 92 /* VecSetSizes(ini,PETSC_DECIDE,N_factor); */ 93 /* VecSetSizes(final,PETSC_DECIDE,N_factor); */ 94 VecSetFromOptions(ini); 95 VecSetFromOptions(final); 96 97 if (N2%2==0) NM=N2+2; 98 else NM=N2+1; 99 100 ierr = VecGetOwnershipRange(fin,&low,NULL);CHKERRQ(ierr); 101 printf("The local index is %d from %d\n",low,rank); 102 ierr = PetscMalloc1(local_n0*N1*N2,&indx3);CHKERRQ(ierr); 103 ierr = PetscMalloc1(local_n0*N1*N2,&indx4);CHKERRQ(ierr); 104 for (i=0; i<local_n0; i++) { 105 for (j=0;j<N1;j++) { 106 for (k=0;k<N2;k++) { 107 tempindx = i*N1*N2 + j*N2 + k; 108 tempindx1 = i*N1*NM + j*NM + k; 109 110 indx3[tempindx]=local_0_start*N1*N2+tempindx; 111 indx4[tempindx]=low+tempindx1; 112 } 113 /* printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */ 114 /* printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */ 115 } 116 } 117 VecGetValues(fin,local_n0*N1*N2,indx4,x_arr); 118 VecSetValues(ini,local_n0*N1*N2,indx3,x_arr,INSERT_VALUES); 119 VecAssemblyBegin(ini); 120 VecAssemblyEnd(ini); 121 122 VecGetValues(fout1,local_n0*N1*N2,indx4,y_arr); 123 VecSetValues(final,local_n0*N1*N2,indx3,y_arr,INSERT_VALUES); 124 VecAssemblyBegin(final); 125 VecAssemblyEnd(final); 126 127 printf("The local index value is %ld from %d",local_n0*N1*N2,rank); 128 /* 129 for (i=0;i<N0;i++) { 130 for (j=0;j<N1;j++) { 131 indx=i*N1*NM+j*NM; 132 ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx1); 133 indx=i*N1*N2+j*N2; 134 ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx2); 135 VecScatterCreate(fin,indx1,ini,indx2,&vecscat); 136 VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD); 137 VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD); 138 VecScatterCreate(fout1,indx1,final,indx2,&vecscat1); 139 VecScatterBegin(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD); 140 VecScatterEnd(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD); 141 } 142 } 143 */ 144 a = 1.0/(PetscReal)N_factor; 145 ierr = VecScale(fout1,a);CHKERRQ(ierr); 146 ierr = VecScale(final,a);CHKERRQ(ierr); 147 148 VecAssemblyBegin(ini); 149 VecAssemblyEnd(ini); 150 151 VecAssemblyBegin(final); 152 VecAssemblyEnd(final); 153 154 /* VecView(final,PETSC_VIEWER_STDOUT_WORLD); */ 155 ierr = VecAXPY(final,-1.0,ini);CHKERRQ(ierr); 156 ierr = VecNorm(final,NORM_1,&enorm);CHKERRQ(ierr); 157 ierr = PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z| = %e\n",enorm);CHKERRQ(ierr); 158 fftw_destroy_plan(fplan); 159 fftw_destroy_plan(bplan); 160 fftw_free(in1); ierr = VecDestroy(&fin);CHKERRQ(ierr); 161 fftw_free(out); ierr = VecDestroy(&fout);CHKERRQ(ierr); 162 fftw_free(in2); ierr = VecDestroy(&fout1);CHKERRQ(ierr); 163 164 ierr = PetscFinalize(); 165 return ierr; 166 } 167