Actual source code: ex228.c
1: static char help[] = "Test duplication/destruction of FFTW vecs \n\n";
3: /*
4: Compiling the code:
5: This code uses the FFTW interface.
6: Use one of the options below to configure:
7: --with-fftw-dir=/.... or --download-fftw
8: Usage:
9: mpiexec -np <np> ./ex228
10: */
12: #include <petscmat.h>
13: int main(int argc, char **args)
14: {
15: Mat A; /* FFT Matrix */
16: Vec x, y, z; /* Work vectors */
17: Vec x1, y1, z1; /* Duplicate vectors */
18: PetscRandom rdm; /* for creating random input */
19: PetscScalar a; /* used to scale output */
20: PetscReal enorm; /* norm for sanity check */
21: PetscInt n = 10, N = 1; /* FFT dimension params */
22: PetscInt DIM, dim[5]; /* FFT params */
24: PetscFunctionBeginUser;
25: PetscCall(PetscInitialize(&argc, &args, NULL, help));
26: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
28: /* To create random input vector */
29: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
30: PetscCall(PetscRandomSetFromOptions(rdm));
32: /* Iterate over dimensions, use PETSc-FFTW interface */
33: for (PetscInt i = 1; i < 5; i++) {
34: DIM = i;
35: N = 1;
36: for (PetscInt k = 0; k < i; k++) {
37: dim[k] = n;
38: N *= n;
39: }
41: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N));
43: /* create FFTW object */
44: PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A));
45: /* create vectors of length N */
46: PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
48: PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
49: PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
50: PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
52: /* Test vector duplication*/
53: PetscCall(VecDuplicate(x, &x1));
54: PetscCall(VecDuplicate(y, &y1));
55: PetscCall(VecDuplicate(z, &z1));
57: /* Set values of space vector x, copy to duplicate */
58: PetscCall(VecSetRandom(x, rdm));
59: PetscCall(VecCopy(x, x1));
61: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
62: PetscCall(MatMult(A, x, y));
63: PetscCall(MatMultTranspose(A, y, z));
65: /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
66: PetscCall(MatMult(A, x1, y1));
67: PetscCall(MatMultTranspose(A, y1, z1));
69: /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
70: a = 1.0 / (PetscReal)N;
71: PetscCall(VecScale(z1, a));
72: PetscCall(VecAXPY(z1, -1.0, x));
73: PetscCall(VecNorm(z1, NORM_1, &enorm));
74: if (enorm > 1.e-9) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z1| %g dimension %" PetscInt_FMT "\n", enorm, i));
76: /* free spaces */
77: PetscCall(VecDestroy(&x1));
78: PetscCall(VecDestroy(&y1));
79: PetscCall(VecDestroy(&z1));
80: PetscCall(VecDestroy(&x));
81: PetscCall(VecDestroy(&y));
82: PetscCall(VecDestroy(&z));
83: PetscCall(MatDestroy(&A));
84: }
86: PetscCall(PetscRandomDestroy(&rdm));
87: PetscCall(PetscFinalize());
88: return 0;
89: }
91: /*TEST
93: build:
94: requires: fftw complex
96: test:
97: suffix: 2
98: nsize : 4
99: args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16
101: test:
102: suffix: 3
103: nsize : 2
104: args: -mat_fftw_plannerflags FFTW_MEASURE -n 12
106: test:
107: suffix: 4
108: nsize : 2
109: args: -mat_fftw_plannerflags FFTW_PATIENT -n 10
111: test:
112: suffix: 5
113: nsize : 1
114: args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5
116: TEST*/