Actual source code: ex28.c
1: static char help[] = "Test sequential USFFT interface on a 3-dof field over a uniform DMDA and compares to the result of FFTW acting on a split version of the field\n\n";
3: /*
4: Compiling the code:
5: This code uses the complex numbers version of PETSc and the FFTW package, so configure
6: must be run to enable these.
8: */
10: #define DOF 3
12: #include <petscmat.h>
13: #include <petscdm.h>
14: #include <petscdmda.h>
15: int main(int argc,char **args)
16: {
17: typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
18: const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
19: Mat A, AA;
20: PetscMPIInt size;
21: PetscInt N,i, stencil=1,dof=3;
22: PetscInt dim[3] = {10,10,10}, ndim = 3;
23: Vec coords,x,y,z,xx, yy, zz;
24: Vec xxsplit[DOF], yysplit[DOF], zzsplit[DOF];
25: PetscReal h[3];
26: PetscScalar s;
27: PetscRandom rdm;
28: PetscReal norm, enorm;
29: PetscInt func,ii;
30: FuncType function = TANH;
31: DM da, da1, coordsda;
32: PetscBool view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;
35: PetscInitialize(&argc,&args,(char*)0,help);
36: MPI_Comm_size(PETSC_COMM_WORLD, &size);
38: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
39: PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
40: function = (FuncType) func;
41: PetscOptionsEnd();
42: PetscOptionsGetBool(NULL,NULL,"-view_x",&view_x,NULL);
43: PetscOptionsGetBool(NULL,NULL,"-view_y",&view_y,NULL);
44: PetscOptionsGetBool(NULL,NULL,"-view_z",&view_z,NULL);
45: PetscOptionsGetIntArray(NULL,NULL,"-dim",dim,&ndim,NULL);
47: /* DMDA with the correct fiber dimension */
48: DMDACreate3d(PETSC_COMM_SELF,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,dim[0],dim[1],dim[2],PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
49: dof,stencil,NULL,NULL,NULL,&da);
50: DMSetFromOptions(da);
51: DMSetUp(da);
52: /* DMDA with fiber dimension 1 for split fields */
53: DMDACreate3d(PETSC_COMM_SELF,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,dim[0],dim[1],dim[2],PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
54: 1,stencil,NULL,NULL,NULL,&da1);
55: DMSetFromOptions(da1);
56: DMSetUp(da1);
58: /* Coordinates */
59: DMGetCoordinateDM(da,&coordsda);
60: DMGetGlobalVector(coordsda,&coords);
61: PetscObjectSetName((PetscObject) coords,"Grid coordinates");
62: for (i = 0, N = 1; i < 3; i++) {
63: h[i] = 1.0/dim[i];
64: PetscScalar *a;
65: VecGetArray(coords, &a);
66: PetscInt j,k,n = 0;
67: for (i = 0; i < 3; ++i) {
68: for (j = 0; j < dim[i]; ++j) {
69: for (k = 0; k < 3; ++k) {
70: a[n] = j*h[i]; /* coordinate along the j-th point in the i-th dimension */
71: ++n;
72: }
73: }
74: }
75: VecRestoreArray(coords, &a);
77: }
78: DMSetCoordinates(da, coords);
79: VecDestroy(&coords);
81: /* Work vectors */
82: DMGetGlobalVector(da, &x);
83: PetscObjectSetName((PetscObject) x, "Real space vector");
84: DMGetGlobalVector(da, &xx);
85: PetscObjectSetName((PetscObject) xx, "Real space vector");
86: DMGetGlobalVector(da, &y);
87: PetscObjectSetName((PetscObject) y, "USFFT frequency space vector");
88: DMGetGlobalVector(da, &yy);
89: PetscObjectSetName((PetscObject) yy, "FFTW frequency space vector");
90: DMGetGlobalVector(da, &z);
91: PetscObjectSetName((PetscObject) z, "USFFT reconstructed vector");
92: DMGetGlobalVector(da, &zz);
93: PetscObjectSetName((PetscObject) zz, "FFTW reconstructed vector");
94: /* Split vectors for FFTW */
95: for (ii = 0; ii < 3; ++ii) {
96: DMGetGlobalVector(da1, &xxsplit[ii]);
97: PetscObjectSetName((PetscObject) xxsplit[ii], "Real space split vector");
98: DMGetGlobalVector(da1, &yysplit[ii]);
99: PetscObjectSetName((PetscObject) yysplit[ii], "FFTW frequency space split vector");
100: DMGetGlobalVector(da1, &zzsplit[ii]);
101: PetscObjectSetName((PetscObject) zzsplit[ii], "FFTW reconstructed split vector");
102: }
104: PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");
105: for (i = 0, N = 1; i < 3; i++) {
106: PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);
107: N *= dim[i];
108: }
109: PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);
111: if (function == RANDOM) {
112: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
113: PetscRandomSetFromOptions(rdm);
114: VecSetRandom(x, rdm);
115: PetscRandomDestroy(&rdm);
116: } else if (function == CONSTANT) {
117: VecSet(x, 1.0);
118: } else if (function == TANH) {
119: PetscScalar *a;
120: VecGetArray(x, &a);
121: PetscInt j,k = 0;
122: for (i = 0; i < 3; ++i) {
123: for (j = 0; j < dim[i]; ++j) {
124: a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
125: ++k;
126: }
127: }
128: VecRestoreArray(x, &a);
129: }
130: if (view_x) {
131: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
132: }
133: VecCopy(x,xx);
134: /* Split xx */
135: VecStrideGatherAll(xx,xxsplit, INSERT_VALUES); /*YES! 'Gather' means 'split' (or maybe 'scatter'?)! */
137: VecNorm(x,NORM_2,&norm);
138: PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);
140: /* create USFFT object */
141: MatCreateSeqUSFFT(da,da,&A);
142: /* create FFTW object */
143: MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);
145: /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
146: MatMult(A,x,z);
147: for (ii = 0; ii < 3; ++ii) {
148: MatMult(AA,xxsplit[ii],zzsplit[ii]);
149: }
150: /* Now apply USFFT and FFTW forward several (3) times */
151: for (i=0; i<3; ++i) {
152: MatMult(A,x,y);
153: for (ii = 0; ii < 3; ++ii) {
154: MatMult(AA,xxsplit[ii],yysplit[ii]);
155: }
156: MatMultTranspose(A,y,z);
157: for (ii = 0; ii < 3; ++ii) {
158: MatMult(AA,yysplit[ii],zzsplit[ii]);
159: }
160: }
161: /* Unsplit yy */
162: VecStrideScatterAll(yysplit, yy, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
163: /* Unsplit zz */
164: VecStrideScatterAll(zzsplit, zz, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
166: if (view_y) {
167: PetscPrintf(PETSC_COMM_WORLD, "y = \n");
168: VecView(y, PETSC_VIEWER_STDOUT_WORLD);
169: PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
170: VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
171: }
173: if (view_z) {
174: PetscPrintf(PETSC_COMM_WORLD, "z = \n");
175: VecView(z, PETSC_VIEWER_STDOUT_WORLD);
176: PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
177: VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
178: }
180: /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
181: s = 1.0/(PetscReal)N;
182: VecScale(z,s);
183: VecAXPY(x,-1.0,z);
184: VecNorm(x,NORM_1,&enorm);
185: PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);
187: /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
188: s = 1.0/(PetscReal)N;
189: VecScale(zz,s);
190: VecAXPY(xx,-1.0,zz);
191: VecNorm(xx,NORM_1,&enorm);
192: PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);
194: /* compare y and yy: USFFT and FFTW results*/
195: VecNorm(y,NORM_2,&norm);
196: VecAXPY(y,-1.0,yy);
197: VecNorm(y,NORM_1,&enorm);
198: PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
199: PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);
201: /* compare z and zz: USFFT and FFTW results*/
202: VecNorm(z,NORM_2,&norm);
203: VecAXPY(z,-1.0,zz);
204: VecNorm(z,NORM_1,&enorm);
205: PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
206: PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);
208: /* free spaces */
209: DMRestoreGlobalVector(da,&x);
210: DMRestoreGlobalVector(da,&xx);
211: DMRestoreGlobalVector(da,&y);
212: DMRestoreGlobalVector(da,&yy);
213: DMRestoreGlobalVector(da,&z);
214: DMRestoreGlobalVector(da,&zz);
216: PetscFinalize();
217: return 0;
218: }