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