Actual source code: ex142.c
petsc-3.13.6 2020-09-29
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 {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: PetscErrorCode ierr;
25: PetscScalar *x_array,*y_array,*z_array;
26: fftw_plan fplan,bplan;
28: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
29: #if defined(PETSC_USE_COMPLEX)
30: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
31: #endif
33: MPI_Comm_size(PETSC_COMM_WORLD, &size);
34: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
35: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex142");
36: PetscOptionsEList("-function", "Function type", "ex142", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
37: PetscOptionsBool("-vec_view draw", "View the functions", "ex142", view, &view, NULL);
38: function = (FuncType) func;
39: PetscOptionsEnd();
41: for (DIM = 0; DIM < ndim; DIM++) {
42: dim[DIM] = n; /* size of real space vector in DIM-dimension */
43: }
44: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
45: PetscRandomSetFromOptions(rdm);
47: for (DIM = 1; DIM < 5; DIM++) {
48: /* create vectors of length N=dim[0]*dim[1]* ...*dim[DIM-1] */
49: /*----------------------------------------------------------*/
50: N = Ny = 1;
51: for (i = 0; i < DIM-1; i++) {
52: N *= dim[i];
53: }
54: Ny = N; Ny *= 2*(dim[DIM-1]/2 + 1); /* add padding elements to output vector y */
55: N *= dim[DIM-1];
58: PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
59: VecCreateSeq(PETSC_COMM_SELF,N,&x);
60: PetscObjectSetName((PetscObject) x, "Real space vector");
62: VecCreateSeq(PETSC_COMM_SELF,Ny,&y);
63: PetscObjectSetName((PetscObject) y, "Frequency space vector");
65: VecDuplicate(x,&z);
66: PetscObjectSetName((PetscObject) z, "Reconstructed vector");
68: /* Set fftw plan */
69: /*----------------------------------*/
70: VecGetArray(x,&x_array);
71: VecGetArray(y,&y_array);
72: VecGetArray(z,&z_array);
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: PetscPrintf(PETSC_COMM_SELF,"DIM: %d, N %d, Ny %d\n",DIM,N,Ny);
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: }
98: VecRestoreArray(x,&x_array);
99: VecRestoreArray(y,&y_array);
100: VecRestoreArray(z,&z_array);
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: VecSetRandom(x, rdm);
108: } else if (function == CONSTANT) {
109: VecSet(x, 1.0);
110: } else if (function == TANH) {
111: VecGetArray(x, &x_array);
112: for (i = 0; i < N; ++i) {
113: x_array[i] = tanh((i - N/2.0)*(10.0/N));
114: }
115: VecRestoreArray(x, &x_array);
116: }
117: if (view) {
118: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
119: }
121: /* FFT - also test repeated transformation */
122: /*-------------------------------------------*/
123: VecGetArray(x,&x_array);
124: VecGetArray(y,&y_array);
125: VecGetArray(z,&z_array);
126: for (i=0; i<4; i++) {
127: /* FFTW_FORWARD */
128: fftw_execute(fplan);
130: /* FFTW_BACKWARD: destroys its input array 'y_array' even for out-of-place transforms! */
131: fftw_execute(bplan);
132: }
133: VecRestoreArray(x,&x_array);
134: VecRestoreArray(y,&y_array);
135: VecRestoreArray(z,&z_array);
137: /* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
138: /*------------------------------------------------------------------*/
139: s = 1.0/(PetscReal)N;
140: VecScale(z,s);
141: if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}
142: if (view) {VecView(z, PETSC_VIEWER_DRAW_WORLD);}
143: VecAXPY(z,-1.0,x);
144: VecNorm(z,NORM_1,&enorm);
145: if (enorm > 1.e-11) {
146: PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %g\n",(double)enorm);
147: }
149: /* free spaces */
150: fftw_destroy_plan(fplan);
151: fftw_destroy_plan(bplan);
152: VecDestroy(&x);
153: VecDestroy(&y);
154: VecDestroy(&z);
155: }
156: PetscRandomDestroy(&rdm);
157: PetscFinalize();
158: return ierr;
159: }
162: /*TEST
164: build:
165: requires: fftw !complex
167: test:
168: output_file: output/ex142.out
170: TEST*/