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: PetscInt i, k; /* for iterating over dimensions */
19: PetscRandom rdm; /* for creating random input */
20: PetscScalar a; /* used to scale output */
21: PetscReal enorm; /* norm for sanity check */
22: PetscInt n = 10, N = 1; /* FFT dimension params */
23: PetscInt DIM, dim[5]; /* FFT params */
25: PetscFunctionBeginUser;
26: PetscCall(PetscInitialize(&argc, &args, NULL, help));
27: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
29: /* To create random input vector */
30: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
31: PetscCall(PetscRandomSetFromOptions(rdm));
33: /* Iterate over dimensions, use PETSc-FFTW interface */
34: for (i = 1; i < 5; i++) {
35: DIM = i;
36: N = 1;
37: for (k = 0; k < i; k++) {
38: dim[k] = n;
39: N *= n;
40: }
42: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N));
44: /* create FFTW object */
45: PetscCall(MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A));
46: /* create vectors of length N */
47: PetscCall(MatCreateVecsFFTW(A, &x, &y, &z));
49: PetscCall(PetscObjectSetName((PetscObject)x, "Real space vector"));
50: PetscCall(PetscObjectSetName((PetscObject)y, "Frequency space vector"));
51: PetscCall(PetscObjectSetName((PetscObject)z, "Reconstructed vector"));
53: /* Test vector duplication*/
54: PetscCall(VecDuplicate(x, &x1));
55: PetscCall(VecDuplicate(y, &y1));
56: PetscCall(VecDuplicate(z, &z1));
58: /* Set values of space vector x, copy to duplicate */
59: PetscCall(VecSetRandom(x, rdm));
60: PetscCall(VecCopy(x, x1));
62: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
63: PetscCall(MatMult(A, x, y));
64: PetscCall(MatMultTranspose(A, y, z));
66: /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
67: PetscCall(MatMult(A, x1, y1));
68: PetscCall(MatMultTranspose(A, y1, z1));
70: /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
71: a = 1.0 / (PetscReal)N;
72: PetscCall(VecScale(z1, a));
73: PetscCall(VecAXPY(z1, -1.0, x));
74: PetscCall(VecNorm(z1, NORM_1, &enorm));
75: if (enorm > 1.e-9) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z1| %g\n", enorm));
77: /* free spaces */
78: PetscCall(VecDestroy(&x1));
79: PetscCall(VecDestroy(&y1));
80: PetscCall(VecDestroy(&z1));
81: PetscCall(VecDestroy(&x));
82: PetscCall(VecDestroy(&y));
83: PetscCall(VecDestroy(&z));
84: PetscCall(MatDestroy(&A));
85: }
87: PetscCall(PetscRandomDestroy(&rdm));
88: PetscCall(PetscFinalize());
89: return 0;
90: }
92: /*TEST
94: build:
95: requires: fftw complex
97: test:
98: suffix: 2
99: nsize : 4
100: args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16
102: test:
103: suffix: 3
104: nsize : 2
105: args: -mat_fftw_plannerflags FFTW_MEASURE -n 12
107: test:
108: suffix: 4
109: nsize : 2
110: args: -mat_fftw_plannerflags FFTW_PATIENT -n 10
112: test:
113: suffix: 5
114: nsize : 1
115: args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5
117: TEST*/