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