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 {
18: RANDOM,
19: CONSTANT,
20: TANH,
21: NUM_FUNCS
22: } FuncType;
23: const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
24: Mat A, AA;
25: PetscMPIInt size;
26: PetscInt N, i, stencil = 1, dof = 3;
27: PetscInt dim[3] = {10, 10, 10}, ndim = 3;
28: Vec coords, x, y, z, xx, yy, zz;
29: Vec xxsplit[DOF], yysplit[DOF], zzsplit[DOF];
30: PetscReal h[3];
31: PetscScalar s;
32: PetscRandom rdm;
33: PetscReal norm, enorm;
34: PetscInt func, ii;
35: FuncType function = TANH;
36: DM da, da1, coordsda;
37: PetscBool view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;
40: PetscInitialize(&argc, &args, (char *)0, help);
41: MPI_Comm_size(PETSC_COMM_WORLD, &size);
43: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
44: PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
45: function = (FuncType)func;
46: PetscOptionsEnd();
47: PetscOptionsGetBool(NULL, NULL, "-view_x", &view_x, NULL);
48: PetscOptionsGetBool(NULL, NULL, "-view_y", &view_y, NULL);
49: PetscOptionsGetBool(NULL, NULL, "-view_z", &view_z, NULL);
50: PetscOptionsGetIntArray(NULL, NULL, "-dim", dim, &ndim, NULL);
52: /* DMDA with the correct fiber dimension */
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, dof, stencil, NULL, NULL, NULL, &da);
54: DMSetFromOptions(da);
55: DMSetUp(da);
56: /* DMDA with fiber dimension 1 for split fields */
57: 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, 1, stencil, NULL, NULL, NULL, &da1);
58: DMSetFromOptions(da1);
59: DMSetUp(da1);
61: /* Coordinates */
62: DMGetCoordinateDM(da, &coordsda);
63: DMGetGlobalVector(coordsda, &coords);
64: PetscObjectSetName((PetscObject)coords, "Grid coordinates");
65: for (i = 0, N = 1; i < 3; i++) {
66: h[i] = 1.0 / dim[i];
67: PetscScalar *a;
68: VecGetArray(coords, &a);
69: PetscInt j, k, n = 0;
70: for (i = 0; i < 3; ++i) {
71: for (j = 0; j < dim[i]; ++j) {
72: for (k = 0; k < 3; ++k) {
73: a[n] = j * h[i]; /* coordinate along the j-th point in the i-th dimension */
74: ++n;
75: }
76: }
77: }
78: VecRestoreArray(coords, &a);
79: }
80: DMSetCoordinates(da, coords);
81: VecDestroy(&coords);
83: /* Work vectors */
84: DMGetGlobalVector(da, &x);
85: PetscObjectSetName((PetscObject)x, "Real space vector");
86: DMGetGlobalVector(da, &xx);
87: PetscObjectSetName((PetscObject)xx, "Real space vector");
88: DMGetGlobalVector(da, &y);
89: PetscObjectSetName((PetscObject)y, "USFFT frequency space vector");
90: DMGetGlobalVector(da, &yy);
91: PetscObjectSetName((PetscObject)yy, "FFTW frequency space vector");
92: DMGetGlobalVector(da, &z);
93: PetscObjectSetName((PetscObject)z, "USFFT reconstructed vector");
94: DMGetGlobalVector(da, &zz);
95: PetscObjectSetName((PetscObject)zz, "FFTW reconstructed vector");
96: /* Split vectors for FFTW */
97: for (ii = 0; ii < 3; ++ii) {
98: DMGetGlobalVector(da1, &xxsplit[ii]);
99: PetscObjectSetName((PetscObject)xxsplit[ii], "Real space split vector");
100: DMGetGlobalVector(da1, &yysplit[ii]);
101: PetscObjectSetName((PetscObject)yysplit[ii], "FFTW frequency space split vector");
102: DMGetGlobalVector(da1, &zzsplit[ii]);
103: PetscObjectSetName((PetscObject)zzsplit[ii], "FFTW reconstructed split vector");
104: }
106: PetscPrintf(PETSC_COMM_SELF, "%3-" PetscInt_FMT ": USFFT on vector of ");
107: for (i = 0, N = 1; i < 3; i++) {
108: PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ", i, dim[i]);
109: N *= dim[i];
110: }
111: PetscPrintf(PETSC_COMM_SELF, "; total size %d \n", N);
113: if (function == RANDOM) {
114: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
115: PetscRandomSetFromOptions(rdm);
116: VecSetRandom(x, rdm);
117: PetscRandomDestroy(&rdm);
118: } else if (function == CONSTANT) {
119: VecSet(x, 1.0);
120: } else if (function == TANH) {
121: PetscScalar *a;
122: VecGetArray(x, &a);
123: PetscInt j, k = 0;
124: for (i = 0; i < 3; ++i) {
125: for (j = 0; j < dim[i]; ++j) {
126: a[k] = tanh((j - dim[i] / 2.0) * (10.0 / dim[i]));
127: ++k;
128: }
129: }
130: VecRestoreArray(x, &a);
131: }
132: if (view_x) VecView(x, PETSC_VIEWER_STDOUT_WORLD);
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) MatMult(AA, xxsplit[ii], zzsplit[ii]);
148: /* Now apply USFFT and FFTW forward several (3) times */
149: for (i = 0; i < 3; ++i) {
150: MatMult(A, x, y);
151: for (ii = 0; ii < 3; ++ii) MatMult(AA, xxsplit[ii], yysplit[ii]);
152: MatMultTranspose(A, y, z);
153: for (ii = 0; ii < 3; ++ii) MatMult(AA, yysplit[ii], zzsplit[ii]);
154: }
155: /* Unsplit yy */
156: VecStrideScatterAll(yysplit, yy, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
157: /* Unsplit zz */
158: VecStrideScatterAll(zzsplit, zz, INSERT_VALUES); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
160: if (view_y) {
161: PetscPrintf(PETSC_COMM_WORLD, "y = \n");
162: VecView(y, PETSC_VIEWER_STDOUT_WORLD);
163: PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
164: VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
165: }
167: if (view_z) {
168: PetscPrintf(PETSC_COMM_WORLD, "z = \n");
169: VecView(z, PETSC_VIEWER_STDOUT_WORLD);
170: PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
171: VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
172: }
174: /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
175: s = 1.0 / (PetscReal)N;
176: VecScale(z, s);
177: VecAXPY(x, -1.0, z);
178: VecNorm(x, NORM_1, &enorm);
179: PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n", enorm);
181: /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
182: s = 1.0 / (PetscReal)N;
183: VecScale(zz, s);
184: VecAXPY(xx, -1.0, zz);
185: VecNorm(xx, NORM_1, &enorm);
186: PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n", enorm);
188: /* compare y and yy: USFFT and FFTW results*/
189: VecNorm(y, NORM_2, &norm);
190: VecAXPY(y, -1.0, yy);
191: VecNorm(y, NORM_1, &enorm);
192: PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n", norm);
193: PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n", enorm);
195: /* compare z and zz: USFFT and FFTW results*/
196: VecNorm(z, NORM_2, &norm);
197: VecAXPY(z, -1.0, zz);
198: VecNorm(z, NORM_1, &enorm);
199: PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n", norm);
200: PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n", enorm);
202: /* free spaces */
203: DMRestoreGlobalVector(da, &x);
204: DMRestoreGlobalVector(da, &xx);
205: DMRestoreGlobalVector(da, &y);
206: DMRestoreGlobalVector(da, &yy);
207: DMRestoreGlobalVector(da, &z);
208: DMRestoreGlobalVector(da, &zz);
210: PetscFinalize();
211: return 0;
212: }