Actual source code: ex142.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: static char help[] = "Test sequential r2c/c2r FFTW without PETSc 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>

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

 30:   PetscInitialize(&argc,&args,(char*)0,help);
 31: #if defined(PETSC_USE_COMPLEX)
 32:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
 33: #endif

 35:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 36:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
 37:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex142");
 38:   PetscOptionsEList("-function", "Function type", "ex142", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
 39:   PetscOptionsBool("-vec_view draw", "View the functions", "ex142", view, &view, NULL);
 40:   function = (FuncType) func;
 41:   PetscOptionsEnd();

 43:   for (DIM = 0; DIM < ndim; DIM++) {
 44:     dim[DIM] = n;  /* size of real space vector in DIM-dimension */
 45:   }
 46:   PetscRandomCreate(PETSC_COMM_SELF, &rdm);
 47:   PetscRandomSetFromOptions(rdm);

 49:   for (DIM = 1; DIM < 5; DIM++) {
 50:     /* create vectors of length N=dim[0]*dim[1]* ...*dim[DIM-1] */
 51:     /*----------------------------------------------------------*/
 52:     N = Ny = 1;
 53:     for (i = 0; i < DIM-1; i++) {
 54:       N *= dim[i];
 55:     }
 56:     Ny = N; Ny *= 2*(dim[DIM-1]/2 + 1); /* add padding elements to output vector y */
 57:     N *= dim[DIM-1];


 60:     PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
 61:     VecCreateSeq(PETSC_COMM_SELF,N,&x);
 62:     PetscObjectSetName((PetscObject) x, "Real space vector");

 64:     VecCreateSeq(PETSC_COMM_SELF,Ny,&y);
 65:     PetscObjectSetName((PetscObject) y, "Frequency space vector");

 67:     VecDuplicate(x,&z);
 68:     PetscObjectSetName((PetscObject) z, "Reconstructed vector");

 70:     /* Set fftw plan                    */
 71:     /*----------------------------------*/
 72:     VecGetArray(x,&x_array);
 73:     VecGetArray(y,&y_array);
 74:     VecGetArray(z,&z_array);

 76:     unsigned int flags = FFTW_ESTIMATE; /*or FFTW_MEASURE */
 77:     /* The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so such planning
 78:      should be done before the input is initialized by the user. */
 79:     PetscPrintf(PETSC_COMM_SELF,"DIM: %d, N %d, Ny %d\n",DIM,N,Ny);

 81:     switch (DIM) {
 82:     case 1:
 83:       fplan = fftw_plan_dft_r2c_1d(dim[0], (double*)x_array, (fftw_complex*)y_array, flags);
 84:       bplan = fftw_plan_dft_c2r_1d(dim[0], (fftw_complex*)y_array, (double*)z_array, flags);
 85:       break;
 86:     case 2:
 87:       fplan = fftw_plan_dft_r2c_2d(dim[0],dim[1],(double*)x_array, (fftw_complex*)y_array,flags);
 88:       bplan = fftw_plan_dft_c2r_2d(dim[0],dim[1],(fftw_complex*)y_array,(double*)z_array,flags);
 89:       break;
 90:     case 3:
 91:       fplan = fftw_plan_dft_r2c_3d(dim[0],dim[1],dim[2],(double*)x_array, (fftw_complex*)y_array,flags);
 92:       bplan = fftw_plan_dft_c2r_3d(dim[0],dim[1],dim[2],(fftw_complex*)y_array,(double*)z_array,flags);
 93:       break;
 94:     default:
 95:       fplan = fftw_plan_dft_r2c(DIM,(int*)dim,(double*)x_array, (fftw_complex*)y_array,flags);
 96:       bplan = fftw_plan_dft_c2r(DIM,(int*)dim,(fftw_complex*)y_array,(double*)z_array,flags);
 97:       break;
 98:     }

100:     VecRestoreArray(x,&x_array);
101:     VecRestoreArray(y,&y_array);
102:     VecRestoreArray(z,&z_array);

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

123:     /* FFT - also test repeated transformation   */
124:     /*-------------------------------------------*/
125:     VecGetArray(x,&x_array);
126:     VecGetArray(y,&y_array);
127:     VecGetArray(z,&z_array);
128:     for (i=0; i<4; i++) {
129:       /* FFTW_FORWARD */
130:       fftw_execute(fplan);
131: 
132:       /* FFTW_BACKWARD: destroys its input array 'y_array' even for out-of-place transforms! */
133:       fftw_execute(bplan);
134:     }
135:     VecRestoreArray(x,&x_array);
136:     VecRestoreArray(y,&y_array);
137:     VecRestoreArray(z,&z_array);

139:     /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
140:     /*------------------------------------------------------------------*/
141:     s    = 1.0/(PetscReal)N;
142:     VecScale(z,s);
143:     if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}
144:     if (view) {VecView(z, PETSC_VIEWER_DRAW_WORLD);}
145:     VecAXPY(z,-1.0,x);
146:     VecNorm(z,NORM_1,&enorm);
147:     if (enorm > 1.e-11) {
148:       PetscPrintf(PETSC_COMM_SELF,"  Error norm of |x - z| %g\n",(double)enorm);
149:     }

151:     /* free spaces */
152:     fftw_destroy_plan(fplan);
153:     fftw_destroy_plan(bplan);
154:     VecDestroy(&x);
155:     VecDestroy(&y);
156:     VecDestroy(&z);
157:   }
158:   PetscRandomDestroy(&rdm);
159:   PetscFinalize();
160:   return 0;
161: }