1 static char help[] = "Test duplication/destruction of FFTW vecs \n\n";
2
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 */
11
12 #include <petscmat.h>
main(int argc,char ** args)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 PetscErrorCode ierr; /* to catch bugs, if any */
23 PetscInt n=10,N=1; /* FFT dimension params */
24 PetscInt DIM,dim[5];/* FFT params */
25
26 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
27 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
28
29 /* To create random input vector */
30 ierr = PetscRandomCreate(PETSC_COMM_SELF, &rdm);CHKERRQ(ierr);
31 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
32
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++){dim[k] = n; N*=n;}
38
39 ierr = PetscPrintf(PETSC_COMM_WORLD, "\n %D dimensions: FFTW on vector of size %D \n",DIM,N);CHKERRQ(ierr);
40
41 /* create FFTW object */
42 ierr = MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);CHKERRQ(ierr);
43 /* create vectors of length N */
44 ierr = MatCreateVecsFFTW(A,&x,&y,&z);CHKERRQ(ierr);
45
46 ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr);
47 ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr);
48 ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr);
49
50 /* Test vector duplication*/
51 ierr = VecDuplicate(x,&x1);CHKERRQ(ierr);
52 ierr = VecDuplicate(y,&y1);CHKERRQ(ierr);
53 ierr = VecDuplicate(z,&z1);CHKERRQ(ierr);
54
55 /* Set values of space vector x, copy to duplicate */
56 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr);
57 ierr = VecCopy(x,x1);CHKERRQ(ierr);
58
59 /* Apply FFTW_FORWARD and FFTW_BACKWARD */
60 ierr = MatMult(A,x,y);CHKERRQ(ierr);
61 ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
62
63 /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
64 ierr = MatMult(A,x1,y1);CHKERRQ(ierr);
65 ierr = MatMultTranspose(A,y1,z1);CHKERRQ(ierr);
66
67 /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
68 a = 1.0/(PetscReal)N;
69 ierr = VecScale(z1,a);CHKERRQ(ierr);
70 ierr = VecAXPY(z1,-1.0,x);CHKERRQ(ierr);
71 ierr = VecNorm(z1,NORM_1,&enorm);CHKERRQ(ierr);
72 if (enorm > 1.e-9){ierr = PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z1| %g\n",enorm);CHKERRQ(ierr);}
73
74 /* free spaces */
75 ierr = VecDestroy(&x1);CHKERRQ(ierr);
76 ierr = VecDestroy(&y1);CHKERRQ(ierr);
77 ierr = VecDestroy(&z1);CHKERRQ(ierr);
78 ierr = VecDestroy(&x);CHKERRQ(ierr);
79 ierr = VecDestroy(&y);CHKERRQ(ierr);
80 ierr = VecDestroy(&z);CHKERRQ(ierr);
81 ierr = MatDestroy(&A);CHKERRQ(ierr);
82 }
83
84 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
85 ierr = PetscFinalize();
86 return ierr;
87 }
88
89 /*TEST
90
91 build:
92 requires: fftw complex
93
94 test:
95 suffix: 2
96 nsize : 4
97 args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16
98
99 test:
100 suffix: 3
101 nsize : 2
102 args: -mat_fftw_plannerflags FFTW_MEASURE -n 12
103
104 test:
105 suffix: 4
106 nsize : 2
107 args: -mat_fftw_plannerflags FFTW_PATIENT -n 10
108
109 test:
110 suffix: 5
111 nsize : 1
112 args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5
113
114 TEST*/
115