1 /* This program illustrates use of parallel real FFT*/ 2 static char help[] = "This program illustrates the use of parallel real 3D 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 = 256, N1 = 256, N2 = 256, N3 = 2, dim[4]; 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, n1; 16 PetscInt size, rank, n, N, *in, N_factor, NM; 17 PetscScalar *data_fin, value1, one = 1.57, zero = 0.0; 18 PetscScalar a, *x_arr, *y_arr, *z_arr, enorm; 19 Vec fin, fout, fout1, ini, final; 20 PetscRandom rnd; 21 VecScatter vecscat, vecscat1; 22 IS indx1, indx2; 23 PetscInt *indx3, k, l, *indx4; 24 PetscInt low, tempindx, tempindx1; 25 26 PetscFunctionBeginUser; 27 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 28 #if defined(PETSC_USE_COMPLEX) 29 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers. Your current scalar type is complex"); 30 #endif 31 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 32 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 33 34 PetscRandomCreate(PETSC_COMM_WORLD, &rnd); 35 36 alloc_local = fftw_mpi_local_size_3d_transposed(N0, N1, N2 / 2 + 1, PETSC_COMM_WORLD, &local_n0, &local_0_start, &local_n1, &local_1_start); 37 38 /* printf("The value alloc_local is %ld from process %d\n",alloc_local,rank); */ 39 printf("The value local_n0 is %ld from process %d\n", local_n0, rank); 40 /* printf("The value local_0_start is %ld from process %d\n",local_0_start,rank);*/ 41 /* printf("The value local_n1 is %ld from process %d\n",local_n1,rank); */ 42 /* printf("The value local_1_start is %ld from process %d\n",local_1_start,rank);*/ 43 44 /* Allocate space for input and output arrays */ 45 46 in1 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2); 47 in2 = (double *)fftw_malloc(sizeof(double) * alloc_local * 2); 48 out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * alloc_local); 49 50 N = 2 * N0 * N1 * (N2 / 2 + 1); 51 N_factor = N0 * N1 * N2; 52 n = 2 * local_n0 * N1 * (N2 / 2 + 1); 53 n1 = local_n1 * N0 * 2 * N1; 54 55 /* printf("The value N is %d from process %d\n",N,rank); */ 56 /* printf("The value n is %d from process %d\n",n,rank); */ 57 /* printf("The value n1 is %d from process %d\n",n1,rank); */ 58 /* Creating data vector and accompanying array with VeccreateMPIWithArray */ 59 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in1, &fin)); 60 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)out, &fout)); 61 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in2, &fout1)); 62 63 /* VecGetSize(fin,&size); */ 64 /* printf("The size is %d\n",size); */ 65 66 VecSet(fin, one); 67 VecSet(fout, zero); 68 VecSet(fout1, zero); 69 70 VecAssemblyBegin(fin); 71 VecAssemblyEnd(fin); 72 /* VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */ 73 74 VecGetArray(fin, &x_arr); 75 VecGetArray(fout1, &z_arr); 76 VecGetArray(fout, &y_arr); 77 78 fplan = fftw_mpi_plan_dft_r2c_3d(N0, N1, N2, (double *)x_arr, (fftw_complex *)y_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE); 79 bplan = fftw_mpi_plan_dft_c2r_3d(N0, N1, N2, (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 VecCreate(PETSC_COMM_WORLD, &ini); 91 VecCreate(PETSC_COMM_WORLD, &final); 92 VecSetSizes(ini, local_n0 * N1 * N2, N_factor); 93 VecSetSizes(final, local_n0 * N1 * N2, N_factor); 94 /* VecSetSizes(ini,PETSC_DECIDE,N_factor); */ 95 /* VecSetSizes(final,PETSC_DECIDE,N_factor); */ 96 VecSetFromOptions(ini); 97 VecSetFromOptions(final); 98 99 if (N2 % 2 == 0) NM = N2 + 2; 100 else NM = N2 + 1; 101 102 PetscCall(VecGetOwnershipRange(fin, &low, NULL)); 103 printf("The local index is %d from %d\n", low, rank); 104 PetscCall(PetscMalloc1(local_n0 * N1 * N2, &indx3)); 105 PetscCall(PetscMalloc1(local_n0 * N1 * N2, &indx4)); 106 for (i = 0; i < local_n0; i++) { 107 for (j = 0; j < N1; j++) { 108 for (k = 0; k < N2; k++) { 109 tempindx = i * N1 * N2 + j * N2 + k; 110 tempindx1 = i * N1 * NM + j * NM + k; 111 112 indx3[tempindx] = local_0_start * N1 * N2 + tempindx; 113 indx4[tempindx] = low + tempindx1; 114 } 115 /* printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */ 116 /* printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */ 117 } 118 } 119 VecGetValues(fin, local_n0 * N1 * N2, indx4, x_arr); 120 VecSetValues(ini, local_n0 * N1 * N2, indx3, x_arr, INSERT_VALUES); 121 VecAssemblyBegin(ini); 122 VecAssemblyEnd(ini); 123 124 VecGetValues(fout1, local_n0 * N1 * N2, indx4, y_arr); 125 VecSetValues(final, local_n0 * N1 * N2, indx3, y_arr, INSERT_VALUES); 126 VecAssemblyBegin(final); 127 VecAssemblyEnd(final); 128 129 printf("The local index value is %ld from %d", local_n0 * N1 * N2, rank); 130 /* 131 for (i=0;i<N0;i++) { 132 for (j=0;j<N1;j++) { 133 indx=i*N1*NM+j*NM; 134 ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx1); 135 indx=i*N1*N2+j*N2; 136 ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx2); 137 VecScatterCreate(fin,indx1,ini,indx2,&vecscat); 138 VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD); 139 VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD); 140 VecScatterCreate(fout1,indx1,final,indx2,&vecscat1); 141 VecScatterBegin(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD); 142 VecScatterEnd(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD); 143 } 144 } 145 */ 146 a = 1.0 / (PetscReal)N_factor; 147 PetscCall(VecScale(fout1, a)); 148 PetscCall(VecScale(final, a)); 149 150 VecAssemblyBegin(ini); 151 VecAssemblyEnd(ini); 152 153 VecAssemblyBegin(final); 154 VecAssemblyEnd(final); 155 156 /* VecView(final,PETSC_VIEWER_STDOUT_WORLD); */ 157 PetscCall(VecAXPY(final, -1.0, ini)); 158 PetscCall(VecNorm(final, NORM_1, &enorm)); 159 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z| = %e\n", enorm)); 160 fftw_destroy_plan(fplan); 161 fftw_destroy_plan(bplan); 162 fftw_free(in1); 163 PetscCall(VecDestroy(&fin)); 164 fftw_free(out); 165 PetscCall(VecDestroy(&fout)); 166 fftw_free(in2); 167 PetscCall(VecDestroy(&fout1)); 168 169 PetscCall(PetscFinalize()); 170 return 0; 171 } 172