Actual source code: ex121.c
1: static char help[] = "Test sequential FFTW convolution\n\n";
3: /*
4: Compiling the code:
5: This code uses the complex numbers, so configure must be given --with-scalar-type=complex to enable this
6: */
7: #include petscmat.h
11: PetscInt main(PetscInt 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,j;
18: Vec w,x,y1,y2,z1,z2;
19: PetscScalar *a, *a2, *a3;
20: PetscScalar s;
21: PetscRandom rdm;
22: PetscReal enorm;
23: PetscInt func;
24: FuncType function = RANDOM;
25: PetscTruth view = PETSC_FALSE;
28: PetscInitialize(&argc,&args,(char *)0,help);
29: #if !defined(PETSC_USE_COMPLEX)
30: SETERRQ(PETSC_ERR_SUP, "This example requires complex numbers");
31: #endif
32: MPI_Comm_size(PETSC_COMM_WORLD, &size);
33: if (size != 1) SETERRQ(PETSC_ERR_SUP, "This is a uniprocessor example only!");
34: PetscOptionsBegin(PETSC_COMM_WORLD, PETSC_NULL, "FFTW Options", "ex112");
35: PetscOptionsEList("-function", "Function type", "ex112", funcNames, NUM_FUNCS, funcNames[function], &func, PETSC_NULL);
36: PetscOptionsTruth("-vec_view_draw", "View the functions", "ex112", view, &view, PETSC_NULL);
37: function = (FuncType) func;
38: PetscOptionsEnd();
40: for(DIM = 0; DIM < ndim; DIM++){
41: dim[DIM] = n; /* size of transformation in DIM-dimension */
42: }
43: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
44: PetscRandomSetFromOptions(rdm);
46: for(DIM = 1; DIM < 5; DIM++){
47: /* create vectors of length N=n^DIM */
48: for(i = 0, N = 1; i < DIM; i++) N *= dim[i];
49: PetscPrintf(PETSC_COMM_SELF, "\n %d-D: FFTW on vector of size %d \n",DIM,N);
50: VecCreateSeq(PETSC_COMM_SELF,N,&x);
51: PetscObjectSetName((PetscObject) x, "Real space vector");
52: VecDuplicate(x,&w);
53: PetscObjectSetName((PetscObject) w, "Window vector");
54: VecDuplicate(x,&y1);
55: PetscObjectSetName((PetscObject) y1, "Frequency space vector");
56: VecDuplicate(x,&y2);
57: PetscObjectSetName((PetscObject) y2, "Frequency space window vector");
58: VecDuplicate(x,&z1);
59: PetscObjectSetName((PetscObject) z1, "Reconstructed convolution");
60: VecDuplicate(x,&z2);
61: PetscObjectSetName((PetscObject) z2, "Real space convolution");
62: if (function == RANDOM) {
63: VecSetRandom(x, rdm);
64: } else if (function == CONSTANT) {
65: VecSet(x, 1.0);
66: } else if (function == TANH) {
67: VecGetArray(x, &a);
68: for(i = 0; i < N; ++i) {
69: a[i] = tanh((i - N/2.0)*(10.0/N));
70: }
71: VecRestoreArray(x, &a);
72: }
73: if (view) {VecView(x, PETSC_VIEWER_DRAW_WORLD);}
75: /* Create window function */
76: VecGetArray(w, &a);
77: for(i = 0; i < N; ++i) {
78: // Step Function
79: a[i] = (i > N/4 && i < 3*N/4)? 1.0: 0.0;
80: // Delta Function
81: //a[i] = (i == N/2)? 1.0: 0.0;
82: }
83: VecRestoreArray(w, &a);
84: if (view) {VecView(w, PETSC_VIEWER_DRAW_WORLD);}
86: /* create FFTW object */
87: //MatCreateSeqFFTW(PETSC_COMM_SELF,DIM,dim,&A);
88: MatCreateSeqFFTW(PETSC_COMM_SELF,1,&N,&A);
90: /* Convolve x with w*/
91: MatMult(A,x,y1);
92: MatMult(A,w,y2);
93: VecPointwiseMult(y1, y1, y2);
94: if (view && i == 0) {VecView(y1, PETSC_VIEWER_DRAW_WORLD);}
95: MatMultTranspose(A,y1,z1);
96:
97: /* Compute the real space convolution */
98: VecGetArray(x, &a);
99: VecGetArray(w, &a2);
100: VecGetArray(z2, &a3);
101: for(i = 0; i < N; ++i) {
102: /* PetscInt checkInd = (i > N/2-1)? i-N/2: i+N/2;*/
104: //if (!(i%100)) PetscPrintf(PETSC_COMM_WORLD, "Finished convolution row %d\n", i);
105: a3[i] = 0.0;
106: for(j = -N/2+1; j < N/2; ++j) {
107: PetscInt xpInd = (j < 0)? N+j: j;
108: PetscInt diffInd = (i-j < 0)? N-(j-i): (i-j > N-1)? i-j-N: i-j;
110: a3[i] += a[xpInd]*a2[diffInd];
111: }
112: //if (PetscAbsScalar(a3[i]) > PetscAbsScalar(a[checkInd])+0.1) PetscPrintf(PETSC_COMM_WORLD, "Invalid convolution at row %d\n", i);
113: }
114: VecRestoreArray(x, &a);
115: VecRestoreArray(w, &a2);
116: VecRestoreArray(z2, &a3);
118: /* compare z1 and z2. FFTW computes an unnormalized DFT, thus z1 = N*z2 */
119: s = 1.0/(PetscReal)N;
120: VecScale(z1,s);
121: if (view) {VecView(z1, PETSC_VIEWER_DRAW_WORLD);}
122: if (view) {VecView(z2, PETSC_VIEWER_DRAW_WORLD);}
123: VecAXPY(z1,-1.0,z2);
124: VecNorm(z1,NORM_1,&enorm);
125: if (enorm > 1.e-11){
126: PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %A\n",enorm);
127: }
129: /* free spaces */
130: VecDestroy(x);
131: VecDestroy(y1);
132: VecDestroy(y2);
133: VecDestroy(z1);
134: VecDestroy(z2);
135: VecDestroy(w);
136: MatDestroy(A);
137: }
138: PetscRandomDestroy(rdm);
139: PetscFinalize();
140: return 0;
141: }