1 /* This program illustrates use of parallel real FFT */ 2 static char help[]="This program illustrates the use of parallel real multi-dimensional 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=2,N1=2,N2=2,N3=2,dim[4],N,D; 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[100],n1; 16 PetscInt size,rank,n,*in,N_factor; 17 PetscScalar *data_fin,value1,one=1.0,zero=0.0; 18 PetscScalar a,*x_arr,*y_arr,*z_arr,enorm; 19 Vec fin,fout,fout1,x,y; 20 PetscRandom rnd; 21 PetscErrorCode ierr; 22 23 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 24 #if defined(PETSC_USE_COMPLEX) 25 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex"); 26 #endif 27 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); 28 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); 29 30 PetscRandomCreate(PETSC_COMM_WORLD,&rnd); 31 D =4; 32 dim[0]=N0;dim[1]=N1;dim[2]=N2;dim[3]=N3/2+1; 33 34 alloc_local = fftw_mpi_local_size_transposed(D,dim,PETSC_COMM_WORLD,&local_n0,&local_0_start,&local_n1,&local_1_start); 35 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 42 /* Allocate space for input and output arrays */ 43 44 in1=(double*)fftw_malloc(sizeof(double)*alloc_local*2); 45 in2=(double*)fftw_malloc(sizeof(double)*alloc_local*2); 46 out=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local); 47 48 N=2*N0*N1*N2*(N3/2+1);N_factor=N0*N1*N2*N3; 49 n=2*local_n0*N1*N2*(N3/2+1);n1=local_n1*N0*2*N1*N2; 50 51 /* printf("The value N is %d from process %d\n",N,rank); */ 52 /* printf("The value n is %d from process %d\n",n,rank); */ 53 /* printf("The value n1 is %d from process %d\n",n1,rank); */ 54 /* Creating data vector and accompanying array with VeccreateMPIWithArray */ 55 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in1,&fin);CHKERRQ(ierr); 56 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)out,&fout);CHKERRQ(ierr); 57 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,(PetscScalar*)in2,&fout1);CHKERRQ(ierr); 58 59 /* VecGetSize(fin,&size); */ 60 /* printf("The size is %d\n",size); */ 61 62 VecSet(fin,one); 63 /* VecAssemblyBegin(fin); */ 64 /* VecAssemblyEnd(fin); */ 65 /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */ 66 67 VecGetArray(fin,&x_arr); 68 VecGetArray(fout1,&z_arr); 69 VecGetArray(fout,&y_arr); 70 71 dim[3]=N3; 72 73 fplan=fftw_mpi_plan_dft_r2c(D,dim,(double*)x_arr,(fftw_complex*)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE); 74 bplan=fftw_mpi_plan_dft_c2r(D,dim,(fftw_complex*)y_arr,(double*)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE); 75 76 fftw_execute(fplan); 77 fftw_execute(bplan); 78 79 VecRestoreArray(fin,&x_arr); 80 VecRestoreArray(fout1,&z_arr); 81 VecRestoreArray(fout,&y_arr); 82 83 /* a = 1.0/(PetscReal)N_factor; */ 84 /* ierr = VecScale(fout1,a);CHKERRQ(ierr); */ 85 86 VecAssemblyBegin(fout1); 87 VecAssemblyEnd(fout1); 88 89 VecView(fout1,PETSC_VIEWER_STDOUT_WORLD); 90 91 fftw_destroy_plan(fplan); 92 fftw_destroy_plan(bplan); 93 fftw_free(in1); ierr = VecDestroy(&fin);CHKERRQ(ierr); 94 fftw_free(out); ierr = VecDestroy(&fout);CHKERRQ(ierr); 95 fftw_free(in2); ierr = VecDestroy(&fout1);CHKERRQ(ierr); 96 97 ierr = PetscFinalize(); 98 return ierr; 99 } 100