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 7d71ae5a4SJacob Faibussowitsch int main(int argc, char **args) 8d71ae5a4SJacob Faibussowitsch { 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 22327415f7SBarry Smith PetscFunctionBeginUser; 23*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &args, NULL, help)); 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 279566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 289566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 29c4762a1bSJed Brown 30c4762a1bSJed Brown PetscRandomCreate(PETSC_COMM_WORLD, &rnd); 31c4762a1bSJed Brown D = 4; 329371c9d4SSatish Balay dim[0] = N0; 339371c9d4SSatish Balay dim[1] = N1; 349371c9d4SSatish Balay dim[2] = N2; 359371c9d4SSatish Balay dim[3] = N3 / 2 + 1; 36c4762a1bSJed Brown 37c4762a1bSJed Brown alloc_local = fftw_mpi_local_size_transposed(D, dim, PETSC_COMM_WORLD, &local_n0, &local_0_start, &local_n1, &local_1_start); 38c4762a1bSJed Brown 39c4762a1bSJed Brown printf("The value alloc_local is %ld from process %d\n", alloc_local, rank); 40c4762a1bSJed Brown printf("The value local_n0 is %ld from process %d\n", local_n0, rank); 41c4762a1bSJed Brown printf("The value local_0_start is %ld from process %d\n", local_0_start, rank); 42c4762a1bSJed Brown printf("The value local_n1 is %ld from process %d\n", local_n1, rank); 43c4762a1bSJed Brown printf("The value local_1_start is %ld from process %d\n", local_1_start, rank); 44c4762a1bSJed Brown 45c4762a1bSJed Brown /* Allocate space for input and output arrays */ 46c4762a1bSJed Brown 47c4762a1bSJed Brown in1 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2); 48c4762a1bSJed Brown in2 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2); 49c4762a1bSJed Brown out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local); 50c4762a1bSJed Brown 519371c9d4SSatish Balay N = 2 * N0 * N1 * N2 * (N3 / 2 + 1); 529371c9d4SSatish Balay N_factor = N0 * N1 * N2 * N3; 539371c9d4SSatish Balay n = 2 * local_n0 * N1 * N2 * (N3 / 2 + 1); 549371c9d4SSatish Balay n1 = local_n1 * N0 * 2 * N1 * N2; 55c4762a1bSJed Brown 56c4762a1bSJed Brown /* printf("The value N is %d from process %d\n",N,rank); */ 57c4762a1bSJed Brown /* printf("The value n is %d from process %d\n",n,rank); */ 58c4762a1bSJed Brown /* printf("The value n1 is %d from process %d\n",n1,rank); */ 59c4762a1bSJed Brown /* Creating data vector and accompanying array with VeccreateMPIWithArray */ 609566063dSJacob Faibussowitsch PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in1, &fin)); 619566063dSJacob Faibussowitsch PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)out, &fout)); 629566063dSJacob Faibussowitsch PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in2, &fout1)); 63c4762a1bSJed Brown 64c4762a1bSJed Brown /* VecGetSize(fin,&size); */ 65c4762a1bSJed Brown /* printf("The size is %d\n",size); */ 66c4762a1bSJed Brown 67c4762a1bSJed Brown VecSet(fin, one); 68c4762a1bSJed Brown /* VecAssemblyBegin(fin); */ 69c4762a1bSJed Brown /* VecAssemblyEnd(fin); */ 70c4762a1bSJed Brown /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */ 71c4762a1bSJed Brown 72c4762a1bSJed Brown VecGetArray(fin, &x_arr); 73c4762a1bSJed Brown VecGetArray(fout1, &z_arr); 74c4762a1bSJed Brown VecGetArray(fout, &y_arr); 75c4762a1bSJed Brown 76c4762a1bSJed Brown dim[3] = N3; 77c4762a1bSJed Brown 78c4762a1bSJed Brown fplan = fftw_mpi_plan_dft_r2c(D, dim, (double *)x_arr, (fftw_complex *)y_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE); 79c4762a1bSJed Brown bplan = fftw_mpi_plan_dft_c2r(D, dim, (fftw_complex *)y_arr, (double *)z_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE); 80c4762a1bSJed Brown 81c4762a1bSJed Brown fftw_execute(fplan); 82c4762a1bSJed Brown fftw_execute(bplan); 83c4762a1bSJed Brown 84c4762a1bSJed Brown VecRestoreArray(fin, &x_arr); 85c4762a1bSJed Brown VecRestoreArray(fout1, &z_arr); 86c4762a1bSJed Brown VecRestoreArray(fout, &y_arr); 87c4762a1bSJed Brown 88c4762a1bSJed Brown /* a = 1.0/(PetscReal)N_factor; */ 899566063dSJacob Faibussowitsch /* PetscCall(VecScale(fout1,a)); */ 90c4762a1bSJed Brown 91c4762a1bSJed Brown VecAssemblyBegin(fout1); 92c4762a1bSJed Brown VecAssemblyEnd(fout1); 93c4762a1bSJed Brown 94c4762a1bSJed Brown VecView(fout1, PETSC_VIEWER_STDOUT_WORLD); 95c4762a1bSJed Brown 96c4762a1bSJed Brown fftw_destroy_plan(fplan); 97c4762a1bSJed Brown fftw_destroy_plan(bplan); 989371c9d4SSatish Balay fftw_free(in1); 999371c9d4SSatish Balay PetscCall(VecDestroy(&fin)); 1009371c9d4SSatish Balay fftw_free(out); 1019371c9d4SSatish Balay PetscCall(VecDestroy(&fout)); 1029371c9d4SSatish Balay fftw_free(in2); 1039371c9d4SSatish Balay PetscCall(VecDestroy(&fout1)); 104c4762a1bSJed Brown 1059566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 106b122ec5aSJacob Faibussowitsch return 0; 107c4762a1bSJed Brown } 108