1 static char help[] = "Test sequential USFFT interface on a uniform DMDA and compares the result to FFTW\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 #include <petscmat.h>
11 #include <petscdm.h>
12 #include <petscdmda.h>
main(int argc,char ** args)13 int main(int argc, char **args)
14 {
15 typedef enum {
16 RANDOM,
17 CONSTANT,
18 TANH,
19 NUM_FUNCS
20 } FuncType;
21 const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
22 Mat A, AA;
23 PetscMPIInt size;
24 PetscInt N, i, stencil = 1, dof = 1;
25 PetscInt dim[3] = {10, 10, 10}, ndim = 3;
26 Vec coords, x, y, z, xx, yy, zz;
27 PetscReal h[3];
28 PetscScalar s;
29 PetscRandom rdm;
30 PetscReal norm, enorm;
31 PetscInt func;
32 FuncType function = TANH;
33 DM da, coordsda;
34 PetscBool view_x = PETSC_FALSE, view_y = PETSC_FALSE, view_z = PETSC_FALSE;
35
36 PetscFunctionBeginUser;
37 PetscCall(PetscInitialize(&argc, &args, NULL, help));
38 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
39 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
40 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
41 PetscCall(PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL));
42 function = (FuncType)func;
43 PetscOptionsEnd();
44 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_x", &view_x, NULL));
45 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_y", &view_y, NULL));
46 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_z", &view_z, NULL));
47 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dim", dim, &ndim, NULL));
48
49 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));
50 PetscCall(DMSetFromOptions(da));
51 PetscCall(DMSetUp(da));
52
53 /* Coordinates */
54 PetscCall(DMGetCoordinateDM(da, &coordsda));
55 PetscCall(DMGetGlobalVector(coordsda, &coords));
56 PetscCall(PetscObjectSetName((PetscObject)coords, "Grid coordinates"));
57 for (i = 0, N = 1; i < 3; i++) {
58 h[i] = 1.0 / dim[i];
59 PetscScalar *a;
60 PetscCall(VecGetArray(coords, &a));
61 PetscInt j, k, n = 0;
62 for (i = 0; i < 3; ++i) {
63 for (j = 0; j < dim[i]; ++j) {
64 for (k = 0; k < 3; ++k) {
65 a[n] = j * h[i]; /* coordinate along the j-th point in the i-th dimension */
66 ++n;
67 }
68 }
69 }
70 PetscCall(VecRestoreArray(coords, &a));
71 }
72 PetscCall(DMSetCoordinates(da, coords));
73
74 /* Work vectors */
75 PetscCall(DMGetGlobalVector(da, &x));
76 PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
77 PetscCall(DMGetGlobalVector(da, &xx));
78 PetscCall(PetscObjectSetName((PetscObject)xx, "Real space vector"));
79 PetscCall(DMGetGlobalVector(da, &y));
80 PetscCall(PetscObjectSetName((PetscObject)y, "USFFT frequency space vector"));
81 PetscCall(DMGetGlobalVector(da, &yy));
82 PetscCall(PetscObjectSetName((PetscObject)yy, "FFTW frequency space vector"));
83 PetscCall(DMGetGlobalVector(da, &z));
84 PetscCall(PetscObjectSetName((PetscObject)z, "USFFT reconstructed vector"));
85 PetscCall(DMGetGlobalVector(da, &zz));
86 PetscCall(PetscObjectSetName((PetscObject)zz, "FFTW reconstructed vector"));
87
88 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%3-" PetscInt_FMT ": USFFT on vector of "));
89 for (i = 0, N = 1; i < 3; i++) {
90 PetscCall(PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ", i, dim[i]));
91 N *= dim[i];
92 }
93 PetscCall(PetscPrintf(PETSC_COMM_SELF, "; total size %d \n", N));
94
95 if (function == RANDOM) {
96 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
97 PetscCall(PetscRandomSetFromOptions(rdm));
98 PetscCall(VecSetRandom(x, rdm));
99 PetscCall(PetscRandomDestroy(&rdm));
100 } else if (function == CONSTANT) {
101 PetscCall(VecSet(x, 1.0));
102 } else if (function == TANH) {
103 PetscScalar *a;
104 PetscCall(VecGetArray(x, &a));
105 PetscInt j, k = 0;
106 for (i = 0; i < 3; ++i) {
107 for (j = 0; j < dim[i]; ++j) {
108 a[k] = tanh((j - dim[i] / 2.0) * (10.0 / dim[i]));
109 ++k;
110 }
111 }
112 PetscCall(VecRestoreArray(x, &a));
113 }
114 if (view_x) PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
115 PetscCall(VecCopy(x, xx));
116
117 PetscCall(VecNorm(x, NORM_2, &norm));
118 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n", norm));
119
120 /* create USFFT object */
121 PetscCall(MatCreateSeqUSFFT(coords, da, &A));
122 /* create FFTW object */
123 PetscCall(MatCreateSeqFFTW(PETSC_COMM_SELF, 3, dim, &AA));
124
125 /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
126 PetscCall(MatMult(A, x, z));
127 PetscCall(MatMult(AA, xx, zz));
128 /* Now apply USFFT and FFTW forward several (3) times */
129 for (i = 0; i < 3; ++i) {
130 PetscCall(MatMult(A, x, y));
131 PetscCall(MatMult(AA, xx, yy));
132 PetscCall(MatMultTranspose(A, y, z));
133 PetscCall(MatMultTranspose(AA, yy, zz));
134 }
135
136 if (view_y) {
137 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "y = \n"));
138 PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
139 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "yy = \n"));
140 PetscCall(VecView(yy, PETSC_VIEWER_STDOUT_WORLD));
141 }
142
143 if (view_z) {
144 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "z = \n"));
145 PetscCall(VecView(z, PETSC_VIEWER_STDOUT_WORLD));
146 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "zz = \n"));
147 PetscCall(VecView(zz, PETSC_VIEWER_STDOUT_WORLD));
148 }
149
150 /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
151 s = 1.0 / (PetscReal)N;
152 PetscCall(VecScale(z, s));
153 PetscCall(VecAXPY(x, -1.0, z));
154 PetscCall(VecNorm(x, NORM_1, &enorm));
155 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n", enorm));
156
157 /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
158 s = 1.0 / (PetscReal)N;
159 PetscCall(VecScale(zz, s));
160 PetscCall(VecAXPY(xx, -1.0, zz));
161 PetscCall(VecNorm(xx, NORM_1, &enorm));
162 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n", enorm));
163
164 /* compare y and yy: USFFT and FFTW results*/
165 PetscCall(VecNorm(y, NORM_2, &norm));
166 PetscCall(VecAXPY(y, -1.0, yy));
167 PetscCall(VecNorm(y, NORM_1, &enorm));
168 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n", norm));
169 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n", enorm));
170
171 /* compare z and zz: USFFT and FFTW results*/
172 PetscCall(VecNorm(z, NORM_2, &norm));
173 PetscCall(VecAXPY(z, -1.0, zz));
174 PetscCall(VecNorm(z, NORM_1, &enorm));
175 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n", norm));
176 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n", enorm));
177
178 /* free spaces */
179 PetscCall(DMRestoreGlobalVector(da, &x));
180 PetscCall(DMRestoreGlobalVector(da, &xx));
181 PetscCall(DMRestoreGlobalVector(da, &y));
182 PetscCall(DMRestoreGlobalVector(da, &yy));
183 PetscCall(DMRestoreGlobalVector(da, &z));
184 PetscCall(DMRestoreGlobalVector(da, &zz));
185 PetscCall(VecDestroy(&coords));
186 PetscCall(DMDestroy(&da));
187 PetscCall(PetscFinalize());
188 return 0;
189 }
190