1 static char help[] = "Test sequential r2c/c2r FFTW without PETSc interface \n\n"; 2 3 /* 4 Compiling the code: 5 This code uses the real numbers version of PETSc 6 */ 7 8 #include <petscmat.h> 9 #include <fftw3.h> 10 11 int main(int argc, char **args) { 12 typedef enum { 13 RANDOM, 14 CONSTANT, 15 TANH, 16 NUM_FUNCS 17 } FuncType; 18 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"}; 19 PetscMPIInt size; 20 int n = 10, N, Ny, ndim = 4, i, dim[4], DIM; 21 Vec x, y, z; 22 PetscScalar s; 23 PetscRandom rdm; 24 PetscReal enorm; 25 PetscInt func = RANDOM; 26 FuncType function = RANDOM; 27 PetscBool view = PETSC_FALSE; 28 PetscScalar *x_array, *y_array, *z_array; 29 fftw_plan fplan, bplan; 30 31 PetscFunctionBeginUser; 32 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 33 #if defined(PETSC_USE_COMPLEX) 34 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires real numbers"); 35 #endif 36 37 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 38 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 39 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex142"); 40 PetscCall(PetscOptionsEList("-function", "Function type", "ex142", funcNames, NUM_FUNCS, funcNames[function], &func, NULL)); 41 PetscCall(PetscOptionsBool("-vec_view draw", "View the functions", "ex142", view, &view, NULL)); 42 function = (FuncType)func; 43 PetscOptionsEnd(); 44 45 for (DIM = 0; DIM < ndim; DIM++) { dim[DIM] = n; /* size of real space vector in DIM-dimension */ } 46 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm)); 47 PetscCall(PetscRandomSetFromOptions(rdm)); 48 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++) { N *= dim[i]; } 54 Ny = N; 55 Ny *= 2 * (dim[DIM - 1] / 2 + 1); /* add padding elements to output vector y */ 56 N *= dim[DIM - 1]; 57 58 PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N)); 59 PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x)); 60 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector")); 61 62 PetscCall(VecCreateSeq(PETSC_COMM_SELF, Ny, &y)); 63 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector")); 64 65 PetscCall(VecDuplicate(x, &z)); 66 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector")); 67 68 /* Set fftw plan */ 69 /*----------------------------------*/ 70 PetscCall(VecGetArray(x, &x_array)); 71 PetscCall(VecGetArray(y, &y_array)); 72 PetscCall(VecGetArray(z, &z_array)); 73 74 unsigned int flags = FFTW_ESTIMATE; /*or FFTW_MEASURE */ 75 /* The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so such planning 76 should be done before the input is initialized by the user. */ 77 PetscCall(PetscPrintf(PETSC_COMM_SELF, "DIM: %d, N %d, Ny %d\n", DIM, N, Ny)); 78 79 switch (DIM) { 80 case 1: 81 fplan = fftw_plan_dft_r2c_1d(dim[0], (double *)x_array, (fftw_complex *)y_array, flags); 82 bplan = fftw_plan_dft_c2r_1d(dim[0], (fftw_complex *)y_array, (double *)z_array, flags); 83 break; 84 case 2: 85 fplan = fftw_plan_dft_r2c_2d(dim[0], dim[1], (double *)x_array, (fftw_complex *)y_array, flags); 86 bplan = fftw_plan_dft_c2r_2d(dim[0], dim[1], (fftw_complex *)y_array, (double *)z_array, flags); 87 break; 88 case 3: 89 fplan = fftw_plan_dft_r2c_3d(dim[0], dim[1], dim[2], (double *)x_array, (fftw_complex *)y_array, flags); 90 bplan = fftw_plan_dft_c2r_3d(dim[0], dim[1], dim[2], (fftw_complex *)y_array, (double *)z_array, flags); 91 break; 92 default: 93 fplan = fftw_plan_dft_r2c(DIM, (int *)dim, (double *)x_array, (fftw_complex *)y_array, flags); 94 bplan = fftw_plan_dft_c2r(DIM, (int *)dim, (fftw_complex *)y_array, (double *)z_array, flags); 95 break; 96 } 97 98 PetscCall(VecRestoreArray(x, &x_array)); 99 PetscCall(VecRestoreArray(y, &y_array)); 100 PetscCall(VecRestoreArray(z, &z_array)); 101 102 /* Initialize Real space vector x: 103 The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so planning 104 should be done before the input is initialized by the user. 105 --------------------------------------------------------*/ 106 if (function == RANDOM) { 107 PetscCall(VecSetRandom(x, rdm)); 108 } else if (function == CONSTANT) { 109 PetscCall(VecSet(x, 1.0)); 110 } else if (function == TANH) { 111 PetscCall(VecGetArray(x, &x_array)); 112 for (i = 0; i < N; ++i) { x_array[i] = tanh((i - N / 2.0) * (10.0 / N)); } 113 PetscCall(VecRestoreArray(x, &x_array)); 114 } 115 if (view) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 116 117 /* FFT - also test repeated transformation */ 118 /*-------------------------------------------*/ 119 PetscCall(VecGetArray(x, &x_array)); 120 PetscCall(VecGetArray(y, &y_array)); 121 PetscCall(VecGetArray(z, &z_array)); 122 for (i = 0; i < 4; i++) { 123 /* FFTW_FORWARD */ 124 fftw_execute(fplan); 125 126 /* FFTW_BACKWARD: destroys its input array 'y_array' even for out-of-place transforms! */ 127 fftw_execute(bplan); 128 } 129 PetscCall(VecRestoreArray(x, &x_array)); 130 PetscCall(VecRestoreArray(y, &y_array)); 131 PetscCall(VecRestoreArray(z, &z_array)); 132 133 /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */ 134 /*------------------------------------------------------------------*/ 135 s = 1.0 / (PetscReal)N; 136 PetscCall(VecScale(z, s)); 137 if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD)); 138 if (view) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD)); 139 PetscCall(VecAXPY(z, -1.0, x)); 140 PetscCall(VecNorm(z, NORM_1, &enorm)); 141 if (enorm > 1.e-11) { PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm)); } 142 143 /* free spaces */ 144 fftw_destroy_plan(fplan); 145 fftw_destroy_plan(bplan); 146 PetscCall(VecDestroy(&x)); 147 PetscCall(VecDestroy(&y)); 148 PetscCall(VecDestroy(&z)); 149 } 150 PetscCall(PetscRandomDestroy(&rdm)); 151 PetscCall(PetscFinalize()); 152 return 0; 153 } 154 155 /*TEST 156 157 build: 158 requires: fftw !complex 159 160 test: 161 output_file: output/ex142.out 162 163 TEST*/ 164