Actual source code: ex142.c

  1: static char help[] = "Test sequential r2c/c2r FFTW interface \n\n";

  3: /*
  4:   Compiling the code:
  5:       This code uses the real numbers version of PETSc
  6: */

  8: #include <petscmat.h>
  9: #include <fftw3.h>
 10: #include<fftw3-mpi.h>

 14: PetscInt main(PetscInt argc,char **args)
 15: {
 16:   typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
 17:   const char    *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
 18:   PetscMPIInt    size;
 19:   PetscInt       n = 10,N,Ny,ndim=4,dim[4],DIM,i;
 20:   Vec            x,y,z;
 21:   PetscScalar    s;
 22:   PetscRandom    rdm;
 23:   PetscReal      enorm;
 24:   PetscInt       func=RANDOM;
 25:   FuncType       function = RANDOM;
 26:   PetscBool      view = PETSC_FALSE;
 28:   PetscScalar    *x_array,*y_array,*z_array;
 29:   fftw_plan      fplan,bplan;
 30:   const ptrdiff_t N0 = 20, N1 = 20;

 32:   ptrdiff_t alloc_local, local_n0, local_0_start;
 33:   PetscInitialize(&argc,&args,(char *)0,help);
 34: #if defined(PETSC_USE_COMPLEX)
 35:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
 36: #endif
 37:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 38:   alloc_local=fftw_mpi_local_size_2d(N0, N1, PETSC_COMM_WORLD,
 39:                                               &local_n0, &local_0_start);

 41:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
 42:   PetscOptionsBegin(PETSC_COMM_WORLD, PETSC_NULL, "FFTW Options", "ex142");
 43:     PetscOptionsEList("-function", "Function type", "ex142", funcNames, NUM_FUNCS, funcNames[function], &func, PETSC_NULL);
 44:     PetscOptionsBool("-vec_view_draw", "View the functions", "ex112", view, &view, PETSC_NULL);
 45:     function = (FuncType) func;
 46:   PetscOptionsEnd();

 48:   for(DIM = 0; DIM < ndim; DIM++){
 49:     dim[DIM]  = n; /* size of real space vector in DIM-dimension */
 50:   }
 51:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 52:   PetscRandomSetFromOptions(rdm);

 54:   for(DIM = 1; DIM < 5; DIM++){
 55:     /* create vectors of length N=dim[0]*dim[1]* ...*dim[DIM-1] */
 56:     /*----------------------------------------------------------*/
 57:     N = Ny = 1;
 58:     for(i = 0; i < DIM-1; i++) {
 59:       N *= dim[i];
 60:     }
 61:     Ny = N; Ny *= 2*(dim[DIM-1]/2 + 1); /* add padding elements to output vector y */
 62:     N *= dim[DIM-1];
 63: 

 65:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
 66:     VecCreateSeq(PETSC_COMM_SELF,N,&x);
 67:     PetscObjectSetName((PetscObject) x, "Real space vector");
 68: 
 69:     VecCreateSeq(PETSC_COMM_SELF,Ny,&y);
 70:     PetscObjectSetName((PetscObject) y, "Frequency space vector");

 72:     VecDuplicate(x,&z);
 73:     PetscObjectSetName((PetscObject) z, "Reconstructed vector");
 74: 

 76:     /* Set fftw plan                    */
 77:     /*----------------------------------*/
 78:     VecGetArray(x,&x_array);
 79:     VecGetArray(y,&y_array);
 80:     VecGetArray(z,&z_array);

 82:     unsigned int flags = FFTW_ESTIMATE; //or FFTW_MEASURE
 83:     /* The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so such planning
 84:      should be done before the input is initialized by the user. */
 85:     printf("DIM: %d, N %d, Ny %d\n",DIM,N,Ny);

 87:     switch (DIM){
 88:     case 1:
 89:       fplan = fftw_plan_dft_r2c_1d(dim[0], (double *)x_array, (fftw_complex*)y_array, flags);
 90:       bplan = fftw_plan_dft_c2r_1d(dim[0], (fftw_complex*)y_array, (double *)z_array, flags);
 91:       break;
 92:     case 2:
 93:       fplan = fftw_plan_dft_r2c_2d(dim[0],dim[1],(double *)x_array, (fftw_complex*)y_array,flags);
 94:       bplan = fftw_plan_dft_c2r_2d(dim[0],dim[1],(fftw_complex*)y_array,(double *)z_array,flags);
 95:       break;
 96:     case 3:
 97:       fplan = fftw_plan_dft_r2c_3d(dim[0],dim[1],dim[2],(double *)x_array, (fftw_complex*)y_array,flags);
 98:       bplan = fftw_plan_dft_c2r_3d(dim[0],dim[1],dim[2],(fftw_complex*)y_array,(double *)z_array,flags);
 99:       break;
100:     default:
101:       fplan = fftw_plan_dft_r2c(DIM,dim,(double *)x_array, (fftw_complex*)y_array,flags);
102:       bplan = fftw_plan_dft_c2r(DIM,dim,(fftw_complex*)y_array,(double *)z_array,flags);
103:       break;
104:     }

106:     VecRestoreArray(x,&x_array);
107:     VecRestoreArray(y,&y_array);
108:     VecRestoreArray(z,&z_array);

110:     /* Initialize Real space vector x: 
111:        The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so planning
112:        should be done before the input is initialized by the user.
113:     --------------------------------------------------------*/
114:     if (function == RANDOM) {
115:       VecSetRandom(x, rdm);
116:     } else if (function == CONSTANT) {
117:       VecSet(x, 1.0);
118:     } else if (function == TANH) {
119:       VecGetArray(x, &x_array);
120:       for(i = 0; i < N; ++i) {
121:         x_array[i] = tanh((i - N/2.0)*(10.0/N));
122:       }
123:       VecRestoreArray(x, &x_array);
124:     }
125:     if (view) {
126:       VecView(x, PETSC_VIEWER_STDOUT_WORLD);
127:     }

129:     /* FFT - also test repeated transformation   */
130:     /*-------------------------------------------*/
131:     VecGetArray(x,&x_array);
132:     VecGetArray(y,&y_array);
133:     VecGetArray(z,&z_array);
134:     for (i=0; i<3; i++){
135:       /* FFTW_FORWARD */
136:       fftw_execute(fplan);
137:       //printf("\n fout:\n");
138:       //fftw_complex* fout = (fftw_complex*)y_array;
139:       //for (i=0; i<N/2+1; i++) printf("%d (%g %g)\n",i,fout[i][0],fout[i][1]);
140: 
141:       /* FFTW_BACKWARD: destroys its input array 'y_array' even for out-of-place transforms! */
142:       fftw_execute(bplan);
143:     }
144:     VecRestoreArray(x,&x_array);
145:     VecRestoreArray(y,&y_array);
146:     VecRestoreArray(z,&z_array);

148:     /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
149:     /*------------------------------------------------------------------*/
150:     s = 1.0/(PetscReal)N;
151:     VecScale(z,s);
152:     if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}
153:     if (view) {VecView(z, PETSC_VIEWER_DRAW_WORLD);}
154:     VecAXPY(z,-1.0,x);
155:     VecNorm(z,NORM_1,&enorm);
156:     if (enorm > 1.e-11){
157:       PetscPrintf(PETSC_COMM_SELF,"  Error norm of |x - z| %A\n",enorm);
158:     }

160:     /* free spaces */
161:     fftw_destroy_plan(fplan);
162:     fftw_destroy_plan(bplan);
163:     VecDestroy(&x);
164:     VecDestroy(&y);
165:     VecDestroy(&z);
166:   }
167:   PetscRandomDestroy(&rdm);
168:   PetscFinalize();
169:   return 0;
170: }