Actual source code: ex144.c

  1: /* This program illustrates use of parallel real FFT */
  2: static char help[] = "This program illustrates the use of parallel real 2D fft using fftw without PETSc interface";

  4: #include <petscmat.h>
  5: #include <fftw3.h>
  6: #include <fftw3-mpi.h>

  8: int main(int argc, char **args)
  9: {
 10:   const ptrdiff_t N0 = 2056, N1 = 2056;
 11:   fftw_plan       bplan, fplan;
 12:   fftw_complex   *out;
 13:   double         *in1, *in2;
 14:   ptrdiff_t       alloc_local, local_n0, local_0_start;
 15:   ptrdiff_t       local_n1, local_1_start;
 16:   PetscInt        i, j;
 17:   PetscMPIInt     size, rank;
 18:   int             n, N, N_factor, NM;
 19:   PetscScalar     one = 2.0, zero = 0.5;
 20:   PetscScalar     two = 4.0, three = 8.0, four = 16.0;
 21:   PetscScalar     a, *x_arr, *y_arr, *z_arr;
 22:   PetscReal       enorm;
 23:   Vec             fin, fout, fout1;
 24:   Vec             ini, final;
 25:   PetscRandom     rnd;
 26:   PetscInt       *indx3, tempindx, low, *indx4, tempindx1;

 29:   PetscInitialize(&argc, &args, (char *)0, help);
 30:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 31:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

 33:   PetscRandomCreate(PETSC_COMM_WORLD, &rnd);

 35:   alloc_local = fftw_mpi_local_size_2d_transposed(N0, N1 / 2 + 1, PETSC_COMM_WORLD, &local_n0, &local_0_start, &local_n1, &local_1_start);
 36: #if defined(DEBUGGING)
 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: /*    printf("The value local_n0 is  %ld from process %d\n",local_n0,rank); */
 43: #endif

 45:   /* 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 / 2 + 1);
 51:   N_factor = N0 * N1;
 52:   n        = 2 * local_n0 * (N1 / 2 + 1);

 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:   VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in1, &fin);
 59:   VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)out, &fout);
 60:   VecCreateMPIWithArray(PETSC_COMM_WORLD, 1, n, N, (PetscScalar *)in2, &fout1);

 62:   /* Set the vector with random data */
 63:   VecSet(fin, zero);
 64:   /*    for (i=0;i<N0*N1;i++) */
 65:   /*       { */
 66:   /*       VecSetValues(fin,1,&i,&one,INSERT_VALUES); */
 67:   /*     } */

 69:   /*    VecSet(fin,one); */
 70:   i = 0;
 71:   VecSetValues(fin, 1, &i, &one, INSERT_VALUES);
 72:   i = 1;
 73:   VecSetValues(fin, 1, &i, &two, INSERT_VALUES);
 74:   i = 4;
 75:   VecSetValues(fin, 1, &i, &three, INSERT_VALUES);
 76:   i = 5;
 77:   VecSetValues(fin, 1, &i, &four, INSERT_VALUES);
 78:   VecAssemblyBegin(fin);
 79:   VecAssemblyEnd(fin);

 81:   VecSet(fout, zero);
 82:   VecSet(fout1, zero);

 84:   /* Get the meaningful portion of array */
 85:   VecGetArray(fin, &x_arr);
 86:   VecGetArray(fout1, &z_arr);
 87:   VecGetArray(fout, &y_arr);

 89:   fplan = fftw_mpi_plan_dft_r2c_2d(N0, N1, (double *)x_arr, (fftw_complex *)y_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);
 90:   bplan = fftw_mpi_plan_dft_c2r_2d(N0, N1, (fftw_complex *)y_arr, (double *)z_arr, PETSC_COMM_WORLD, FFTW_ESTIMATE);

 92:   fftw_execute(fplan);
 93:   fftw_execute(bplan);

 95:   VecRestoreArray(fin, &x_arr);
 96:   VecRestoreArray(fout1, &z_arr);
 97:   VecRestoreArray(fout, &y_arr);

 99:   /*    VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */
100:   VecCreate(PETSC_COMM_WORLD, &ini);
101:   VecCreate(PETSC_COMM_WORLD, &final);
102:   VecSetSizes(ini, local_n0 * N1, N0 * N1);
103:   VecSetSizes(final, local_n0 * N1, N0 * N1);
104:   VecSetFromOptions(ini);
105:   VecSetFromOptions(final);

107:   if (N1 % 2 == 0) {
108:     NM = N1 + 2;
109:   } else {
110:     NM = N1 + 1;
111:   }
112:   /*printf("The Value of NM is %d",NM); */
113:   VecGetOwnershipRange(fin, &low, NULL);
114:   /*printf("The local index is %d from %d\n",low,rank); */
115:   PetscMalloc1(local_n0 * N1, &indx3);
116:   PetscMalloc1(local_n0 * N1, &indx4);
117:   for (i = 0; i < local_n0; i++) {
118:     for (j = 0; j < N1; j++) {
119:       tempindx  = i * N1 + j;
120:       tempindx1 = i * NM + j;

122:       indx3[tempindx] = local_0_start * N1 + tempindx;
123:       indx4[tempindx] = low + tempindx1;
124:       /*          printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
125:       /*          printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
126:     }
127:   }

129:   PetscMalloc2(local_n0 * N1, &x_arr, local_n0 * N1, &y_arr); /* arr must be allocated for VecGetValues() */
130:   VecGetValues(fin, local_n0 * N1, indx4, (PetscScalar *)x_arr);
131:   VecSetValues(ini, local_n0 * N1, indx3, x_arr, INSERT_VALUES);

133:   VecAssemblyBegin(ini);
134:   VecAssemblyEnd(ini);

136:   VecGetValues(fout1, local_n0 * N1, indx4, y_arr);
137:   VecSetValues(final, local_n0 * N1, indx3, y_arr, INSERT_VALUES);
138:   VecAssemblyBegin(final);
139:   VecAssemblyEnd(final);
140:   PetscFree2(x_arr, y_arr);

142:   /*
143:     VecScatter      vecscat;
144:     IS              indx1,indx2;
145:     for (i=0;i<N0;i++) {
146:        indx = i*NM;
147:        ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx1);
148:        indx = i*N1;
149:        ISCreateStride(PETSC_COMM_WORLD,N1,indx,1,&indx2);
150:        VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
151:        VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
152:        VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
153:        VecScatterBegin(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
154:        VecScatterEnd(vecscat,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
155:     }
156: */

158:   a = 1.0 / (PetscReal)N_factor;
159:   VecScale(fout1, a);
160:   VecScale(final, a);

162:   /*    VecView(ini,PETSC_VIEWER_STDOUT_WORLD);   */
163:   /*    VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
164:   VecAXPY(final, -1.0, ini);

166:   VecNorm(final, NORM_1, &enorm);
167:   if (enorm > 1.e-10) PetscPrintf(PETSC_COMM_WORLD, "  Error norm of |x - z|  = %e\n", enorm);

169:   /* Execute fftw with function fftw_execute and destroy it after execution */
170:   fftw_destroy_plan(fplan);
171:   fftw_destroy_plan(bplan);
172:   fftw_free(in1);
173:   VecDestroy(&fin);
174:   fftw_free(out);
175:   VecDestroy(&fout);
176:   fftw_free(in2);
177:   VecDestroy(&fout1);

179:   VecDestroy(&ini);
180:   VecDestroy(&final);

182:   PetscRandomDestroy(&rnd);
183:   PetscFree(indx3);
184:   PetscFree(indx4);
185:   PetscFinalize();
186:   return 0;
187: }

189: /*TEST

191:    build:
192:       requires: !mpiuni fftw !complex

194:    test:
195:       output_file: output/ex144.out

197:    test:
198:       suffix: 2
199:       nsize: 3
200:       output_file: output/ex144.out

202: TEST*/