Actual source code: ex147.c

petsc-3.3-p6 2013-02-11
  1: /* This program illustrates use of parallel real FFT */
  2: static char help[]="This program illustrates the use of parallel real 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=2,N1=2,N2=2,N3=2,dim[4],N,D;
 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[100],n1;
 18:     PetscInt  size,rank,n,*in,N_factor;
 19:     PetscScalar *data_fin,value1,one=1.0,zero=0.0;
 20:     PetscScalar a,*x_arr,*y_arr,*z_arr,enorm;
 21:     Vec fin,fout,fout1,x,y;
 22:     PetscRandom    rnd;

 25:     PetscInitialize(&argc,&args,(char *)0,help);
 26:     MPI_Comm_size(PETSC_COMM_WORLD, &size);
 27:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

 29:     PetscRandomCreate(PETSC_COMM_WORLD,&rnd);
 30:     D=4;
 31:     dim[0]=N0;dim[1]=N1;dim[2]=N2;dim[3]=N3/2+1;


 34:     alloc_local = fftw_mpi_local_size_transposed(D,dim,PETSC_COMM_WORLD,&local_n0,&local_0_start,&local_n1,&local_1_start);

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

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

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


 49:     N=2*N0*N1*N2*(N3/2+1);N_factor=N0*N1*N2*N3;
 50:     n=2*local_n0*N1*N2*(N3/2+1);n1=local_n1*N0*2*N1*N2;

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

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

 63:     VecSet(fin,one);
 64: //    VecAssemblyBegin(fin);
 65: //    VecAssemblyEnd(fin);
 66: //    VecView(fin,PETSC_VIEWER_STDOUT_WORLD);


 69:     VecGetArray(fin,&x_arr);
 70:     VecGetArray(fout1,&z_arr);
 71:     VecGetArray(fout,&y_arr);
 72: 
 73:     dim[3]=N3;
 74: 
 75:     fplan=fftw_mpi_plan_dft_r2c(D,dim,(double *)x_arr,(fftw_complex *)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
 76:     bplan=fftw_mpi_plan_dft_c2r(D,dim,(fftw_complex *)y_arr,(double *)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);

 78:     fftw_execute(fplan);
 79:     fftw_execute(bplan);

 81:     VecRestoreArray(fin,&x_arr);
 82:     VecRestoreArray(fout1,&z_arr);
 83:     VecRestoreArray(fout,&y_arr);

 85: //    a = 1.0/(PetscReal)N_factor;
 86: //    VecScale(fout1,a);

 88:     VecAssemblyBegin(fout1);
 89:     VecAssemblyEnd(fout1);

 91:     VecView(fout1,PETSC_VIEWER_STDOUT_WORLD);

 93:     fftw_destroy_plan(fplan);
 94:     fftw_destroy_plan(bplan);
 95:     fftw_free(in1); VecDestroy(&fin) ;
 96:     fftw_free(out); VecDestroy(&fout);
 97:     fftw_free(in2); VecDestroy(&fout1);

 99:     PetscFinalize();
100:     return 0;
101: }