Actual source code: ex142.c

  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>

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

 32:   PetscFunctionBeginUser;
 33:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 34: #if defined(PETSC_USE_COMPLEX)
 35:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers");
 36: #endif

 38:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 39:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
 40:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex142");
 41:   PetscCall(PetscOptionsEList("-function", "Function type", "ex142", funcNames, NUM_FUNCS, funcNames[function], &func, NULL));
 42:   PetscCall(PetscOptionsBool("-vec_view draw", "View the functions", "ex142", view, &view, NULL));
 43:   function = (FuncType)func;
 44:   PetscOptionsEnd();

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

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

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

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

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

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

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

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

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

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

118:     /* FFT - also test repeated transformation   */
119:     /*-------------------------------------------*/
120:     PetscCall(VecGetArray(x, &x_array));
121:     PetscCall(VecGetArray(y, &y_array));
122:     PetscCall(VecGetArray(z, &z_array));
123:     for (i = 0; i < 4; i++) {
124:       /* FFTW_FORWARD */
125:       fftw_execute(fplan);

127:       /* FFTW_BACKWARD: destroys its input array 'y_array' even for out-of-place transforms! */
128:       fftw_execute(bplan);
129:     }
130:     PetscCall(VecRestoreArray(x, &x_array));
131:     PetscCall(VecRestoreArray(y, &y_array));
132:     PetscCall(VecRestoreArray(z, &z_array));

134:     /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
135:     /*------------------------------------------------------------------*/
136:     s = 1.0 / (PetscReal)N;
137:     PetscCall(VecScale(z, s));
138:     if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
139:     if (view) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD));
140:     PetscCall(VecAXPY(z, -1.0, x));
141:     PetscCall(VecNorm(z, NORM_1, &enorm));
142:     if (enorm > 1.e-11) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Error norm of |x - z| %g\n", (double)enorm));

144:     /* free spaces */
145:     fftw_destroy_plan(fplan);
146:     fftw_destroy_plan(bplan);
147:     PetscCall(VecDestroy(&x));
148:     PetscCall(VecDestroy(&y));
149:     PetscCall(VecDestroy(&z));
150:   }
151:   PetscCall(PetscRandomDestroy(&rdm));
152:   PetscCall(PetscFinalize());
153:   return 0;
154: }

156: /*TEST

158:    build:
159:      requires: fftw !complex

161:    test:
162:      output_file: output/ex142.out

164: TEST*/