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 { 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; 31 32 PetscFunctionBeginUser; 33 PetscCall(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 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(); 45 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)); 49 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]; 58 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")); 62 63 PetscCall(VecCreateSeq(PETSC_COMM_SELF, Ny, &y)); 64 PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector")); 65 66 PetscCall(VecDuplicate(x, &z)); 67 PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector")); 68 69 /* Set fftw plan */ 70 /*----------------------------------*/ 71 PetscCall(VecGetArray(x, &x_array)); 72 PetscCall(VecGetArray(y, &y_array)); 73 PetscCall(VecGetArray(z, &z_array)); 74 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)); 79 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 } 98 99 PetscCall(VecRestoreArray(x, &x_array)); 100 PetscCall(VecRestoreArray(y, &y_array)); 101 PetscCall(VecRestoreArray(z, &z_array)); 102 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)); 117 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); 126 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)); 133 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)); 143 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 } 155 156 /*TEST 157 158 build: 159 requires: fftw !complex 160 161 test: 162 output_file: output/ex142.out 163 164 TEST*/ 165