Actual source code: ex146.c

  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>

  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;

 26:   PetscFunctionBeginUser;
 27:   PetscCall(PetscInitialize(&argc, &args, NULL, 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));

 34:   PetscRandomCreate(PETSC_COMM_WORLD, &rnd);

 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);

 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);*/

 44:   /* Allocate space for input and output arrays  */

 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);

 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;

 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));

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

 66:   VecSet(fin, one);
 67:   VecSet(fout, zero);
 68:   VecSet(fout1, zero);

 70:   VecAssemblyBegin(fin);
 71:   VecAssemblyEnd(fin);
 72:   /*    VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */

 74:   VecGetArray(fin, &x_arr);
 75:   VecGetArray(fout1, &z_arr);
 76:   VecGetArray(fout, &y_arr);

 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);

 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:   /*    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);

 99:   if (N2 % 2 == 0) NM = N2 + 2;
100:   else NM = N2 + 1;

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;

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);

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);

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));

150:   VecAssemblyBegin(ini);
151:   VecAssemblyEnd(ini);

153:   VecAssemblyBegin(final);
154:   VecAssemblyEnd(final);

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));

169:   PetscCall(PetscFinalize());
170:   return 0;
171: }