Actual source code: ex112.c
1: static char help[] = "Test sequential FFTW interface \n\n";
3: /*
4: Compiling the code:
5: This code uses the complex numbers version of PETSc, so configure
6: must be run to enable this
8: */
10: #include <petscmat.h>
11: int main(int argc, char **args)
12: {
13: typedef enum {
14: RANDOM,
15: CONSTANT,
16: TANH,
17: NUM_FUNCS
18: } FuncType;
19: const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
20: Mat A;
21: PetscMPIInt size;
22: PetscInt n = 10, N, ndim = 4, dim[4], DIM, i;
23: Vec x, y, z;
24: PetscScalar s;
25: PetscRandom rdm;
26: PetscReal enorm, tol = PETSC_SMALL;
27: PetscInt func;
28: FuncType function = RANDOM;
29: PetscBool view = PETSC_FALSE;
31: PetscFunctionBeginUser;
32: PetscCall(PetscInitialize(&argc, &args, NULL, help));
33: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
34: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
35: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");
36: PetscCall(PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, NULL));
37: PetscCall(PetscOptionsBool("-vec_view_draw", "View the functions", "ex112", view, &view, NULL));
38: function = (FuncType)func;
39: PetscOptionsEnd();
41: for (DIM = 0; DIM < ndim; DIM++) { dim[DIM] = n; /* size of transformation in DIM-dimension */ }
42: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
43: PetscCall(PetscRandomSetFromOptions(rdm));
45: for (DIM = 1; DIM < 5; DIM++) {
46: for (i = 0, N = 1; i < DIM; i++) N *= dim[i];
47: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n", DIM, N));
49: /* create FFTW object */
50: PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A));
52: /* create vectors of length N=n^DIM */
53: PetscCall(MatCreateVecs(A, &x, &y));
54: PetscCall(MatCreateVecs(A, &z, NULL));
55: PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
56: PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
57: PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
59: /* set values of space vector x */
60: if (function == RANDOM) {
61: PetscCall(VecSetRandom(x, rdm));
62: } else if (function == CONSTANT) {
63: PetscCall(VecSet(x, 1.0));
64: } else if (function == TANH) {
65: PetscScalar *a;
66: PetscCall(VecGetArray(x, &a));
67: for (i = 0; i < N; ++i) a[i] = tanh((i - N / 2.0) * (10.0 / N));
68: PetscCall(VecRestoreArray(x, &a));
69: }
70: if (view) PetscCall(VecView(x, PETSC_VIEWER_DRAW_WORLD));
72: /* apply FFTW_FORWARD and FFTW_BACKWARD several times on same x, y, and z */
73: for (i = 0; i < 3; i++) {
74: PetscCall(MatMult(A, x, y));
75: if (view && i == 0) PetscCall(VecView(y, PETSC_VIEWER_DRAW_WORLD));
76: PetscCall(MatMultTranspose(A, y, z));
78: /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
79: s = 1.0 / (PetscReal)N;
80: PetscCall(VecScale(z, s));
81: if (view && i == 0) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD));
82: PetscCall(VecAXPY(z, -1.0, x));
83: PetscCall(VecNorm(z, NORM_1, &enorm));
84: if (enorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of |x - z| %g\n", (double)enorm));
85: }
87: /* apply FFTW_FORWARD and FFTW_BACKWARD several times on different x */
88: for (i = 0; i < 3; i++) {
89: PetscCall(VecDestroy(&x));
90: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
91: PetscCall(VecSetRandom(x, rdm));
93: PetscCall(MatMult(A, x, y));
94: PetscCall(MatMultTranspose(A, y, z));
96: /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
97: s = 1.0 / (PetscReal)N;
98: PetscCall(VecScale(z, s));
99: if (view && i == 0) PetscCall(VecView(z, PETSC_VIEWER_DRAW_WORLD));
100: PetscCall(VecAXPY(z, -1.0, x));
101: PetscCall(VecNorm(z, NORM_1, &enorm));
102: if (enorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Error norm of new |x - z| %g\n", (double)enorm));
103: }
105: /* free spaces */
106: PetscCall(VecDestroy(&x));
107: PetscCall(VecDestroy(&y));
108: PetscCall(VecDestroy(&z));
109: PetscCall(MatDestroy(&A));
110: }
111: PetscCall(PetscRandomDestroy(&rdm));
112: PetscCall(PetscFinalize());
113: return 0;
114: }
116: /*TEST
118: build:
119: requires: fftw complex
121: test:
122: args: -mat_fftw_plannerflags FFTW_ESTIMATE
123: output_file: output/ex112.out
125: test:
126: suffix: 2
127: args: -mat_fftw_plannerflags FFTW_MEASURE
128: output_file: output/ex112.out
129: requires: !defined(PETSC_USE_CXXCOMPLEX)
131: test:
132: suffix: 3
133: args: -mat_fftw_plannerflags FFTW_PATIENT
134: output_file: output/ex112.out
136: test:
137: suffix: 4
138: args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE
139: output_file: output/ex112.out
141: TEST*/