Actual source code: ex112.c
petsc-3.14.6 2021-03-30
1: static char help[] = "Test sequential FFTW interface \n\n";
3: /*
4: Compiling the code:
5: This code uses the complex numbers version of PETSc, so configure
6: must be run to enable this
8: */
10: #include <petscmat.h>
11: int main(int argc,char **args)
12: {
13: typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
14: const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
15: Mat A;
16: PetscMPIInt size;
17: PetscInt n = 10,N,ndim=4,dim[4],DIM,i;
18: Vec x,y,z;
19: PetscScalar s;
20: PetscRandom rdm;
21: PetscReal enorm, tol = PETSC_SMALL;
22: PetscInt func;
23: FuncType function = RANDOM;
24: PetscBool view = PETSC_FALSE;
27: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
28: MPI_Comm_size(PETSC_COMM_WORLD, &size);
29: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This is a uniprocessor example only!");
30: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex112");
31: PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, NULL);
32: PetscOptionsBool("-vec_view_draw", "View the functions", "ex112", view, &view, NULL);
33: function = (FuncType) func;
34: PetscOptionsEnd();
36: for (DIM = 0; DIM < ndim; DIM++) {
37: dim[DIM] = n; /* size of transformation in DIM-dimension */
38: }
39: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
40: PetscRandomSetFromOptions(rdm);
42: for (DIM = 1; DIM < 5; DIM++) {
43: for (i = 0, N = 1; i < DIM; i++) N *= dim[i];
44: PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
46: /* create FFTW object */
47: MatCreateFFT(PETSC_COMM_SELF,DIM,dim,MATFFTW,&A);
49: /* create vectors of length N=n^DIM */
50: MatCreateVecs(A,&x,&y);
51: MatCreateVecs(A,&z,NULL);
52: PetscObjectSetName((PetscObject) x, "Real space vector");
53: PetscObjectSetName((PetscObject) y, "Frequency space vector");
54: PetscObjectSetName((PetscObject) z, "Reconstructed vector");
56: /* set values of space vector x */
57: if (function == RANDOM) {
58: VecSetRandom(x, rdm);
59: } else if (function == CONSTANT) {
60: VecSet(x, 1.0);
61: } else if (function == TANH) {
62: PetscScalar *a;
63: VecGetArray(x, &a);
64: for (i = 0; i < N; ++i) {
65: a[i] = tanh((i - N/2.0)*(10.0/N));
66: }
67: VecRestoreArray(x, &a);
68: }
69: if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}
71: /* apply FFTW_FORWARD and FFTW_BACKWARD several times on same x, y, and z */
72: for (i=0; i<3; i++) {
73: MatMult(A,x,y);
74: if (view && i == 0) {VecView(y, PETSC_VIEWER_DRAW_WORLD);}
75: MatMultTranspose(A,y,z);
77: /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
78: s = 1.0/(PetscReal)N;
79: VecScale(z,s);
80: if (view && i == 0) {VecView(z, PETSC_VIEWER_DRAW_WORLD);}
81: VecAXPY(z,-1.0,x);
82: VecNorm(z,NORM_1,&enorm);
83: if (enorm > tol) {
84: PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %g\n",(double)enorm);
85: }
86: }
88: /* apply FFTW_FORWARD and FFTW_BACKWARD several times on different x */
89: for (i=0; i<3; i++) {
90: VecDestroy(&x);
91: VecCreateSeq(PETSC_COMM_SELF,N,&x);
92: VecSetRandom(x, rdm);
94: MatMult(A,x,y);
95: MatMultTranspose(A,y,z);
97: /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
98: s = 1.0/(PetscReal)N;
99: VecScale(z,s);
100: if (view && i == 0) {VecView(z, PETSC_VIEWER_DRAW_WORLD);}
101: VecAXPY(z,-1.0,x);
102: VecNorm(z,NORM_1,&enorm);
103: if (enorm > tol) {
104: PetscPrintf(PETSC_COMM_SELF," Error norm of new |x - z| %g\n",(double)enorm);
105: }
106: }
108: /* free spaces */
109: VecDestroy(&x);
110: VecDestroy(&y);
111: VecDestroy(&z);
112: MatDestroy(&A);
113: }
114: PetscRandomDestroy(&rdm);
115: PetscFinalize();
116: return ierr;
117: }
120: /*TEST
122: build:
123: requires: fftw complex
125: test:
126: args: -mat_fftw_plannerflags FFTW_ESTIMATE
127: output_file: output/ex112.out
129: test:
130: suffix: 2
131: args: -mat_fftw_plannerflags FFTW_MEASURE
132: output_file: output/ex112.out
133: requires: !define(PETSC_USE_CXXCOMPLEX)
135: test:
136: suffix: 3
137: args: -mat_fftw_plannerflags FFTW_PATIENT
138: output_file: output/ex112.out
140: test:
141: suffix: 4
142: args: -mat_fftw_plannerflags FFTW_EXHAUSTIVE
143: output_file: output/ex112.out
145: TEST*/