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 22 PetscFunctionBeginUser; 23 PetscCall(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 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 28 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 29 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; 36 37 alloc_local = fftw_mpi_local_size_transposed(D, dim, PETSC_COMM_WORLD, &local_n0, &local_0_start, &local_n1, &local_1_start); 38 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); 44 45 /* Allocate space for input and output arrays */ 46 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); 50 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; 55 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 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in1, &fin)); 61 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)out, &fout)); 62 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in2, &fout1)); 63 64 /* VecGetSize(fin,&size); */ 65 /* printf("The size is %d\n",size); */ 66 67 VecSet(fin, one); 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 dim[3] = N3; 77 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); 80 81 fftw_execute(fplan); 82 fftw_execute(bplan); 83 84 VecRestoreArray(fin, &x_arr); 85 VecRestoreArray(fout1, &z_arr); 86 VecRestoreArray(fout, &y_arr); 87 88 /* a = 1.0/(PetscReal)N_factor; */ 89 /* PetscCall(VecScale(fout1,a)); */ 90 91 VecAssemblyBegin(fout1); 92 VecAssemblyEnd(fout1); 93 94 VecView(fout1, PETSC_VIEWER_STDOUT_WORLD); 95 96 fftw_destroy_plan(fplan); 97 fftw_destroy_plan(bplan); 98 fftw_free(in1); 99 PetscCall(VecDestroy(&fin)); 100 fftw_free(out); 101 PetscCall(VecDestroy(&fout)); 102 fftw_free(in2); 103 PetscCall(VecDestroy(&fout1)); 104 105 PetscCall(PetscFinalize()); 106 return 0; 107 } 108