Actual source code: ex112.c
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
13: PetscInt main(PetscInt argc,char **args)
14: {
15: typedef enum {RANDOM, CONSTANT, TANH, NUM_FUNCS} FuncType;
16: const char *funcNames[NUM_FUNCS] = {"random", "constant", "tanh"};
17: Mat A;
18: PetscMPIInt size;
19: PetscInt n = 10,N,ndim=4,dim[4],DIM,i;
20: Vec x,y,z;
21: PetscScalar s;
22: PetscRandom rdm;
23: PetscReal enorm;
24: PetscInt func;
25: FuncType function = RANDOM;
26: PetscTruth view = PETSC_FALSE;
29: PetscInitialize(&argc,&args,(char *)0,help);
30: #if !defined(PETSC_USE_COMPLEX)
31: SETERRQ(PETSC_ERR_SUP, "This example requires complex numbers");
32: #endif
33: MPI_Comm_size(PETSC_COMM_WORLD, &size);
34: if (size != 1) SETERRQ(PETSC_ERR_SUP, "This is a uniprocessor example only!");
35: PetscOptionsBegin(PETSC_COMM_WORLD, PETSC_NULL, "FFTW Options", "ex112");
36: PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, PETSC_NULL);
37: PetscOptionsTruth("-vec_view_draw", "View the functions", "ex112", view, &view, PETSC_NULL);
38: function = (FuncType) func;
39: PetscOptionsEnd();
41: for(DIM = 0; DIM < ndim; DIM++){
42: dim[DIM] = n; /* size of transformation in DIM-dimension */
43: }
44: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
45: PetscRandomSetFromOptions(rdm);
47: for(DIM = 1; DIM < 5; DIM++){
48: /* create vectors of length N=n^DIM */
49: for(i = 0, N = 1; i < DIM; i++) N *= dim[i];
50: PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
51: VecCreateSeq(PETSC_COMM_SELF,N,&x);
52: PetscObjectSetName((PetscObject) x, "Real space vector");
53: VecDuplicate(x,&y);
54: PetscObjectSetName((PetscObject) y, "Frequency space vector");
55: VecDuplicate(x,&z);
56: PetscObjectSetName((PetscObject) z, "Reconstructed vector");
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;
64: VecGetArray(x, &a);
65: for(i = 0; i < N; ++i) {
66: a[i] = tanh((i - N/2.0)*(10.0/N));
67: }
68: VecRestoreArray(x, &a);
69: }
70: if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}
72: /* create FFTW object */
73: MatCreateSeqFFTW(PETSC_COMM_SELF,DIM,dim,&A);
75: /* apply FFTW_FORWARD several times, so the fftw_plan can be reused on different vectors */
76: MatMult(A,x,z);
77: for (i=0; i<3; i++){
78: MatMult(A,x,y);
79: if (view && i == 0) {VecView(y, PETSC_VIEWER_DRAW_WORLD);}
80: /* apply FFTW_BACKWARD several times */
81: MatMultTranspose(A,y,z);
82: }
83:
84: /* compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
85: s = 1.0/(PetscReal)N;
86: VecScale(z,s);
87: if (view) {VecView(z, PETSC_VIEWER_DRAW_WORLD);}
88: VecAXPY(z,-1.0,x);
89: VecNorm(z,NORM_1,&enorm);
90: if (enorm > 1.e-11){
91: PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %A\n",enorm);
92: }
94: /* free spaces */
95: VecDestroy(x);
96: VecDestroy(y);
97: VecDestroy(z);
98: MatDestroy(A);
99: }
100: PetscRandomDestroy(rdm);
101: PetscFinalize();
102: return 0;
103: }