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