Actual source code: ex27.c
1: static char help[] = "Test sequential USFFT interface on a uniform DMDA and compares the result to FFTW\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: #include <petscmat.h>
11: #include <petscdm.h>
12: #include <petscdmda.h>
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;
37: PetscInitialize(&argc, &args, (char *)0, help);
38: MPI_Comm_size(PETSC_COMM_WORLD, &size);
40: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "USFFT Options", "ex27");
41: PetscOptionsEList("-function", "Function type", "ex27", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
42: function = (FuncType)func;
43: PetscOptionsEnd();
44: PetscOptionsGetBool(NULL, NULL, "-view_x", &view_x, NULL);
45: PetscOptionsGetBool(NULL, NULL, "-view_y", &view_y, NULL);
46: PetscOptionsGetBool(NULL, NULL, "-view_z", &view_z, NULL);
47: PetscOptionsGetIntArray(NULL, NULL, "-dim", dim, &ndim, NULL);
49: 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: DMSetFromOptions(da);
51: DMSetUp(da);
53: /* Coordinates */
54: DMGetCoordinateDM(da, &coordsda);
55: DMGetGlobalVector(coordsda, &coords);
56: 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: 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: VecRestoreArray(coords, &a);
71: }
72: DMSetCoordinates(da, coords);
74: /* Work vectors */
75: DMGetGlobalVector(da, &x);
76: PetscObjectSetName((PetscObject)x, "Real space vector");
77: DMGetGlobalVector(da, &xx);
78: PetscObjectSetName((PetscObject)xx, "Real space vector");
79: DMGetGlobalVector(da, &y);
80: PetscObjectSetName((PetscObject)y, "USFFT frequency space vector");
81: DMGetGlobalVector(da, &yy);
82: PetscObjectSetName((PetscObject)yy, "FFTW frequency space vector");
83: DMGetGlobalVector(da, &z);
84: PetscObjectSetName((PetscObject)z, "USFFT reconstructed vector");
85: DMGetGlobalVector(da, &zz);
86: PetscObjectSetName((PetscObject)zz, "FFTW reconstructed vector");
88: PetscPrintf(PETSC_COMM_SELF, "%3-" PetscInt_FMT ": USFFT on vector of ");
89: for (i = 0, N = 1; i < 3; i++) {
90: PetscPrintf(PETSC_COMM_SELF, "dim[%d] = %d ", i, dim[i]);
91: N *= dim[i];
92: }
93: PetscPrintf(PETSC_COMM_SELF, "; total size %d \n", N);
95: if (function == RANDOM) {
96: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
97: PetscRandomSetFromOptions(rdm);
98: VecSetRandom(x, rdm);
99: PetscRandomDestroy(&rdm);
100: } else if (function == CONSTANT) {
101: VecSet(x, 1.0);
102: } else if (function == TANH) {
103: PetscScalar *a;
104: 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: VecRestoreArray(x, &a);
113: }
114: if (view_x) VecView(x, PETSC_VIEWER_STDOUT_WORLD);
115: VecCopy(x, xx);
117: VecNorm(x, NORM_2, &norm);
118: PetscPrintf(PETSC_COMM_SELF, "|x|_2 = %g\n", norm);
120: /* create USFFT object */
121: MatCreateSeqUSFFT(coords, da, &A);
122: /* create FFTW object */
123: MatCreateSeqFFTW(PETSC_COMM_SELF, 3, dim, &AA);
125: /* apply USFFT and FFTW FORWARD "preemptively", so the fftw_plans can be reused on different vectors */
126: MatMult(A, x, z);
127: MatMult(AA, xx, zz);
128: /* Now apply USFFT and FFTW forward several (3) times */
129: for (i = 0; i < 3; ++i) {
130: MatMult(A, x, y);
131: MatMult(AA, xx, yy);
132: MatMultTranspose(A, y, z);
133: MatMultTranspose(AA, yy, zz);
134: }
136: if (view_y) {
137: PetscPrintf(PETSC_COMM_WORLD, "y = \n");
138: VecView(y, PETSC_VIEWER_STDOUT_WORLD);
139: PetscPrintf(PETSC_COMM_WORLD, "yy = \n");
140: VecView(yy, PETSC_VIEWER_STDOUT_WORLD);
141: }
143: if (view_z) {
144: PetscPrintf(PETSC_COMM_WORLD, "z = \n");
145: VecView(z, PETSC_VIEWER_STDOUT_WORLD);
146: PetscPrintf(PETSC_COMM_WORLD, "zz = \n");
147: VecView(zz, PETSC_VIEWER_STDOUT_WORLD);
148: }
150: /* compare x and z. USFFT computes an unnormalized DFT, thus z = N*x */
151: s = 1.0 / (PetscReal)N;
152: VecScale(z, s);
153: VecAXPY(x, -1.0, z);
154: VecNorm(x, NORM_1, &enorm);
155: PetscPrintf(PETSC_COMM_SELF, "|x-z| = %g\n", enorm);
157: /* compare xx and zz. FFTW computes an unnormalized DFT, thus zz = N*x */
158: s = 1.0 / (PetscReal)N;
159: VecScale(zz, s);
160: VecAXPY(xx, -1.0, zz);
161: VecNorm(xx, NORM_1, &enorm);
162: PetscPrintf(PETSC_COMM_SELF, "|xx-zz| = %g\n", enorm);
164: /* compare y and yy: USFFT and FFTW results*/
165: VecNorm(y, NORM_2, &norm);
166: VecAXPY(y, -1.0, yy);
167: VecNorm(y, NORM_1, &enorm);
168: PetscPrintf(PETSC_COMM_SELF, "|y|_2 = %g\n", norm);
169: PetscPrintf(PETSC_COMM_SELF, "|y-yy| = %g\n", enorm);
171: /* compare z and zz: USFFT and FFTW results*/
172: VecNorm(z, NORM_2, &norm);
173: VecAXPY(z, -1.0, zz);
174: VecNorm(z, NORM_1, &enorm);
175: PetscPrintf(PETSC_COMM_SELF, "|z|_2 = %g\n", norm);
176: PetscPrintf(PETSC_COMM_SELF, "|z-zz| = %g\n", enorm);
178: /* free spaces */
179: DMRestoreGlobalVector(da, &x);
180: DMRestoreGlobalVector(da, &xx);
181: DMRestoreGlobalVector(da, &y);
182: DMRestoreGlobalVector(da, &yy);
183: DMRestoreGlobalVector(da, &z);
184: DMRestoreGlobalVector(da, &zz);
185: VecDestroy(&coords);
186: DMDestroy(&da);
187: PetscFinalize();
188: return 0;
189: }