Actual source code: ex147.c

  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>

  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;

 23:   PetscInitialize(&argc, &args, (char *)0, help);
 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:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 28:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

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

 37:   alloc_local = fftw_mpi_local_size_transposed(D, dim, 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);

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

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

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

 67:   VecSet(fin, one);
 68:   /*    VecAssemblyBegin(fin); */
 69:   /*    VecAssemblyEnd(fin); */
 70:   /*    VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */

 72:   VecGetArray(fin, &x_arr);
 73:   VecGetArray(fout1, &z_arr);
 74:   VecGetArray(fout, &y_arr);

 76:   dim[3] = N3;

 78:   fplan = fftw_mpi_plan_dft_r2c(D, dim, (double *)x_arr, (fftw_complex *)y_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);
 79:   bplan = fftw_mpi_plan_dft_c2r(D, dim, (fftw_complex *)y_arr, (double *)z_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);

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

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

 88:   /*    a = 1.0/(PetscReal)N_factor; */
 89:   /*    VecScale(fout1,a); */

 91:   VecAssemblyBegin(fout1);
 92:   VecAssemblyEnd(fout1);

 94:   VecView(fout1, PETSC_VIEWER_STDOUT_WORLD);

 96:   fftw_destroy_plan(fplan);
 97:   fftw_destroy_plan(bplan);
 98:   fftw_free(in1);
 99:   VecDestroy(&fin);
100:   fftw_free(out);
101:   VecDestroy(&fout);
102:   fftw_free(in2);
103:   VecDestroy(&fout1);

105:   PetscFinalize();
106:   return 0;
107: }