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";
2
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.
7
8 */
9
10 #define DOF 3
11
12 #include <petscmat.h>
13 #include <petscdm.h>
14 #include <petscdmda.h>
main(int argc,char ** args)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;
38
39 PetscFunctionBeginUser;
40 PetscCall(PetscInitialize(&argc, &args, NULL, help));
41 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
42 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
43 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
44 PetscCall(PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL));
45 function = (FuncType)func;
46 PetscOptionsEnd();
47 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_x", &view_x, NULL));
48 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_y", &view_y, NULL));
49 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_z", &view_z, NULL));
50 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dim", dim, &ndim, NULL));
51
52 /* DMDA with the correct fiber dimension */
53 PetscCall(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 PetscCall(DMSetFromOptions(da));
55 PetscCall(DMSetUp(da));
56 /* DMDA with fiber dimension 1 for split fields */
57 PetscCall(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 PetscCall(DMSetFromOptions(da1));
59 PetscCall(DMSetUp(da1));
60
61 /* Coordinates */
62 PetscCall(DMGetCoordinateDM(da, &coordsda));
63 PetscCall(DMGetGlobalVector(coordsda, &coords));
64 PetscCall(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 PetscCall(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 PetscCall(VecRestoreArray(coords, &a));
79 }
80 PetscCall(DMSetCoordinates(da, coords));
81 PetscCall(VecDestroy(&coords));
82
83 /* Work vectors */
84 PetscCall(DMGetGlobalVector(da, &x));
85 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
86 PetscCall(DMGetGlobalVector(da, &xx));
87 PetscCall(PetscObjectSetName((PetscObject)xx, "Real space vector"));
88 PetscCall(DMGetGlobalVector(da, &y));
89 PetscCall(PetscObjectSetName((PetscObject)y, "USFFT frequency space vector"));
90 PetscCall(DMGetGlobalVector(da, &yy));
91 PetscCall(PetscObjectSetName((PetscObject)yy, "FFTW frequency space vector"));
92 PetscCall(DMGetGlobalVector(da, &z));
93 PetscCall(PetscObjectSetName((PetscObject)z, "USFFT reconstructed vector"));
94 PetscCall(DMGetGlobalVector(da, &zz));
95 PetscCall(PetscObjectSetName((PetscObject)zz, "FFTW reconstructed vector"));
96 /* Split vectors for FFTW */
97 for (ii = 0; ii < 3; ++ii) {
98 PetscCall(DMGetGlobalVector(da1, &xxsplit[ii]));
99 PetscCall(PetscObjectSetName((PetscObject)xxsplit[ii], "Real space split vector"));
100 PetscCall(DMGetGlobalVector(da1, &yysplit[ii]));
101 PetscCall(PetscObjectSetName((PetscObject)yysplit[ii], "FFTW frequency space split vector"));
102 PetscCall(DMGetGlobalVector(da1, &zzsplit[ii]));
103 PetscCall(PetscObjectSetName((PetscObject)zzsplit[ii], "FFTW reconstructed split vector"));
104 }
105
106 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%3-" PetscInt_FMT ": USFFT on vector of "));
107 for (i = 0, N = 1; i < 3; i++) {
108 PetscCall(PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ", i, dim[i]));
109 N *= dim[i];
110 }
111 PetscCall(PetscPrintf(PETSC_COMM_SELF, "; total size %d \n", N));
112
113 if (function == RANDOM) {
114 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
115 PetscCall(PetscRandomSetFromOptions(rdm));
116 PetscCall(VecSetRandom(x, rdm));
117 PetscCall(PetscRandomDestroy(&rdm));
118 } else if (function == CONSTANT) {
119 PetscCall(VecSet(x, 1.0));
120 } else if (function == TANH) {
121 PetscScalar *a;
122 PetscCall(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 PetscCall(VecRestoreArray(x, &a));
131 }
132 if (view_x) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
133 PetscCall(VecCopy(x, xx));
134 /* Split xx */
135 PetscCall(VecStrideGatherAll(xx, xxsplit, INSERT_VALUES)); /*YES! 'Gather' means 'split' (or maybe 'scatter'?)! */
136
137 PetscCall(VecNorm(x, NORM_2, &norm));
138 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n", norm));
139
140 /* create USFFT object */
141 PetscCall(MatCreateSeqUSFFT(da, da, &A));
142 /* create FFTW object */
143 PetscCall(MatCreateSeqFFTW(PETSC_COMM_SELF, 3, dim, &AA));
144
145 /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
146 PetscCall(MatMult(A, x, z));
147 for (ii = 0; ii < 3; ++ii) PetscCall(MatMult(AA, xxsplit[ii], zzsplit[ii]));
148 /* Now apply USFFT and FFTW forward several (3) times */
149 for (i = 0; i < 3; ++i) {
150 PetscCall(MatMult(A, x, y));
151 for (ii = 0; ii < 3; ++ii) PetscCall(MatMult(AA, xxsplit[ii], yysplit[ii]));
152 PetscCall(MatMultTranspose(A, y, z));
153 for (ii = 0; ii < 3; ++ii) PetscCall(MatMult(AA, yysplit[ii], zzsplit[ii]));
154 }
155 /* Unsplit yy */
156 PetscCall(VecStrideScatterAll(yysplit, yy, INSERT_VALUES)); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
157 /* Unsplit zz */
158 PetscCall(VecStrideScatterAll(zzsplit, zz, INSERT_VALUES)); /*YES! 'Scatter' means 'collect' (or maybe 'gather'?)! */
159
160 if (view_y) {
161 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "y = \n"));
162 PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
163 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "yy = \n"));
164 PetscCall(VecView(yy, PETSC_VIEWER_STDOUT_WORLD));
165 }
166
167 if (view_z) {
168 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "z = \n"));
169 PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD));
170 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "zz = \n"));
171 PetscCall(VecView(zz, PETSC_VIEWER_STDOUT_WORLD));
172 }
173
174 /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
175 s = 1.0 / (PetscReal)N;
176 PetscCall(VecScale(z, s));
177 PetscCall(VecAXPY(x, -1.0, z));
178 PetscCall(VecNorm(x, NORM_1, &enorm));
179 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n", enorm));
180
181 /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
182 s = 1.0 / (PetscReal)N;
183 PetscCall(VecScale(zz, s));
184 PetscCall(VecAXPY(xx, -1.0, zz));
185 PetscCall(VecNorm(xx, NORM_1, &enorm));
186 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n", enorm));
187
188 /* compare y and yy: USFFT and FFTW results*/
189 PetscCall(VecNorm(y, NORM_2, &norm));
190 PetscCall(VecAXPY(y, -1.0, yy));
191 PetscCall(VecNorm(y, NORM_1, &enorm));
192 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n", norm));
193 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n", enorm));
194
195 /* compare z and zz: USFFT and FFTW results*/
196 PetscCall(VecNorm(z, NORM_2, &norm));
197 PetscCall(VecAXPY(z, -1.0, zz));
198 PetscCall(VecNorm(z, NORM_1, &enorm));
199 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n", norm));
200 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n", enorm));
201
202 /* free spaces */
203 PetscCall(DMRestoreGlobalVector(da, &x));
204 PetscCall(DMRestoreGlobalVector(da, &xx));
205 PetscCall(DMRestoreGlobalVector(da, &y));
206 PetscCall(DMRestoreGlobalVector(da, &yy));
207 PetscCall(DMRestoreGlobalVector(da, &z));
208 PetscCall(DMRestoreGlobalVector(da, &zz));
209
210 PetscCall(PetscFinalize());
211 return 0;
212 }
213