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: PetscInitialize(&argc,&args,(char*)0,help);
26: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
28: /* To create random input vector */
29: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
30: PetscRandomSetFromOptions(rdm);
32: /* Iterate over dimensions, use PETSc-FFTW interface */
33: for (i=1; i<5; i++) {
34: DIM = i;
35: N = 1;
36: for (k=0; k<i; k++){dim[k] = n; N*=n;}
38: PetscPrintf(PETSC_COMM_WORLD, "\n %" PetscInt_FMT " dimensions: FFTW on vector of size %" PetscInt_FMT " \n",DIM,N);
40: /* create FFTW object */
41: MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);
42: /* create vectors of length N */
43: MatCreateVecsFFTW(A,&x,&y,&z);
45: PetscObjectSetName((PetscObject) x, "Real space vector");
46: PetscObjectSetName((PetscObject) y, "Frequency space vector");
47: PetscObjectSetName((PetscObject) z, "Reconstructed vector");
49: /* Test vector duplication*/
50: VecDuplicate(x,&x1);
51: VecDuplicate(y,&y1);
52: VecDuplicate(z,&z1);
54: /* Set values of space vector x, copy to duplicate */
55: VecSetRandom(x,rdm);
56: VecCopy(x,x1);
58: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
59: MatMult(A,x,y);
60: MatMultTranspose(A,y,z);
62: /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
63: MatMult(A,x1,y1);
64: MatMultTranspose(A,y1,z1);
66: /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
67: a = 1.0/(PetscReal)N;
68: VecScale(z1,a);
69: VecAXPY(z1,-1.0,x);
70: VecNorm(z1,NORM_1,&enorm);
71: if (enorm > 1.e-9) PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z1| %g\n",enorm);
73: /* free spaces */
74: VecDestroy(&x1);
75: VecDestroy(&y1);
76: VecDestroy(&z1);
77: VecDestroy(&x);
78: VecDestroy(&y);
79: VecDestroy(&z);
80: MatDestroy(&A);
81: }
83: PetscRandomDestroy(&rdm);
84: PetscFinalize();
85: return 0;
86: }
88: /*TEST
90: build:
91: requires: fftw complex
93: test:
94: suffix: 2
95: nsize : 4
96: args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16
98: test:
99: suffix: 3
100: nsize : 2
101: args: -mat_fftw_plannerflags FFTW_MEASURE -n 12
103: test:
104: suffix: 4
105: nsize : 2
106: args: -mat_fftw_plannerflags FFTW_PATIENT -n 10
108: test:
109: suffix: 5
110: nsize : 1
111: args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5
113: TEST*/