Actual source code: ex228.c
petsc-3.14.6 2021-03-30
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 */
23: PetscInt n=10,N=1; /* FFT dimension params */
24: PetscInt DIM,dim[5];/* FFT params */
26: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
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++){dim[k] = n; N*=n;}
39: PetscPrintf(PETSC_COMM_WORLD, "\n %D dimensions: FFTW on vector of size %D \n",DIM,N);
41: /* create FFTW object */
42: MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);
43: /* create vectors of length N */
44: MatCreateVecsFFTW(A,&x,&y,&z);
46: PetscObjectSetName((PetscObject) x, "Real space vector");
47: PetscObjectSetName((PetscObject) y, "Frequency space vector");
48: PetscObjectSetName((PetscObject) z, "Reconstructed vector");
50: /* Test vector duplication*/
51: VecDuplicate(x,&x1);
52: VecDuplicate(y,&y1);
53: VecDuplicate(z,&z1);
55: /* Set values of space vector x, copy to duplicate */
56: VecSetRandom(x,rdm);
57: VecCopy(x,x1);
59: /* Apply FFTW_FORWARD and FFTW_BACKWARD */
60: MatMult(A,x,y);
61: MatMultTranspose(A,y,z);
63: /* Apply FFTW_FORWARD and FFTW_BACKWARD for duplicate vecs */
64: MatMult(A,x1,y1);
65: MatMultTranspose(A,y1,z1);
67: /* Compare x and z1. FFTW computes an unnormalized DFT, thus z1 = N*x */
68: a = 1.0/(PetscReal)N;
69: VecScale(z1,a);
70: VecAXPY(z1,-1.0,x);
71: VecNorm(z1,NORM_1,&enorm);
72: if (enorm > 1.e-9){PetscPrintf(PETSC_COMM_WORLD," Error norm of |x - z1| %g\n",enorm);}
74: /* free spaces */
75: VecDestroy(&x1);
76: VecDestroy(&y1);
77: VecDestroy(&z1);
78: VecDestroy(&x);
79: VecDestroy(&y);
80: VecDestroy(&z);
81: MatDestroy(&A);
82: }
84: PetscRandomDestroy(&rdm);
85: PetscFinalize();
86: return ierr;
87: }
89: /*TEST
91: build:
92: requires: fftw complex
94: test:
95: suffix: 2
96: nsize : 4
97: args: -mat_fftw_plannerflags FFTW_ESTIMATE -n 16
99: test:
100: suffix: 3
101: nsize : 2
102: args: -mat_fftw_plannerflags FFTW_MEASURE -n 12
104: test:
105: suffix: 4
106: nsize : 2
107: args: -mat_fftw_plannerflags FFTW_PATIENT -n 10
109: test:
110: suffix: 5
111: nsize : 1
112: args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE -n 5
114: TEST*/