Actual source code: ex28.c
1: static char help[] = "Test sequential USFFT interface on a 3-dof field over a uniform DA 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 petscda.h
16: PetscInt main(PetscInt argc,char **args)
17: {
18: typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
19: const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
20: Mat A, AA;
21: PetscMPIInt size;
22: PetscInt N,i, stencil=1,dof=3;
23: PetscInt dim[3] = {10,10,10}, ndim = 3;
24: Vec coords,x,y,z,xx, yy, zz;
25: Vec xxsplit[DOF], yysplit[DOF], zzsplit[DOF];
26: PetscReal h[3];
27: PetscScalar s;
28: PetscRandom rdm;
29: PetscReal norm, enorm;
30: PetscInt func;
31: FuncType function = TANH;
32: DA da, da1, coordsda;
33: PetscTruth view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;
36: PetscInitialize(&argc,&args,(char *)0,help);
37: #if !defined(PETSC_USE_COMPLEX)
38: SETERRQ(PETSC_ERR_SUP, "This example requires complex numbers");
39: #endif
40: MPI_Comm_size(PETSC_COMM_WORLD, &size);
41: if (size != 1) SETERRQ(PETSC_ERR_SUP, "This is a uniprocessor example only!");
42: PetscOptionsBegin(PETSC_COMM_WORLD, PETSC_NULL, "USFFT Options", "ex27");
43: PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, PETSC_NULL);
44: function = (FuncType) func;
45: PetscOptionsEnd();
46: PetscOptionsGetTruth(PETSC_NULL,"-view_x",&view_x,PETSC_NULL);
47: PetscOptionsGetTruth(PETSC_NULL,"-view_y",&view_y,PETSC_NULL);
48: PetscOptionsGetTruth(PETSC_NULL,"-view_z",&view_z,PETSC_NULL);
49: PetscOptionsGetIntArray(PETSC_NULL,"-dim",dim,&ndim,PETSC_NULL);
51: // DA with the correct fiber dimension
52: DACreate3d(PETSC_COMM_SELF,DA_NONPERIODIC,DA_STENCIL_STAR,
53: dim[0], dim[1], dim[2],
54: PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
55: dof, stencil,
56: PETSC_NULL, PETSC_NULL, PETSC_NULL,
57: &da);
58:
59: // DA with fiber dimension 1 for split fields
60: DACreate3d(PETSC_COMM_SELF,DA_NONPERIODIC,DA_STENCIL_STAR,
61: dim[0], dim[1], dim[2],
62: PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE,
63: 1, stencil,
64: PETSC_NULL, PETSC_NULL, PETSC_NULL,
65: &da1);
66:
67:
68: // Coordinates
69: DAGetCoordinateDA(da, &coordsda);
70: DAGetGlobalVector(coordsda, &coords);
71: PetscObjectSetName((PetscObject) coords, "Grid coordinates");
72: for(i = 0, N = 1; i < 3; i++) {
73: h[i] = 1.0/dim[i];
74: PetscScalar *a;
75: VecGetArray(coords, &a);
76: PetscInt j,k,n = 0;
77: for(i = 0; i < 3; ++i) {
78: for(j = 0; j < dim[i]; ++j){
79: for(k = 0; k < 3; ++k) {
80: a[n] = j*h[i]; // coordinate along the j-th point in the i-th dimension
81: ++n;
82: }
83: }
84: }
85: VecRestoreArray(coords, &a);
87: }
88: DASetCoordinates(da, coords);
89: VecDestroy(coords);
91: // Work vectors
92: DAGetGlobalVector(da, &x);
93: PetscObjectSetName((PetscObject) x, "Real space vector");
94: DAGetGlobalVector(da, &xx);
95: PetscObjectSetName((PetscObject) xx, "Real space vector");
96: DAGetGlobalVector(da, &y);
97: PetscObjectSetName((PetscObject) y, "USFFT frequency space vector");
98: DAGetGlobalVector(da, &yy);
99: PetscObjectSetName((PetscObject) yy, "FFTW frequency space vector");
100: DAGetGlobalVector(da, &z);
101: PetscObjectSetName((PetscObject) z, "USFFT reconstructed vector");
102: DAGetGlobalVector(da, &zz);
103: PetscObjectSetName((PetscObject) zz, "FFTW reconstructed vector");
104: // Split vectors for FFTW
105: for(int ii = 0; ii < 3; ++ii) {
106: DAGetGlobalVector(da1, &xxsplit[ii]);
107: PetscObjectSetName((PetscObject) xxsplit[ii], "Real space split vector");
108: DAGetGlobalVector(da1, &yysplit[ii]);
109: PetscObjectSetName((PetscObject) yysplit[ii], "FFTW frequency space split vector");
110: DAGetGlobalVector(da1, &zzsplit[ii]);
111: PetscObjectSetName((PetscObject) zzsplit[ii], "FFTW reconstructed split vector");
112: }
115: PetscPrintf(PETSC_COMM_SELF, "%3-D: USFFT on vector of ");
116: for(i = 0, N = 1; i < 3; i++) {
117: PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ",i,dim[i]);
118: N *= dim[i];
119: }
120: PetscPrintf(PETSC_COMM_SELF, "; total size %d \n",N);
122:
123: if (function == RANDOM) {
124: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
125: PetscRandomSetFromOptions(rdm);
126: VecSetRandom(x, rdm);
127: PetscRandomDestroy(rdm);
128: }
129: else if (function == CONSTANT) {
130: VecSet(x, 1.0);
131: }
132: else if (function == TANH) {
133: PetscScalar *a;
134: VecGetArray(x, &a);
135: PetscInt j,k = 0;
136: for(i = 0; i < 3; ++i) {
137: for(j = 0; j < dim[i]; ++j) {
138: a[k] = tanh((j - dim[i]/2.0)*(10.0/dim[i]));
139: ++k;
140: }
141: }
142: VecRestoreArray(x, &a);
143: }
144: if(view_x) {
145: VecView(x, PETSC_VIEWER_STDOUT_WORLD);
146: }
147: VecCopy(x,xx);
148: // Split xx
149: VecStrideGatherAll(xx,xxsplit, INSERT_VALUES); //YES! 'Gather' means 'split' (or maybe 'scatter'?)!
151: VecNorm(x,NORM_2,&norm);
152: PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n",norm);
153:
154: /* create USFFT object */
155: MatCreateSeqUSFFT(da,da,&A);
156: /* create FFTW object */
157: MatCreateSeqFFTW(PETSC_COMM_SELF,3,dim,&AA);
158:
159: /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
160: MatMult(A,x,z);
161: for(int ii = 0; ii < 3; ++ii) {
162: MatMult(AA,xxsplit[ii],zzsplit[ii]);
163: }
164: // Now apply USFFT and FFTW forward several (3) times
165: for (i=0; i<3; ++i){
166: MatMult(A,x,y);
167: for(int ii = 0; ii < 3; ++ii) {
168: MatMult(AA,xxsplit[ii],yysplit[ii]);
169: }
170: MatMultTranspose(A,y,z);
171: for(int ii = 0; ii < 3; ++ii) {
172: MatMult(AA,yysplit[ii],zzsplit[ii]);
173: }
174: }
175: // Unsplit yy
176: VecStrideScatterAll(yysplit, yy, INSERT_VALUES); //YES! 'Scatter' means 'collect' (or maybe 'gather'?)!
177: // Unsplit zz
178: VecStrideScatterAll(zzsplit, zz, INSERT_VALUES); //YES! 'Scatter' means 'collect' (or maybe 'gather'?)!
180: if(view_y) {
181: PetscPrintf(PETSC_COMM_WORLD, "y = \n");
182: VecView(y, PETSC_VIEWER_STDOUT_WORLD);
183: PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
184: VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
185: }
186:
187: if(view_z) {
188: PetscPrintf(PETSC_COMM_WORLD, "z = \n");
189: VecView(z, PETSC_VIEWER_STDOUT_WORLD);
190: PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
191: VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
192: }
193:
194: /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
195: s = 1.0/(PetscReal)N;
196: VecScale(z,s);
197: VecAXPY(x,-1.0,z);
198: VecNorm(x,NORM_1,&enorm);
199: PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n",enorm);
201: /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
202: s = 1.0/(PetscReal)N;
203: VecScale(zz,s);
204: VecAXPY(xx,-1.0,zz);
205: VecNorm(xx,NORM_1,&enorm);
206: PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n",enorm);
208: /* compare y and yy: USFFT and FFTW results*/
209: VecNorm(y,NORM_2,&norm);
210: VecAXPY(y,-1.0,yy);
211: VecNorm(y,NORM_1,&enorm);
212: PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n",norm);
213: PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n",enorm);
214:
215: /* compare z and zz: USFFT and FFTW results*/
216: VecNorm(z,NORM_2,&norm);
217: VecAXPY(z,-1.0,zz);
218: VecNorm(z,NORM_1,&enorm);
219: PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n",norm);
220: PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n",enorm);
221:
223: /* free spaces */
224: DARestoreGlobalVector(da,&x);
225: DARestoreGlobalVector(da,&xx);
226: DARestoreGlobalVector(da,&y);
227: DARestoreGlobalVector(da,&yy);
228: DARestoreGlobalVector(da,&z);
229: DARestoreGlobalVector(da,&zz);
231: PetscFinalize();
232: return 0;
233: }