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 */
26: PetscInitialize(&argc, &args, (char *)0, help);
27: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
29: /* To create random input vector */
30: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
31: 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: PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n", DIM, N);
44: /* create FFTW object */
45: MatCreateFFT(PETSC_COMM_SELF, DIM, dim, MATFFTW, &A);
46: /* create vectors of length N */
47: MatCreateVecsFFTW(A, &x, &y, &z);
49: PetscObjectSetName((PetscObject)x, "Real space vector");
50: PetscObjectSetName((PetscObject)y, "Frequency space vector");
51: PetscObjectSetName((PetscObject)z, "Reconstructed vector");
53: /* Test vector duplication*/
54: VecDuplicate(x, &x1);
55: VecDuplicate(y, &y1);
56: VecDuplicate(z, &z1);
58: /* Set values of space vector x, copy to duplicate */
59: VecSetRandom(x, rdm);
60: VecCopy(x, x1);
62: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
63: MatMult(A, x, y);
64: MatMultTranspose(A, y, z);
66: /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
67: MatMult(A, x1, y1);
68: 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: VecScale(z1, a);
73: VecAXPY(z1, -1.0, x);
74: VecNorm(z1, NORM_1, &enorm);
75: if (enorm > 1.e-9) PetscPrintf(PETSC_COMM_WORLD, " Error norm of |x - z1| %g\n", enorm);
77: /* free spaces */
78: VecDestroy(&x1);
79: VecDestroy(&y1);
80: VecDestroy(&z1);
81: VecDestroy(&x);
82: VecDestroy(&y);
83: VecDestroy(&z);
84: MatDestroy(&A);
85: }
87: PetscRandomDestroy(&rdm);
88: 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*/