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