1
2
3
4
5
6
7
8
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <time.h>
12 #include <fftw3.h>
13 #include <complex.h>
14 #include <math.h>
15 #include <float.h>
16 #include <assert.h>
17 #include <stdint.h>
18 #include <inttypes.h>
19
20
21 #define PI acos(-1)
22 #define MILLION 1000000
23 #define BILLION 1000000000L
24 #define NUMBER_OF_LENGTHS 7
25 #define NUMBER_OF_ITERATIONS 50
26 #define OFFSET 6
27
28
29 double randf(void)
30 {
31 double low=-1.0, high=1.0;
32 return (rand()/(double)(RAND_MAX))*abs(low-high)+low;
33 }
34
35
36 double get_op_time()
37 {
38 complex double a,b;
39 struct timespec start,end;
40 double ns, sum = 0.;
41 int i;
42 for (i=0; i<MILLION; i++)
43 {
44 a = randf()+randf()*I;
45 b = randf()+randf()*I;
46 clock_gettime(CLOCK_MONOTONIC, &start);
47 a*b;
48 clock_gettime(CLOCK_MONOTONIC, &end);
49 ns = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
50
51 sum = sum+ns;
52
53 }
54 return (sum/MILLION);
55 }
56
57
58 int n_close(double A, double B, double maxDiff, double maxRelDiff)
59 {
60 float diff = fabs(A - B);
61 if (diff <= maxDiff)
62 {return 0;}
63 A = fabs(A);
64 B = fabs(B);
65 double largest = (double)(B > A) ? B : A;
66 if (diff <= largest * maxRelDiff)
67 {return 0;}
68 assert(0);
69 }
70
71
72 void bracewell_buneman(complex double *x, int m, int N)
73 {
74 int M=1, i, T, k, r[N], muplus, q, p, upper_range, nu;
75 int nprime, rprime, rprimeprime;
76 complex double temp;
77 muplus = (m+1)>>1;
78 for (i=0; i < N; i++)
79 {r[i] = 0;}
80 upper_range = muplus + 1;
81 for (nu=1; nu < upper_range; nu++)
82 {
83 for (k=0; k<M; k++)
84 {
85 T = r[k] << 1;
86 r[k] = T;
87 r[k+M] = T + 1;
88 }
89 M = M << 1;
90 }
91 if (m & 0x01) {M = M >> 1;};
92 for (q=1; q < M; q++)
93 {
94 nprime = q - M;
95 rprimeprime = r[q] * M;
96 for (p=0; p < r[q]; p++)
97 {
98 nprime = nprime + M;
99 rprime = rprimeprime + r[p];
100 temp = x[nprime];
101 x[nprime] = x[rprime];
102 x[rprime] = temp;
103 }
104 }
105 return;
106 }
107
108
109 void recursion_dif (int BaseE, int N, complex double *x, complex double *twiddle, int tss)
110 {
111 int Nprime=0, BaseO=0, n=0, twiddle_factor;
112 complex double e, o;
113 if (N==1) return;
114 else
115 {
116 Nprime = N/2;
117 BaseO = BaseE + Nprime;
118 for (n=0; n < Nprime; n++)
119 {
120 twiddle_factor = n*tss;
121 e = x[BaseE+n] + x[BaseO+n];
122 o = (x[BaseE+n] - x[BaseO+n]) * twiddle[twiddle_factor];
123 x[BaseE+n] = e;
124 x[BaseO+n] = o;
125 }
126 tss = 2*tss;
127 recursion_dif(BaseE, Nprime, x, twiddle, tss);
128 recursion_dif(BaseO, Nprime, x, twiddle, tss);
129 }
130 }
131
132
133 void iteration_dif(int N, complex double *x2, complex double *twiddle, int p)
134 {
135 int Nprime = 0, BaseE = 0, BaseO = 0, n = 0;
136 int twiddle_factor = 0, Bp = 1, tss = 1, P=0, b=0;
137 complex double e, o;
138 for (P=0; P<p; P++)
139 {
140 Nprime = N/2;
141 BaseE = 0;
142 for (b=0; b<Bp; b++)
143 {
144 BaseO = BaseE+Nprime;
145 for (n=0; n<Nprime; n++)
146 {
147 twiddle_factor = n*tss;
148 e = x2[BaseE+n]+x2[BaseO+n];
149 o = (x2[BaseE+n]-x2[BaseO+n])*twiddle[twiddle_factor];
150 x2[BaseE+n] = e;
151 x2[BaseO+n] = o;
152 }
153 BaseE = BaseE+N;
154 }
155 Bp = Bp*2;
156 N = N/2;
157 tss = 2*tss;
158 }
159 }
160
161
162 int main()
163 {
164 int p;
165 int N;
166 int n, u, c, k, i, nk, tss=1, baseE=0;
167 double dft_mean[NUMBER_OF_LENGTHS], difr_mean[NUMBER_OF_LENGTHS];
168 double difi_mean[NUMBER_OF_LENGTHS], fftw3_mean[NUMBER_OF_LENGTHS];
169 uint64_t dft_mean_int[NUMBER_OF_LENGTHS], difr_mean_int[NUMBER_OF_LENGTHS];
170 uint64_t difi_mean_int[NUMBER_OF_LENGTHS], fftw3_mean_int[NUMBER_OF_LENGTHS];
171 double dft_sum, difr_sum, difi_sum, fftw3_sum;
172 double dft_diff, difr_diff, difi_diff, fftw3_diff;
173 struct timespec start, end;
174 FILE *f = fopen("file.txt","w");
175 double complex_mult_time = get_op_time();
176 double iterations_complex_mult_time_product = NUMBER_OF_ITERATIONS * complex_mult_time;
177 complex double *x, *dft, *twiddle, *complex_exponential;
178 fftw_complex *in, *out;
179 fftw_plan p2;
180 for (u=0; u<NUMBER_OF_LENGTHS; u++)
181 {
182 printf(" \n%d sequence length(s) left to compute performance...\n",NUMBER_OF_LENGTHS-u);
183 dft_sum = 0.;
184 difr_sum = 0.;
185 difi_sum = 0.;
186 fftw3_sum = 0.;
187 dft_diff = 0.;
188 difr_diff = 0.;
189 difi_diff = 0.;
190 fftw3_diff = 0.;
191 p = u + OFFSET;
192 N = (int)powf(2.0, (float)(p));
193 dft = malloc(sizeof(complex double) * N);
194 x = malloc(sizeof(complex double) * N);
195 complex_exponential = malloc(sizeof(complex double) * N);
196 twiddle = malloc(sizeof(complex double) * N/2);
197 for ( nk = 0; nk < N; nk++ )
198 complex_exponential[nk]=cexp(-2*PI*I*nk/N);
199 for (c=0; c<NUMBER_OF_ITERATIONS; c++)
200 {
201 srand(1);
202 for ( i = 0; i < N; i++ )
203 x[i] = randf()+randf()*I;
204 for ( i = 0; i < N; i++ )
205 dft[i] = 0.0 + 0.0*I;
206 clock_gettime(CLOCK_MONOTONIC, &start);
207 for (k=0; k<N; k++)
208 for (n=0; n<N; n++)
209 dft[k] = dft[k] + (x[n] * complex_exponential[(n*k)%N]);
210 clock_gettime(CLOCK_MONOTONIC, &end);
211 dft_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
212 dft_sum = dft_sum + dft_diff;
213 for ( i = 0; i < N/2; i++ )
214 twiddle[i]=cexp(-2*PI*I*i/N);
215 clock_gettime(CLOCK_MONOTONIC, &start);
216 iteration_dif(N,x,twiddle,p);
217 bracewell_buneman(x,p,N);
218 clock_gettime(CLOCK_MONOTONIC, &end);
219 difi_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
220 difi_sum = difi_sum + difi_diff;
221 for (i=0;i<N;i++)
222 {
223 n_close(creal(x[i])/N,creal(dft[i])/N,0.01,DBL_EPSILON);
224 n_close(cimag(x[i])/N,cimag(dft[i])/N,0.01,DBL_EPSILON);
225 }
226 srand(1);
227 for ( i = 0; i < N; i++ )
228 x[i] = randf()+randf()*I;
229 clock_gettime(CLOCK_MONOTONIC, &start);
230 recursion_dif(baseE,N,x,twiddle,tss);
231 bracewell_buneman(x,p,N);
232 clock_gettime(CLOCK_MONOTONIC, &end);
233 difr_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
234 difr_sum = difr_sum + difr_diff;
235 for (i=0;i<N;i++)
236 {
237 n_close(creal(x[i])/N,creal(dft[i])/N,0.01,DBL_EPSILON);
238 n_close(cimag(x[i])/N,cimag(dft[i])/N,0.01,DBL_EPSILON);
239 }
240 in = fftw_malloc(sizeof(fftw_complex)*N);
241 out = fftw_malloc(sizeof(fftw_complex)*N );
242 srand(1);
243 for ( i = 0; i < N; i++ )
244 {
245 x[i] = randf()+randf()*I;
246 in[i][0] = creal(x[i]);
247 in[i][1] = cimag(x[i]);
248 }
249 p2 = fftw_plan_dft_1d ( N, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
250 clock_gettime(CLOCK_MONOTONIC, &start);
251 fftw_execute(p2);
252 clock_gettime(CLOCK_MONOTONIC,&end);
253 fftw3_diff = BILLION * (end.tv_sec - start.tv_sec) + end.tv_nsec - start.tv_nsec;
254 fftw3_sum = fftw3_sum + fftw3_diff;
255 fftw_destroy_plan(p2);
256 }
257 free(x);
258 free(dft);
259 free(twiddle);
260 free(complex_exponential);
261 fftw_free(in);
262 fftw_free(out);
263 dft_mean[u] = dft_sum/iterations_complex_mult_time_product;
264 difr_mean[u] = difr_sum/iterations_complex_mult_time_product;
265 difi_mean[u] = difi_sum/iterations_complex_mult_time_product;
266 fftw3_mean[u] = fftw3_sum/iterations_complex_mult_time_product;
267 dft_mean_int[u] = (uint64_t)round(dft_mean[u]);
268 difr_mean_int[u] = (uint64_t)round(difr_mean[u]);
269 difi_mean_int[u] = (uint64_t)round(difi_mean[u]);
270 fftw3_mean_int[u] = (uint64_t)round(fftw3_mean[u]);
271
272 }
273 for ( i = 0; i < NUMBER_OF_LENGTHS; i++ )
274 fprintf (f, " %3d " "%"PRIu64 "\n", (int)powf(2.0, (float)(i+OFFSET)) , difr_mean_int[i] );
275 for ( i = 0; i< NUMBER_OF_LENGTHS; i++ )
276 fprintf (f, " %3d " "%"PRIu64 "\n", (int)powf(2.0, (float)(i+OFFSET)) , difi_mean_int[i] );
277 for ( i = 0; i < NUMBER_OF_LENGTHS; i++ )
278 fprintf (f, " %3d " "%"PRIu64 "\n", (int)powf(2.0, (float)(i+OFFSET)) , dft_mean_int[i] );
279 for ( i = 0; i < NUMBER_OF_LENGTHS; i++ )
280 fprintf (f, " %3d " "%"PRIu64 "\n", (int)powf(2.0, (float)(i+OFFSET)) , fftw3_mean_int[i] );
281
282 fclose(f);
283 system("python plot11.py");
284 exit(0);
285 }