1 | /*************************************************************************/ |
---|
2 | /* */ |
---|
3 | /* Copyright (c) 1994 Stanford University */ |
---|
4 | /* */ |
---|
5 | /* All rights reserved. */ |
---|
6 | /* */ |
---|
7 | /* Permission is given to use, copy, and modify this software for any */ |
---|
8 | /* non-commercial purpose as long as this copyright notice is not */ |
---|
9 | /* removed. All other uses, including redistribution in whole or in */ |
---|
10 | /* part, are forbidden without prior written permission. */ |
---|
11 | /* */ |
---|
12 | /* This software is provided with absolutely no warranty and no */ |
---|
13 | /* support. */ |
---|
14 | /* */ |
---|
15 | /*************************************************************************/ |
---|
16 | |
---|
17 | /////////////////////////////////////////////////////////////////////////// |
---|
18 | // This port of the SPLASH FFT benchmark on the ALMOS-MKH OS has been |
---|
19 | // done by Alain Greiner (august 2018). |
---|
20 | // |
---|
21 | // This application performs the 1D fast Fourier transfom for an array |
---|
22 | // of N complex points, using the Cooley-Tuckey FFT method. |
---|
23 | // The N data points are seen as a 2D array (rootN rows * rootN columns). |
---|
24 | // Each thread handle (rootN / nthreads) rows. |
---|
25 | // The N input data points can be initialised in three different modes: |
---|
26 | // - CONSTANT : all data points have the same [1,0] value |
---|
27 | // - COSIN : data point n has [cos(n/N) , sin(n/N)] values |
---|
28 | // - RANDOM : data points have pseudo random values |
---|
29 | // |
---|
30 | // The main parameters for this generic application are the following: |
---|
31 | // - M : N = 2**M = number of data points / M must be an even number. |
---|
32 | // - T : nthreads = ncores defined by the hardware / must be power of 2. |
---|
33 | // The number of threads cannot be larger than the number of rows. |
---|
34 | // |
---|
35 | // This application uses 3 shared data arrays, that are dynamically |
---|
36 | // allocated and distributed in clusters, with one sub-buffer per cluster: |
---|
37 | // - data[N] contains N input data points, |
---|
38 | // - trans[N] contains N intermediate data points, |
---|
39 | // - twid[N] contains N coefs : exp(2*pi*i*j/N) / i and j in [0,rootN-1] |
---|
40 | // Each sub-buffer contains (N/nclusters) entries, with 2 double per entry. |
---|
41 | // These distributed buffers are allocated and initialised in parallel |
---|
42 | // by the working threads running on core 0 in each cluster. |
---|
43 | // |
---|
44 | // Each working thread allocates also a private coefs[rootN-1] buffer, |
---|
45 | // that contains all coefs required for a rootN points FFT. |
---|
46 | // |
---|
47 | // There is one working thread per core. |
---|
48 | // The actual number of cores and cluster in a given hardware architecture |
---|
49 | // is obtained by the get_config() syscall (x_size, y_size, ncores). |
---|
50 | // The max number of clusters is bounded by (X_MAX * Y_MAX). |
---|
51 | // The max number of cores per cluster is bounded by CORES_MAX. |
---|
52 | // |
---|
53 | // Several configuration parameters can be defined below: |
---|
54 | // - PRINT_ARRAY : Print out complex data points arrays. |
---|
55 | // - CHECK : Perform both FFT and inverse FFT to check output/input. |
---|
56 | // - DEBUG_MAIN : Display intermediate results in main() |
---|
57 | // - DEBUG_FFT1D : Display intermediate results in FFT1D() |
---|
58 | // - DEBUG_ROW : Display intermedite results in FFTrow() |
---|
59 | // |
---|
60 | // Regarding final instrumentation: |
---|
61 | // - the sequencial initialisation time (init_time) is computed |
---|
62 | // by the main thread in the main() function. |
---|
63 | // - The parallel execution time (parallel_time[i]) is computed by each |
---|
64 | // working thread(i) in the work() function. |
---|
65 | // - The synchronisation time related to the barriers (sync_time[i]) |
---|
66 | // is computed by each thread(i) in the work() function. |
---|
67 | // The results are displayed on the TXT terminal, and registered on disk. |
---|
68 | /////////////////////////////////////////////////////////////////////////// |
---|
69 | |
---|
70 | #include <math.h> |
---|
71 | #include <stdio.h> |
---|
72 | #include <stdlib.h> |
---|
73 | #include <fcntl.h> |
---|
74 | #include <unistd.h> |
---|
75 | #include <pthread.h> |
---|
76 | #include <almosmkh.h> |
---|
77 | #include <hal_macros.h> |
---|
78 | |
---|
79 | // constants |
---|
80 | |
---|
81 | #define PI 3.14159265359 |
---|
82 | #define PAGE_SIZE 4096 |
---|
83 | #define X_MAX 16 // max number of clusters in a row |
---|
84 | #define Y_MAX 16 // max number of clusters in a column |
---|
85 | #define CORES_MAX 4 // max number of cores in a cluster |
---|
86 | #define CLUSTERS_MAX X_MAX * Y_MAX |
---|
87 | #define THREADS_MAX CLUSTERS_MAX * CORES_MAX |
---|
88 | #define RANDOM 0 |
---|
89 | #define COSIN 1 |
---|
90 | #define CONSTANT 2 |
---|
91 | |
---|
92 | // parameters |
---|
93 | |
---|
94 | #define DEFAULT_M 14 // 16384 data points |
---|
95 | #define USE_DQT_BARRIER 1 // use DDT barrier if non zero |
---|
96 | #define MODE COSIN // DATA array initialisation mode |
---|
97 | #define CHECK 0 |
---|
98 | #define DEBUG_MAIN 0 // trace main() function (detailed if odd) |
---|
99 | #define DEBUG_WORK 0 // trace work() function (detailed if odd) |
---|
100 | #define DEBUG_FFT1D 0 // trace FFT1D() function (detailed if odd) |
---|
101 | #define DEBUG_ROW 0 // trace FFTRow() function (detailed if odd) |
---|
102 | #define PRINT_ARRAY 0 |
---|
103 | #define DISPLAY_SCHED_AND_VMM 0 // display final VMM state in all clusters |
---|
104 | |
---|
105 | // macro to swap two variables |
---|
106 | #define SWAP(a,b) { double tmp; tmp = a; a = b; b = tmp; } |
---|
107 | |
---|
108 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
109 | // FFT specific global variables |
---|
110 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
111 | |
---|
112 | // work function arguments |
---|
113 | typedef struct work_args_s |
---|
114 | { |
---|
115 | unsigned int tid; // thread continuous index |
---|
116 | unsigned int lid; // core local index |
---|
117 | unsigned int cid; // cluster continuous index |
---|
118 | pthread_barrier_t * parent_barrier; // parent barrier to signal completion |
---|
119 | } |
---|
120 | work_args_t; |
---|
121 | |
---|
122 | unsigned int nthreads; // total number of threads (one thread per core) |
---|
123 | unsigned int nclusters; // total number of clusters |
---|
124 | unsigned int M = DEFAULT_M; // log2(number of points) |
---|
125 | unsigned int N; // number of points (N = 2^M) |
---|
126 | unsigned int rootN; // rootN = 2^M/2 |
---|
127 | unsigned int rows_per_thread; // number of data "rows" handled by a single thread |
---|
128 | unsigned int points_per_cluster; // number of data points per cluster |
---|
129 | |
---|
130 | // arrays of pointers on distributed buffers (one sub-buffer per cluster) |
---|
131 | double * data[CLUSTERS_MAX]; // original time-domain data |
---|
132 | double * trans[CLUSTERS_MAX]; // used as auxiliary space for transpose |
---|
133 | double * twid[CLUSTERS_MAX]; // twiddle factor : exp(-2iPI*k*n/N) |
---|
134 | double * bloup[CLUSTERS_MAX]; // used as auxiliary space for DFT |
---|
135 | |
---|
136 | // instrumentation counters |
---|
137 | unsigned int parallel_time[THREADS_MAX]; // total computation time (per thread) |
---|
138 | unsigned int sync_time[THREADS_MAX]; // cumulated waiting time in barriers (per thread) |
---|
139 | unsigned int init_time; // initialisation time (in main) |
---|
140 | |
---|
141 | // synchronisation barrier (all threads) |
---|
142 | pthread_barrier_t barrier; |
---|
143 | pthread_barrierattr_t barrier_attr; |
---|
144 | |
---|
145 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
146 | // Global variables required by parallel_pthread_create() |
---|
147 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
148 | |
---|
149 | // 2D arrays of input arguments for the <work> threads |
---|
150 | // These arrays are initialised by the application main thread |
---|
151 | |
---|
152 | work_args_t work_args[CLUSTERS_MAX][CORES_MAX]; // work function arguments |
---|
153 | work_args_t * work_ptrs[CLUSTERS_MAX][CORES_MAX]; // pointers on arguments |
---|
154 | |
---|
155 | // 1D array of barriers to allow the <work> threads to signal termination |
---|
156 | // this array is initialised in each cluster by the <build[cxy][0]> thread |
---|
157 | |
---|
158 | pthread_barrier_t parent_barriers[CLUSTERS_MAX]; // termination barrier |
---|
159 | |
---|
160 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
161 | // functions declaration |
---|
162 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
163 | |
---|
164 | void work( work_args_t * args ); |
---|
165 | |
---|
166 | double CheckSum( void ); |
---|
167 | |
---|
168 | void InitD( double ** data , |
---|
169 | unsigned int mode, |
---|
170 | unsigned int tid ); |
---|
171 | |
---|
172 | void InitT( double ** twid, |
---|
173 | unsigned int tid ); |
---|
174 | |
---|
175 | void InitU( double * coefs ); |
---|
176 | |
---|
177 | unsigned int BitReverse( unsigned int k ); |
---|
178 | |
---|
179 | void FFT1D( int direction, |
---|
180 | double ** x, |
---|
181 | double ** tmp, |
---|
182 | double * upriv, |
---|
183 | double ** twid, |
---|
184 | unsigned int tid, |
---|
185 | unsigned int MyFirst, |
---|
186 | unsigned int MyLast ); |
---|
187 | |
---|
188 | void TwiddleOneCol( int direction, |
---|
189 | unsigned int j, |
---|
190 | double ** u, |
---|
191 | double ** x, |
---|
192 | unsigned int offset_x ); |
---|
193 | |
---|
194 | void Scale( double ** x, |
---|
195 | unsigned int offset_x ); |
---|
196 | |
---|
197 | void Transpose( double ** src, |
---|
198 | double ** dest, |
---|
199 | unsigned int MyFirst, |
---|
200 | unsigned int MyLast ); |
---|
201 | |
---|
202 | void Copy( double ** src, |
---|
203 | double ** dest, |
---|
204 | unsigned int MyFirst, |
---|
205 | unsigned int MyLast ); |
---|
206 | |
---|
207 | void Reverse( double ** x, |
---|
208 | unsigned int offset_x ); |
---|
209 | |
---|
210 | void FFTRow( int direction, |
---|
211 | double * u, |
---|
212 | double ** x, |
---|
213 | unsigned int offset_x ); |
---|
214 | |
---|
215 | void PrintArray( double ** x, |
---|
216 | unsigned int size ); |
---|
217 | |
---|
218 | void SimpleDft( int direction, |
---|
219 | unsigned int size, |
---|
220 | double ** src, |
---|
221 | unsigned int src_offset, |
---|
222 | double ** dst, |
---|
223 | unsigned int dst_offset ); |
---|
224 | |
---|
225 | /////////////////////////////////////////////////////////////////// |
---|
226 | // This main() function execute the sequencial initialisation |
---|
227 | // launch the parallel execution, and makes the instrumentation. |
---|
228 | /////////////////////////////////////////////////////////////////// |
---|
229 | void main ( void ) |
---|
230 | { |
---|
231 | int error; |
---|
232 | |
---|
233 | unsigned int x_size; // number of clusters per row |
---|
234 | unsigned int y_size; // number of clusters per column |
---|
235 | unsigned int ncores; // max number of cores per cluster |
---|
236 | |
---|
237 | |
---|
238 | unsigned int x; // current index for cluster X coordinate |
---|
239 | unsigned int y; // current index for cluster Y coordinate |
---|
240 | unsigned int lid; // current index for core in a cluster |
---|
241 | unsigned int tid; // continuous thread index |
---|
242 | unsigned int cid; // cluster continuous index |
---|
243 | unsigned int cxy; // hardware specific cluster identifier |
---|
244 | |
---|
245 | char name[64]; // instrumentation file name |
---|
246 | char path[128]; // instrumentation path name |
---|
247 | char string[256]; |
---|
248 | int ret; |
---|
249 | |
---|
250 | unsigned long long start_init_cycle; |
---|
251 | unsigned long long end_init_cycle; |
---|
252 | |
---|
253 | #if DEBUG_MAIN |
---|
254 | unsigned long long debug_cycle; |
---|
255 | #endif |
---|
256 | |
---|
257 | #if CHECK |
---|
258 | double ck1; // for input/output checking |
---|
259 | double ck3; // for input/output checking |
---|
260 | #endif |
---|
261 | |
---|
262 | int pid = getpid(); |
---|
263 | |
---|
264 | // get FFT application start cycle |
---|
265 | get_cycle( &start_init_cycle ); |
---|
266 | |
---|
267 | // get platform parameters |
---|
268 | if( get_config( &x_size , &y_size , &ncores ) ) |
---|
269 | { |
---|
270 | printf("\n[fft error] cannot get hardware configuration\n"); |
---|
271 | exit( 0 ); |
---|
272 | } |
---|
273 | |
---|
274 | // check ncores |
---|
275 | if( (ncores != 1) && (ncores != 2) && (ncores != 4) ) |
---|
276 | { |
---|
277 | printf("\n[fft error] number of cores per cluster must be 1/2/4\n"); |
---|
278 | exit( 0 ); |
---|
279 | } |
---|
280 | |
---|
281 | // check x_size |
---|
282 | if( (x_size != 1) && (x_size != 2) && (x_size != 4) && (x_size != 8) && (x_size != 16) ) |
---|
283 | { |
---|
284 | printf("\n[fft error] x_size must be 1/2/4/8/16\n"); |
---|
285 | exit( 0 ); |
---|
286 | } |
---|
287 | |
---|
288 | // check y_size |
---|
289 | if( (y_size != 1) && (y_size != 2) && (y_size != 4) && (y_size != 8) && (y_size != 16) ) |
---|
290 | { |
---|
291 | printf("\n[fft error] y_size must be 1/2/4/8/16\n"); |
---|
292 | exit( 0 ); |
---|
293 | } |
---|
294 | |
---|
295 | // compute nthreads and nclusters |
---|
296 | nthreads = x_size * y_size * ncores; |
---|
297 | nclusters = x_size * y_size; |
---|
298 | |
---|
299 | // compute covering DQT size an level |
---|
300 | unsigned int z = (x_size > y_size) ? x_size : y_size; |
---|
301 | unsigned int root_level = (z == 1) ? 0 : (z == 2) ? 1 : (z == 4) ? 2 : (z == 8) ? 3 : 4; |
---|
302 | |
---|
303 | // compute various constants depending on N and T |
---|
304 | N = 1 << M; |
---|
305 | rootN = 1 << (M / 2); |
---|
306 | rows_per_thread = rootN / nthreads; |
---|
307 | points_per_cluster = N / nclusters; |
---|
308 | |
---|
309 | // check N versus T |
---|
310 | if( rootN < nthreads ) |
---|
311 | { |
---|
312 | printf("\n[fft error] sqrt(N) must be larger than T\n"); |
---|
313 | exit( 0 ); |
---|
314 | } |
---|
315 | |
---|
316 | printf("\n[fft] starts / %d points / %d thread(s) / PID %x / cycle %d\n", |
---|
317 | N, nthreads, pid, (unsigned int)start_init_cycle ); |
---|
318 | |
---|
319 | // build instrumentation file name |
---|
320 | if( USE_DQT_BARRIER ) |
---|
321 | snprintf( name , 64 , "p_fft_dqt_%d_%d_%d", N , x_size * y_size , ncores ); |
---|
322 | else |
---|
323 | snprintf( name , 64 , "p_fft_smp_%d_%d_%d", N , x_size * y_size , ncores ); |
---|
324 | |
---|
325 | // build pathname |
---|
326 | snprintf( path , 128 , "/home/%s", name ); |
---|
327 | |
---|
328 | // open instrumentation file |
---|
329 | FILE * f = fopen( path , NULL ); |
---|
330 | if ( f == NULL ) |
---|
331 | { |
---|
332 | printf("\n[fft error] cannot open instrumentation file <%s>\n", path ); |
---|
333 | exit( 0 ); |
---|
334 | } |
---|
335 | |
---|
336 | #if DEBUG_MAIN |
---|
337 | get_cycle( &debug_cycle ); |
---|
338 | printf("\n[fft] main open file <%s> at cycle %d\n", |
---|
339 | path, (unsigned int)debug_cycle ); |
---|
340 | #endif |
---|
341 | |
---|
342 | #if CHECK |
---|
343 | ck1 = CheckSum(); |
---|
344 | #endif |
---|
345 | |
---|
346 | #if PRINT_ARRAY |
---|
347 | printf("\nData values / base = %x\n", &data[0][0] ); |
---|
348 | PrintArray( data , N ); |
---|
349 | |
---|
350 | printf("\nTwiddle values / base = %x\n", &twid[0][0] ); |
---|
351 | PrintArray( twid , N ); |
---|
352 | |
---|
353 | SimpleDft( 1 , N , data , 0 , bloup , 0 ); |
---|
354 | |
---|
355 | printf("\nExpected results / base = %x\n", &bloup[0][0] ); |
---|
356 | PrintArray( bloup , N ); |
---|
357 | #endif |
---|
358 | |
---|
359 | // initialise barrier synchronizing all <work> threads |
---|
360 | if( USE_DQT_BARRIER ) |
---|
361 | { |
---|
362 | barrier_attr.x_size = x_size; |
---|
363 | barrier_attr.y_size = y_size; |
---|
364 | barrier_attr.nthreads = ncores; |
---|
365 | error = pthread_barrier_init( &barrier, &barrier_attr , nthreads ); |
---|
366 | } |
---|
367 | else |
---|
368 | { |
---|
369 | error = pthread_barrier_init( &barrier, NULL , nthreads ); |
---|
370 | } |
---|
371 | |
---|
372 | if( error ) |
---|
373 | { |
---|
374 | printf("\n[fft error] cannot initialize barrier\n"); |
---|
375 | exit( 0 ); |
---|
376 | } |
---|
377 | |
---|
378 | #if DEBUG_MAIN |
---|
379 | get_cycle( &debug_cycle ); |
---|
380 | printf("\n[fft] main completes barrier init at cycle %d\n", |
---|
381 | (unsigned int)debug_cycle ); |
---|
382 | #endif |
---|
383 | |
---|
384 | // build array of arguments for the <work> threads |
---|
385 | for (x = 0 ; x < x_size ; x++) |
---|
386 | { |
---|
387 | for (y = 0 ; y < y_size ; y++) |
---|
388 | { |
---|
389 | // compute cluster identifier |
---|
390 | cxy = HAL_CXY_FROM_XY( x , y ); |
---|
391 | |
---|
392 | for ( lid = 0 ; lid < ncores ; lid++ ) |
---|
393 | { |
---|
394 | // compute cluster continuous index |
---|
395 | cid = (x * y_size) + y; |
---|
396 | |
---|
397 | // compute work thread continuous index |
---|
398 | tid = (cid * ncores) + lid; |
---|
399 | |
---|
400 | // initialize 2D array of arguments |
---|
401 | work_args[cxy][lid].tid = tid; |
---|
402 | work_args[cxy][lid].lid = lid; |
---|
403 | work_args[cxy][lid].cid = cid; |
---|
404 | work_args[cxy][lid].parent_barrier = &parent_barriers[cxy]; |
---|
405 | |
---|
406 | // initialize 2D array of pointers |
---|
407 | work_ptrs[cxy][lid] = &work_args[cxy][lid]; |
---|
408 | } |
---|
409 | } |
---|
410 | } |
---|
411 | |
---|
412 | // register sequencial time |
---|
413 | get_cycle( &end_init_cycle ); |
---|
414 | init_time = (unsigned int)(end_init_cycle - start_init_cycle); |
---|
415 | |
---|
416 | #if DEBUG_MAIN |
---|
417 | printf("\n[fft] main completes <work> threads arguments at cycle %d\n", |
---|
418 | (unsigned int)end_init_cycle ); |
---|
419 | #endif |
---|
420 | |
---|
421 | // create and execute the working threads |
---|
422 | if( pthread_parallel_create( root_level, |
---|
423 | &work, |
---|
424 | &work_ptrs[0][0], |
---|
425 | &parent_barriers[0] ) ) |
---|
426 | { |
---|
427 | printf("\n[fft error] creating threads\n"); |
---|
428 | exit( 0 ); |
---|
429 | } |
---|
430 | |
---|
431 | #if DEBUG_MAIN |
---|
432 | get_cycle( &debug_cycle ); |
---|
433 | printf("\n[fft] main resume for instrumentation at cycle %d\n", |
---|
434 | (unsigned int)debug_cycle) ; |
---|
435 | #endif |
---|
436 | |
---|
437 | #if PRINT_ARRAY |
---|
438 | printf("\nData values after FFT:\n"); |
---|
439 | PrintArray( data , N ); |
---|
440 | #endif |
---|
441 | |
---|
442 | #if CHECK |
---|
443 | ck3 = CheckSum(); |
---|
444 | printf("\n*** Results ***\n"); |
---|
445 | printf("Checksum difference is %f (%f, %f)\n", ck1 - ck3, ck1, ck3); |
---|
446 | if (fabs(ck1 - ck3) < 0.001) printf("Results OK\n"); |
---|
447 | else printf("Results KO\n"); |
---|
448 | #endif |
---|
449 | |
---|
450 | // display header on terminal, and save to file |
---|
451 | printf("\n----- %s -----\n", name ); |
---|
452 | |
---|
453 | ret = fprintf( f , "\n----- %s -----\n", name ); |
---|
454 | if( ret < 0 ) |
---|
455 | { |
---|
456 | printf("\n[fft error] cannot write header to file <%s>\n", path ); |
---|
457 | exit(0); |
---|
458 | } |
---|
459 | |
---|
460 | // get instrumentation results for each thread |
---|
461 | for (tid = 0 ; tid < nthreads ; tid++) |
---|
462 | { |
---|
463 | snprintf( string , 256 , "- tid %d : Sequencial %d / Parallel %d / Barrier %d\n", |
---|
464 | tid, init_time, parallel_time[tid], sync_time[tid] ); |
---|
465 | |
---|
466 | // save to instrumentation file |
---|
467 | fprintf( f , "%s" , string ); |
---|
468 | if( ret < 0 ) |
---|
469 | { |
---|
470 | printf("\n[fft error] cannot write thread %d to file <%s>\n", tid, path ); |
---|
471 | printf("%s", string ); |
---|
472 | exit(0); |
---|
473 | } |
---|
474 | } |
---|
475 | |
---|
476 | // compute min/max values |
---|
477 | unsigned int min_para = parallel_time[0]; |
---|
478 | unsigned int max_para = parallel_time[0]; |
---|
479 | unsigned int min_sync = sync_time[0]; |
---|
480 | unsigned int max_sync = sync_time[0]; |
---|
481 | |
---|
482 | for (tid = 0 ; tid < nthreads ; tid++) |
---|
483 | { |
---|
484 | if (parallel_time[tid] > max_para) max_para = parallel_time[tid]; |
---|
485 | if (parallel_time[tid] < min_para) min_para = parallel_time[tid]; |
---|
486 | if (sync_time[tid] > max_sync) max_sync = sync_time[tid]; |
---|
487 | if (sync_time[tid] < min_sync) min_sync = sync_time[tid]; |
---|
488 | } |
---|
489 | |
---|
490 | // display MIN/MAX values on terminal and save to file |
---|
491 | snprintf( string , 256 , "\n Sequencial Parallel Barrier\n" |
---|
492 | "MIN : %d\t | %d\t | %d\t (cycles)\n" |
---|
493 | "MAX : %d\t | %d\t | %d\t (cycles)\n", |
---|
494 | (int)init_time, (int)min_para, (int)min_sync, |
---|
495 | (int)init_time, (int)max_para, (int)max_sync ); |
---|
496 | printf("%s", string ); |
---|
497 | ret = fprintf( f , "%s", string ); |
---|
498 | if( ret < 0 ) |
---|
499 | { |
---|
500 | printf("\n[fft error] cannot write MIN/MAX to file <%s>\n", path ); |
---|
501 | exit(0); |
---|
502 | } |
---|
503 | |
---|
504 | #if DEBUG_MAIN |
---|
505 | get_cycle( &debug_cycle ); |
---|
506 | printf("\n[fft] main close file <%s> at cycle %d\n", |
---|
507 | path, (unsigned int)debug_cycle ); |
---|
508 | #endif |
---|
509 | |
---|
510 | exit( 0 ); |
---|
511 | |
---|
512 | } // end main() |
---|
513 | |
---|
514 | ///////////////////////////////////////////////////////////////// |
---|
515 | // This function is executed in parallel by all <work> threads. |
---|
516 | ///////////////////////////////////////////////////////////////// |
---|
517 | void work( work_args_t * args ) |
---|
518 | { |
---|
519 | unsigned int tid; // this thread continuous index |
---|
520 | unsigned int lid; // core local index |
---|
521 | unsigned int cid; // cluster continuous index |
---|
522 | pthread_barrier_t * parent_barrier; // pointer on parent barrier |
---|
523 | |
---|
524 | unsigned int MyFirst; // index first row allocated to thread |
---|
525 | unsigned int MyLast; // index last row allocated to thread |
---|
526 | double * upriv; // private array of FFT coefs |
---|
527 | |
---|
528 | unsigned long long parallel_start; |
---|
529 | unsigned long long parallel_stop; |
---|
530 | unsigned long long barrier_start; |
---|
531 | unsigned long long barrier_stop; |
---|
532 | |
---|
533 | // get thread arguments |
---|
534 | tid = args->tid; |
---|
535 | lid = args->lid; |
---|
536 | cid = args->cid; |
---|
537 | parent_barrier = args->parent_barrier; |
---|
538 | |
---|
539 | get_cycle( ¶llel_start ); |
---|
540 | |
---|
541 | #if DEBUG_WORK |
---|
542 | printf("\n[fft] %s : thread %d enter / cycle %d\n", |
---|
543 | __FUNCTION__, tid, (unsigned int)parallel_start ); |
---|
544 | #endif |
---|
545 | |
---|
546 | // thread on core 0 allocates memory from the local cluster |
---|
547 | // for the distributed data[], trans[], twid[] buffers |
---|
548 | if( lid == 0 ) |
---|
549 | { |
---|
550 | unsigned int data_size = (N / nclusters) * 2 * sizeof(double); |
---|
551 | |
---|
552 | data[cid] = (double *)malloc( data_size ); |
---|
553 | if( data[cid] == NULL ) |
---|
554 | { |
---|
555 | printf("\n[fft_error] in work : cannot allocate data[%d] buffer\n", cid ); |
---|
556 | pthread_barrier_wait( parent_barrier ); |
---|
557 | pthread_exit( NULL ); |
---|
558 | } |
---|
559 | |
---|
560 | trans[cid] = (double *)malloc( data_size ); |
---|
561 | if( trans[cid] == NULL ) |
---|
562 | { |
---|
563 | printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid ); |
---|
564 | pthread_barrier_wait( parent_barrier ); |
---|
565 | pthread_exit( NULL ); |
---|
566 | } |
---|
567 | |
---|
568 | twid[cid] = (double *)malloc( data_size ); |
---|
569 | if( twid[cid] == NULL ) |
---|
570 | { |
---|
571 | printf("\n[fft_error] in work : cannot allocate twid[%d] buffer\n", cid ); |
---|
572 | pthread_barrier_wait( parent_barrier ); |
---|
573 | pthread_exit( NULL ); |
---|
574 | } |
---|
575 | } |
---|
576 | |
---|
577 | // BARRIER to wait distributed buffers allocation |
---|
578 | get_cycle( &barrier_start ); |
---|
579 | pthread_barrier_wait( &barrier ); |
---|
580 | get_cycle( &barrier_stop ); |
---|
581 | sync_time[tid] += (unsigned int)(barrier_stop - barrier_start); |
---|
582 | |
---|
583 | #if DEBUG_WORK |
---|
584 | printf("\n[fft] %s : thread %d exit barrier for buffer allocation / cycle %d\n", |
---|
585 | __FUNCTION__, tid, (unsigned int)barrier_stop ); |
---|
586 | #endif |
---|
587 | |
---|
588 | #if DISPLAY_SCHED_AND_VMM |
---|
589 | unsigned int x_size; |
---|
590 | unsigned int y_size; |
---|
591 | unsigned int ncores; |
---|
592 | get_config( &x_size , &y_size , &ncores ); |
---|
593 | unsigned int x = cid / y_size; |
---|
594 | unsigned int y = cid % y_size; |
---|
595 | unsigned int cxy = HAL_CXY_FROM_XY( x , y ); |
---|
596 | display_sched( cxy , lid ); |
---|
597 | if( lid == 0 ) display_vmm( cxy , getpid() , 0 ); |
---|
598 | #endif |
---|
599 | |
---|
600 | // all threads contribute to data[] local array initialisation |
---|
601 | InitD( data , MODE , tid ); |
---|
602 | |
---|
603 | // all threads contribute to data[] local array initialisation |
---|
604 | InitT( twid , tid ); |
---|
605 | |
---|
606 | // BARRIER to wait distributed buffers initialisation |
---|
607 | get_cycle( &barrier_start ); |
---|
608 | pthread_barrier_wait( &barrier ); |
---|
609 | get_cycle( &barrier_stop ); |
---|
610 | sync_time[tid] += (unsigned int)(barrier_stop - barrier_start); |
---|
611 | |
---|
612 | #if DEBUG_WORK |
---|
613 | printf("\n[fft] %s : thread %d exit barrier for buffer initialisation / cycle %d\n", |
---|
614 | __FUNCTION__, tid, (unsigned int)barrier_stop ); |
---|
615 | #endif |
---|
616 | |
---|
617 | // all threads allocate memory from the local cluster |
---|
618 | // for the private upriv[] buffer |
---|
619 | upriv = (double *)malloc( (rootN - 1) * 2 * sizeof(double) ); |
---|
620 | if( upriv == NULL ) |
---|
621 | { |
---|
622 | printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid ); |
---|
623 | pthread_barrier_wait( parent_barrier ); |
---|
624 | pthread_exit( NULL ); |
---|
625 | } |
---|
626 | |
---|
627 | // all threads initialise the private upriv[] array |
---|
628 | InitU( upriv ); |
---|
629 | |
---|
630 | // all threads compute first and last rows handled by the thread |
---|
631 | MyFirst = rootN * tid / nthreads; |
---|
632 | MyLast = rootN * (tid + 1) / nthreads; |
---|
633 | |
---|
634 | // all threads perform forward FFT |
---|
635 | FFT1D( 1 , data , trans , upriv , twid , tid , MyFirst , MyLast ); |
---|
636 | |
---|
637 | #if CHECK |
---|
638 | get_cycle( &barrier_start ); |
---|
639 | pthread_barrier_wait( &barrier ); |
---|
640 | get_cycle( &barrier_stop ); |
---|
641 | sync_time[tid] += (unsigned int)(barrier_stop - barrier_start); |
---|
642 | FFT1D( -1 , data , trans , upriv , twid , tid , MyFirst , MyLast ); |
---|
643 | #endif |
---|
644 | |
---|
645 | get_cycle( ¶llel_stop ); |
---|
646 | |
---|
647 | // register parallel time |
---|
648 | parallel_time[tid] = (unsigned int)(parallel_stop - parallel_start); |
---|
649 | |
---|
650 | #if DEBUG_WORK |
---|
651 | printf("\n[fft] %s : thread %d completes fft / p_start %d / p_stop %d\n", |
---|
652 | __FUNCTION__, tid, (unsigned int)parallel_start, (unsigned int)parallel_stop ); |
---|
653 | #endif |
---|
654 | |
---|
655 | // work thread signals completion to main |
---|
656 | pthread_barrier_wait( parent_barrier ); |
---|
657 | |
---|
658 | #if DEBUG_WORK |
---|
659 | printf("\n[fft] %s : thread %d exit\n", |
---|
660 | __FUNCTION__, tid ); |
---|
661 | #endif |
---|
662 | |
---|
663 | #if DISPLAY_SCHED_AND_VMM |
---|
664 | printf("\n[fft] %s : thread %d exit\n", __FUNCTION__, tid ); |
---|
665 | if( lid == 0 ) display_vmm( cxy , getpid() , 0 ); |
---|
666 | #endif |
---|
667 | |
---|
668 | // work thread exit |
---|
669 | pthread_exit( NULL ); |
---|
670 | |
---|
671 | } // end work() |
---|
672 | |
---|
673 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
674 | // This function makes the DFT from the src[nclusters][points_per_cluster] distributed |
---|
675 | // buffer, to the dst[nclusters][points_per_cluster] distributed buffer. |
---|
676 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
677 | void SimpleDft( int direction, // 1 direct / -1 reverse |
---|
678 | unsigned int size, // number of points |
---|
679 | double ** src, // source distributed buffer |
---|
680 | unsigned int src_offset, // offset in source array |
---|
681 | double ** dst, // destination distributed buffer |
---|
682 | unsigned int dst_offset ) // offset in destination array |
---|
683 | { |
---|
684 | unsigned int n , k; |
---|
685 | double phi; // 2*PI*n*k/N |
---|
686 | double u_r; // cos( phi ) |
---|
687 | double u_c; // sin( phi ) |
---|
688 | double d_r; // Re(data[n]) |
---|
689 | double d_c; // Im(data[n]) |
---|
690 | double accu_r; // Re(accu) |
---|
691 | double accu_c; // Im(accu) |
---|
692 | unsigned int c_id; // distributed buffer cluster index |
---|
693 | unsigned int c_offset; // offset in distributed buffer |
---|
694 | |
---|
695 | for ( k = 0 ; k < size ; k++ ) // loop on the output data points |
---|
696 | { |
---|
697 | // initialise accu |
---|
698 | accu_r = 0; |
---|
699 | accu_c = 0; |
---|
700 | |
---|
701 | for ( n = 0 ; n < size ; n++ ) // loop on the input data points |
---|
702 | { |
---|
703 | // compute coef |
---|
704 | phi = (double)(2*PI*n*k) / size; |
---|
705 | u_r = cos( phi ); |
---|
706 | u_c = -sin( phi ) * direction; |
---|
707 | |
---|
708 | // get input data point |
---|
709 | c_id = (src_offset + n) / (points_per_cluster); |
---|
710 | c_offset = (src_offset + n) % (points_per_cluster); |
---|
711 | d_r = src[c_id][2*c_offset]; |
---|
712 | d_c = src[c_id][2*c_offset+1]; |
---|
713 | |
---|
714 | // increment accu |
---|
715 | accu_r += ((u_r*d_r) - (u_c*d_c)); |
---|
716 | accu_c += ((u_r*d_c) + (u_c*d_r)); |
---|
717 | } |
---|
718 | |
---|
719 | // scale for inverse DFT |
---|
720 | if ( direction == -1 ) |
---|
721 | { |
---|
722 | accu_r /= size; |
---|
723 | accu_c /= size; |
---|
724 | } |
---|
725 | |
---|
726 | // set output data point |
---|
727 | c_id = (dst_offset + k) / (points_per_cluster); |
---|
728 | c_offset = (dst_offset + k) % (points_per_cluster); |
---|
729 | dst[c_id][2*c_offset] = accu_r; |
---|
730 | dst[c_id][2*c_offset+1] = accu_c; |
---|
731 | } |
---|
732 | |
---|
733 | } // end SimpleDft() |
---|
734 | |
---|
735 | /////////////////////// |
---|
736 | double CheckSum( void ) |
---|
737 | { |
---|
738 | unsigned int i , j; |
---|
739 | unsigned int c_id; |
---|
740 | unsigned int c_offset; |
---|
741 | double cks; |
---|
742 | |
---|
743 | cks = 0.0; |
---|
744 | for (j = 0; j < rootN ; j++) |
---|
745 | { |
---|
746 | for (i = 0; i < rootN ; i++) |
---|
747 | { |
---|
748 | c_id = (rootN * j + i) / (points_per_cluster); |
---|
749 | c_offset = (rootN * j + i) % (points_per_cluster); |
---|
750 | |
---|
751 | cks += data[c_id][2*c_offset] + data[c_id][2*c_offset+1]; |
---|
752 | } |
---|
753 | } |
---|
754 | return(cks); |
---|
755 | } |
---|
756 | |
---|
757 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
758 | // Each working thread <tid> contributes to initialize (rootN / nthreads) rows, |
---|
759 | // in the shared - and distributed - <data> array. |
---|
760 | ////////////////////////////////////////////////////////////////////////////////////// |
---|
761 | void InitD(double ** data, |
---|
762 | unsigned int mode, |
---|
763 | unsigned int tid ) |
---|
764 | { |
---|
765 | unsigned int i , j; |
---|
766 | unsigned int c_id; |
---|
767 | unsigned int c_offset; |
---|
768 | unsigned int index; |
---|
769 | |
---|
770 | // compute row_min and row_max |
---|
771 | unsigned int row_min = tid * rows_per_thread; |
---|
772 | unsigned int row_max = row_min + rows_per_thread; |
---|
773 | |
---|
774 | for ( j = row_min ; j < row_max ; j++ ) // loop on rows |
---|
775 | { |
---|
776 | for ( i = 0 ; i < rootN ; i++ ) // loop on points in a row |
---|
777 | { |
---|
778 | index = j * rootN + i; |
---|
779 | c_id = index / (points_per_cluster); |
---|
780 | c_offset = index % (points_per_cluster); |
---|
781 | |
---|
782 | // complex input signal is random |
---|
783 | if ( mode == RANDOM ) |
---|
784 | { |
---|
785 | data[c_id][2*c_offset] = ( (double)rand() ) / 65536; |
---|
786 | data[c_id][2*c_offset+1] = ( (double)rand() ) / 65536; |
---|
787 | } |
---|
788 | |
---|
789 | |
---|
790 | // complex input signal is cos(n/N) / sin(n/N) |
---|
791 | if ( mode == COSIN ) |
---|
792 | { |
---|
793 | double phi = (double)( 2 * PI * index) / N; |
---|
794 | data[c_id][2*c_offset] = cos( phi ); |
---|
795 | data[c_id][2*c_offset+1] = sin( phi ); |
---|
796 | } |
---|
797 | |
---|
798 | // complex input signal is constant |
---|
799 | if ( mode == CONSTANT ) |
---|
800 | { |
---|
801 | data[c_id][2*c_offset] = 1.0; |
---|
802 | data[c_id][2*c_offset+1] = 0.0; |
---|
803 | } |
---|
804 | } |
---|
805 | } |
---|
806 | } |
---|
807 | |
---|
808 | /////////////////////////////////////////////////////////////////////////////////////// |
---|
809 | // Each working thread <tid> contributes to initialize (rootN / nthreads) rows, |
---|
810 | // in the shared - and distributed - <twiddle> array. |
---|
811 | /////////////////////////////////////////////////////////////////////////////////////// |
---|
812 | void InitT( double ** twid, |
---|
813 | unsigned int tid ) |
---|
814 | { |
---|
815 | unsigned int i, j; |
---|
816 | unsigned int index; |
---|
817 | unsigned int c_id; |
---|
818 | unsigned int c_offset; |
---|
819 | double phi; |
---|
820 | |
---|
821 | // compute row_min and row_max |
---|
822 | unsigned int row_min = tid * rows_per_thread; |
---|
823 | unsigned int row_max = row_min + rows_per_thread; |
---|
824 | |
---|
825 | for ( j = row_min ; j < row_max ; j++ ) // loop on rows |
---|
826 | { |
---|
827 | for ( i = 0 ; i < rootN ; i++ ) // loop on points in a row |
---|
828 | { |
---|
829 | index = j * rootN + i; |
---|
830 | c_id = index / (points_per_cluster); |
---|
831 | c_offset = index % (points_per_cluster); |
---|
832 | |
---|
833 | phi = (double)(2.0 * PI * i * j) / N; |
---|
834 | twid[c_id][2*c_offset] = cos( phi ); |
---|
835 | twid[c_id][2*c_offset+1] = -sin( phi ); |
---|
836 | } |
---|
837 | } |
---|
838 | } |
---|
839 | |
---|
840 | /////////////////////////////////////////////////////////////////////////////////////// |
---|
841 | // Each working thread initialize the private <upriv> array / (rootN - 1) entries. |
---|
842 | /////////////////////////////////////////////////////////////////////////////////////// |
---|
843 | void InitU( double * upriv ) |
---|
844 | { |
---|
845 | unsigned int q; |
---|
846 | unsigned int j; |
---|
847 | unsigned int base; |
---|
848 | unsigned int n1; |
---|
849 | double phi; |
---|
850 | |
---|
851 | for (q = 0 ; ((unsigned int)(1 << q) < N) ; q++) |
---|
852 | { |
---|
853 | n1 = 1 << q; // n1 == 2**q |
---|
854 | base = n1 - 1; |
---|
855 | for (j = 0; (j < n1) ; j++) |
---|
856 | { |
---|
857 | if (base + j > rootN - 1) return; |
---|
858 | |
---|
859 | phi = (double)(2.0 * PI * j) / (2 * n1); |
---|
860 | upriv[2*(base+j)] = cos( phi ); |
---|
861 | upriv[2*(base+j)+1] = -sin( phi ); |
---|
862 | } |
---|
863 | } |
---|
864 | } |
---|
865 | |
---|
866 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
867 | // This function returns an index value that is the bit reverse of the input value. |
---|
868 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
869 | unsigned int BitReverse( unsigned int k ) |
---|
870 | { |
---|
871 | unsigned int i; |
---|
872 | unsigned int j; |
---|
873 | unsigned int tmp; |
---|
874 | |
---|
875 | j = 0; |
---|
876 | tmp = k; |
---|
877 | for (i = 0; i < M/2 ; i++) |
---|
878 | { |
---|
879 | j = 2 * j + (tmp & 0x1); |
---|
880 | tmp = tmp >> 1; |
---|
881 | } |
---|
882 | return j; |
---|
883 | } |
---|
884 | |
---|
885 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
886 | // This function perform the in place (direct or inverse) FFT on the N data points |
---|
887 | // contained in the distributed buffers x[nclusters][points_per_cluster]. |
---|
888 | // It handles the (N) points 1D array as a (rootN*rootN) points 2D array. |
---|
889 | // 1) it transpose (rootN/nthreads ) rows from x to tmp. |
---|
890 | // 2) it make (rootN/nthreads) FFT on the tmp rows and apply the twiddle factor. |
---|
891 | // 3) it transpose (rootN/nthreads) columns from tmp to x. |
---|
892 | // 4) it make (rootN/nthreads) FFT on the x rows. |
---|
893 | // It calls the FFTRow() 2*(rootN/nthreads) times to perform the in place FFT |
---|
894 | // on the rootN points contained in a row. |
---|
895 | //////////////////////////////////////////////////////////////////////////////////////// |
---|
896 | void FFT1D( int direction, // direct 1 / inverse -1 |
---|
897 | double ** x, // input & output distributed data points array |
---|
898 | double ** tmp, // auxiliary distributed data points array |
---|
899 | double * upriv, // local array containing coefs for rootN FFT |
---|
900 | double ** twid, // distributed arrays containing N twiddle factors |
---|
901 | unsigned int tid, // thread continuous index |
---|
902 | unsigned int MyFirst, |
---|
903 | unsigned int MyLast ) |
---|
904 | { |
---|
905 | unsigned int j; |
---|
906 | unsigned long long barrier_start; |
---|
907 | unsigned long long barrier_stop; |
---|
908 | |
---|
909 | #if DEBUG_FFT1D |
---|
910 | unsigned long long cycle; |
---|
911 | get_cycle( &cycle ); |
---|
912 | printf("\n[fft] %s : thread %d enter / first %d / last %d / cycle %d\n", |
---|
913 | __FUNCTION__, tid, MyFirst, MyLast, (unsigned int)cycle ); |
---|
914 | #endif |
---|
915 | |
---|
916 | // transpose (rootN/nthreads) rows from x to tmp |
---|
917 | Transpose( x , tmp , MyFirst , MyLast ); |
---|
918 | |
---|
919 | #if( DEBUG_FFT1D & 1 ) |
---|
920 | get_cycle( &cycle ); |
---|
921 | printf("\n[fft] %s : thread %d after first transpose / cycle %d\n", |
---|
922 | __FUNCTION__, tid, (unsigned int)cycle ); |
---|
923 | if( PRINT_ARRAY ) PrintArray( tmp , N ); |
---|
924 | #endif |
---|
925 | |
---|
926 | // BARRIER |
---|
927 | get_cycle( &barrier_start ); |
---|
928 | pthread_barrier_wait( &barrier ); |
---|
929 | get_cycle( &barrier_stop ); |
---|
930 | sync_time[tid] = (unsigned int)(barrier_stop - barrier_start); |
---|
931 | |
---|
932 | #if( DEBUG_FFT1D & 1 ) |
---|
933 | get_cycle( &cycle ); |
---|
934 | printf("\n[fft] %s : thread %d exit barrier after first transpose / cycle %d\n", |
---|
935 | __FUNCTION__, tid, (unsigned int)cycle ); |
---|
936 | #endif |
---|
937 | |
---|
938 | // do FFTs on rows of tmp (i.e. columns of x) and apply twiddle factor |
---|
939 | for (j = MyFirst; j < MyLast; j++) |
---|
940 | { |
---|
941 | FFTRow( direction , upriv , tmp , j * rootN ); |
---|
942 | |
---|
943 | TwiddleOneCol( direction , j , twid , tmp , j * rootN ); |
---|
944 | } |
---|
945 | |
---|
946 | #if( DEBUG_FFT1D & 1 ) |
---|
947 | printf("\n[fft] %s : thread %d after first twiddle\n", __FUNCTION__, tid); |
---|
948 | if( PRINT_ARRAY ) PrintArray( tmp , N ); |
---|
949 | #endif |
---|
950 | |
---|
951 | // BARRIER |
---|
952 | get_cycle( &barrier_start ); |
---|
953 | pthread_barrier_wait( &barrier ); |
---|
954 | get_cycle( &barrier_stop ); |
---|
955 | |
---|
956 | #if( DEBUG_FFT1D & 1 ) |
---|
957 | printf("\n[fft] %s : thread %d exit barrier after first twiddle\n", __FUNCTION__, tid); |
---|
958 | #endif |
---|
959 | |
---|
960 | sync_time[tid] += (unsigned int)(barrier_stop - barrier_start); |
---|
961 | |
---|
962 | // transpose tmp to x |
---|
963 | Transpose( tmp , x , MyFirst , MyLast ); |
---|
964 | |
---|
965 | #if( DEBUG_FFT1D & 1 ) |
---|
966 | printf("\n[fft] %s : thread %d after second transpose\n", __FUNCTION__, tid); |
---|
967 | if( PRINT_ARRAY ) PrintArray( x , N ); |
---|
968 | #endif |
---|
969 | |
---|
970 | // BARRIER |
---|
971 | get_cycle( &barrier_start ); |
---|
972 | pthread_barrier_wait( &barrier ); |
---|
973 | get_cycle( &barrier_stop ); |
---|
974 | |
---|
975 | #if( DEBUG_FFT1D & 1 ) |
---|
976 | printf("\n[fft] %s : thread %d exit barrier after second transpose\n", __FUNCTION__, tid); |
---|
977 | #endif |
---|
978 | |
---|
979 | sync_time[tid] += (unsigned int)(barrier_stop - barrier_start); |
---|
980 | |
---|
981 | // do FFTs on rows of x and apply the scaling factor |
---|
982 | for (j = MyFirst; j < MyLast; j++) |
---|
983 | { |
---|
984 | FFTRow( direction , upriv , x , j * rootN ); |
---|
985 | if (direction == -1) Scale( x , j * rootN ); |
---|
986 | } |
---|
987 | |
---|
988 | #if( DEBUG_FFT1D & 1 ) |
---|
989 | printf("\n[fft] %s : thread %d after FFT on rows\n", __FUNCTION__, tid); |
---|
990 | if( PRINT_ARRAY ) PrintArray( x , N ); |
---|
991 | #endif |
---|
992 | |
---|
993 | // BARRIER |
---|
994 | get_cycle( &barrier_start ); |
---|
995 | pthread_barrier_wait( &barrier ); |
---|
996 | get_cycle( &barrier_stop ); |
---|
997 | |
---|
998 | #if( DEBUG_FFT1D & 1 ) |
---|
999 | printf("\n[fft] %s : thread %d exit barrier after FFT on rows\n", __FUNCTION__, tid); |
---|
1000 | #endif |
---|
1001 | sync_time[tid] += (unsigned int)(barrier_stop - barrier_start); |
---|
1002 | |
---|
1003 | // transpose x to tmp |
---|
1004 | Transpose( x , tmp , MyFirst , MyLast ); |
---|
1005 | |
---|
1006 | #if( DEBUG_FFT1D & 1 ) |
---|
1007 | printf("\n[fft] %s : thread %x after third transpose\n", __FUNCTION__, tid); |
---|
1008 | if( PRINT_ARRAY ) PrintArray( x , N ); |
---|
1009 | #endif |
---|
1010 | |
---|
1011 | // BARRIER |
---|
1012 | get_cycle( &barrier_start ); |
---|
1013 | pthread_barrier_wait( &barrier ); |
---|
1014 | get_cycle( &barrier_stop ); |
---|
1015 | |
---|
1016 | #if( DEBUG_FFT1D & 1 ) |
---|
1017 | printf("\n[fft] %s : thread %d exit barrier after third transpose\n", __FUNCTION__, tid); |
---|
1018 | #endif |
---|
1019 | |
---|
1020 | sync_time[tid] += (unsigned int)(barrier_stop - barrier_start); |
---|
1021 | sync_time[tid] += (long)(barrier_stop - barrier_start); |
---|
1022 | |
---|
1023 | // copy tmp to x |
---|
1024 | Copy( tmp , x , MyFirst , MyLast ); |
---|
1025 | |
---|
1026 | #if DEBUG_FFT1D |
---|
1027 | printf("\n[fft] %s : thread %d completed\n", __FUNCTION__, tid); |
---|
1028 | if( PRINT_ARRAY ) PrintArray( x , N ); |
---|
1029 | #endif |
---|
1030 | |
---|
1031 | } // end FFT1D() |
---|
1032 | |
---|
1033 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
1034 | // This function multiply all points contained in a row (rootN points) of the |
---|
1035 | // x[] array by the corresponding twiddle factor, contained in the u[] array. |
---|
1036 | ///////////////////////////////////////////////////////////////////////////////////// |
---|
1037 | void TwiddleOneCol( int direction, |
---|
1038 | unsigned int j, // y coordinate in 2D view of coef array |
---|
1039 | double ** u, // coef array base address |
---|
1040 | double ** x, // data array base address |
---|
1041 | unsigned int offset_x ) // first point in N points data array |
---|
1042 | { |
---|
1043 | unsigned int i; |
---|
1044 | double omega_r; |
---|
1045 | double omega_c; |
---|
1046 | double x_r; |
---|
1047 | double x_c; |
---|
1048 | unsigned int c_id; |
---|
1049 | unsigned int c_offset; |
---|
1050 | |
---|
1051 | for (i = 0; i < rootN ; i++) // loop on the rootN points |
---|
1052 | { |
---|
1053 | // get coef |
---|
1054 | c_id = (j * rootN + i) / (points_per_cluster); |
---|
1055 | c_offset = (j * rootN + i) % (points_per_cluster); |
---|
1056 | omega_r = u[c_id][2*c_offset]; |
---|
1057 | omega_c = direction * u[c_id][2*c_offset+1]; |
---|
1058 | |
---|
1059 | // access data |
---|
1060 | c_id = (offset_x + i) / (points_per_cluster); |
---|
1061 | c_offset = (offset_x + i) % (points_per_cluster); |
---|
1062 | x_r = x[c_id][2*c_offset]; |
---|
1063 | x_c = x[c_id][2*c_offset+1]; |
---|
1064 | |
---|
1065 | x[c_id][2*c_offset] = omega_r*x_r - omega_c * x_c; |
---|
1066 | x[c_id][2*c_offset+1] = omega_r*x_c + omega_c * x_r; |
---|
1067 | } |
---|
1068 | } // end TwiddleOneCol() |
---|
1069 | |
---|
1070 | //////////////////////////// |
---|
1071 | void Scale( double ** x, // data array base address |
---|
1072 | unsigned int offset_x ) // first point of the row to be scaled |
---|
1073 | { |
---|
1074 | unsigned int i; |
---|
1075 | unsigned int c_id; |
---|
1076 | unsigned int c_offset; |
---|
1077 | |
---|
1078 | for (i = 0; i < rootN ; i++) |
---|
1079 | { |
---|
1080 | c_id = (offset_x + i) / (points_per_cluster); |
---|
1081 | c_offset = (offset_x + i) % (points_per_cluster); |
---|
1082 | x[c_id][2*c_offset] /= N; |
---|
1083 | x[c_id][2*c_offset + 1] /= N; |
---|
1084 | } |
---|
1085 | } |
---|
1086 | |
---|
1087 | /////////////////////////////////// |
---|
1088 | void Transpose( double ** src, // source buffer (array of pointers) |
---|
1089 | double ** dest, // destination buffer (array of pointers) |
---|
1090 | unsigned int MyFirst, // first row allocated to the thread |
---|
1091 | unsigned int MyLast ) // last row allocated to the thread |
---|
1092 | { |
---|
1093 | unsigned int row; // row index |
---|
1094 | unsigned int point; // data point index in a row |
---|
1095 | |
---|
1096 | unsigned int index_src; // absolute index in the source N points array |
---|
1097 | unsigned int c_id_src; // cluster for the source buffer |
---|
1098 | unsigned int c_offset_src; // offset in the source buffer |
---|
1099 | |
---|
1100 | unsigned int index_dst; // absolute index in the dest N points array |
---|
1101 | unsigned int c_id_dst; // cluster for the dest buffer |
---|
1102 | unsigned int c_offset_dst; // offset in the dest buffer |
---|
1103 | |
---|
1104 | |
---|
1105 | // scan all data points allocated to the thread |
---|
1106 | // (between MyFirst row and MyLast row) from the source buffer |
---|
1107 | // and write these points to the destination buffer |
---|
1108 | for ( row = MyFirst ; row < MyLast ; row++ ) // loop on the rows |
---|
1109 | { |
---|
1110 | for ( point = 0 ; point < rootN ; point++ ) // loop on points in row |
---|
1111 | { |
---|
1112 | index_src = row * rootN + point; |
---|
1113 | c_id_src = index_src / (points_per_cluster); |
---|
1114 | c_offset_src = index_src % (points_per_cluster); |
---|
1115 | |
---|
1116 | index_dst = point * rootN + row; |
---|
1117 | c_id_dst = index_dst / (points_per_cluster); |
---|
1118 | c_offset_dst = index_dst % (points_per_cluster); |
---|
1119 | |
---|
1120 | dest[c_id_dst][2*c_offset_dst] = src[c_id_src][2*c_offset_src]; |
---|
1121 | dest[c_id_dst][2*c_offset_dst+1] = src[c_id_src][2*c_offset_src+1]; |
---|
1122 | } |
---|
1123 | } |
---|
1124 | } // end Transpose() |
---|
1125 | |
---|
1126 | ////////////////////////////// |
---|
1127 | void Copy( double ** src, // source buffer (array of pointers) |
---|
1128 | double ** dest, // destination buffer (array of pointers) |
---|
1129 | unsigned int MyFirst, // first row allocated to the thread |
---|
1130 | unsigned int MyLast ) // last row allocated to the thread |
---|
1131 | { |
---|
1132 | unsigned int row; // row index |
---|
1133 | unsigned int point; // data point index in a row |
---|
1134 | |
---|
1135 | unsigned int index; // absolute index in the N points array |
---|
1136 | unsigned int c_id; // cluster index |
---|
1137 | unsigned int c_offset; // offset in local buffer |
---|
1138 | |
---|
1139 | // scan all data points allocated to the thread |
---|
1140 | for ( row = MyFirst ; row < MyLast ; row++ ) // loop on the rows |
---|
1141 | { |
---|
1142 | for ( point = 0 ; point < rootN ; point++ ) // loop on points in row |
---|
1143 | { |
---|
1144 | index = row * rootN + point; |
---|
1145 | c_id = index / (points_per_cluster); |
---|
1146 | c_offset = index % (points_per_cluster); |
---|
1147 | |
---|
1148 | dest[c_id][2*c_offset] = src[c_id][2*c_offset]; |
---|
1149 | dest[c_id][2*c_offset+1] = src[c_id][2*c_offset+1]; |
---|
1150 | } |
---|
1151 | } |
---|
1152 | } // end Copy() |
---|
1153 | |
---|
1154 | /////////////////////////////// |
---|
1155 | void Reverse( double ** x, |
---|
1156 | unsigned int offset_x ) |
---|
1157 | { |
---|
1158 | unsigned int j, k; |
---|
1159 | unsigned int c_id_j; |
---|
1160 | unsigned int c_offset_j; |
---|
1161 | unsigned int c_id_k; |
---|
1162 | unsigned int c_offset_k; |
---|
1163 | |
---|
1164 | for (k = 0 ; k < rootN ; k++) |
---|
1165 | { |
---|
1166 | j = BitReverse( k ); |
---|
1167 | if (j > k) |
---|
1168 | { |
---|
1169 | c_id_j = (offset_x + j) / (points_per_cluster); |
---|
1170 | c_offset_j = (offset_x + j) % (points_per_cluster); |
---|
1171 | c_id_k = (offset_x + k) / (points_per_cluster); |
---|
1172 | c_offset_k = (offset_x + k) % (points_per_cluster); |
---|
1173 | |
---|
1174 | SWAP(x[c_id_j][2*c_offset_j] , x[c_id_k][2*c_offset_k]); |
---|
1175 | SWAP(x[c_id_j][2*c_offset_j+1], x[c_id_k][2*c_offset_k+1]); |
---|
1176 | } |
---|
1177 | } |
---|
1178 | } |
---|
1179 | |
---|
1180 | ///////////////////////////////////////////////////////////////////////////// |
---|
1181 | // This function makes the in-place FFT on all points contained in a row |
---|
1182 | // (i.e. rootN points) of the x[nclusters][points_per_cluster] array. |
---|
1183 | ///////////////////////////////////////////////////////////////////////////// |
---|
1184 | void FFTRow( int direction, // 1 direct / -1 inverse |
---|
1185 | double * u, // private coefs array |
---|
1186 | double ** x, // array of pointers on distributed buffers |
---|
1187 | unsigned int offset_x ) // absolute offset in the x array |
---|
1188 | { |
---|
1189 | unsigned int j; |
---|
1190 | unsigned int k; |
---|
1191 | unsigned int q; |
---|
1192 | unsigned int L; |
---|
1193 | unsigned int r; |
---|
1194 | unsigned int Lstar; |
---|
1195 | double * u1; |
---|
1196 | |
---|
1197 | unsigned int offset_x1; // index first butterfly input |
---|
1198 | unsigned int offset_x2; // index second butterfly output |
---|
1199 | |
---|
1200 | double omega_r; // real part butterfy coef |
---|
1201 | double omega_c; // complex part butterfly coef |
---|
1202 | |
---|
1203 | double tau_r; |
---|
1204 | double tau_c; |
---|
1205 | |
---|
1206 | double d1_r; // real part first butterfly input |
---|
1207 | double d1_c; // imag part first butterfly input |
---|
1208 | double d2_r; // real part second butterfly input |
---|
1209 | double d2_c; // imag part second butterfly input |
---|
1210 | |
---|
1211 | unsigned int c_id_1; // cluster index for first butterfly input |
---|
1212 | unsigned int c_offset_1; // offset for first butterfly input |
---|
1213 | unsigned int c_id_2; // cluster index for second butterfly input |
---|
1214 | unsigned int c_offset_2; // offset for second butterfly input |
---|
1215 | |
---|
1216 | #if DEBUG_ROW |
---|
1217 | unsigned int p; |
---|
1218 | printf("\n[fft] ROW data in / %d points / offset = %d\n", rootN , offset_x ); |
---|
1219 | |
---|
1220 | for ( p = 0 ; p < rootN ; p++ ) |
---|
1221 | { |
---|
1222 | unsigned int index = offset_x + p; |
---|
1223 | unsigned int c_id = index / (points_per_cluster); |
---|
1224 | unsigned int c_offset = index % (points_per_cluster); |
---|
1225 | printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] ); |
---|
1226 | } |
---|
1227 | printf("\n"); |
---|
1228 | #endif |
---|
1229 | |
---|
1230 | // This makes the rootN input points reordering |
---|
1231 | Reverse( x , offset_x ); |
---|
1232 | |
---|
1233 | #if DEBUG_ROW |
---|
1234 | printf("\n[fft] ROW data after reverse / %d points / offset = %d\n", rootN , offset_x ); |
---|
1235 | |
---|
1236 | for ( p = 0 ; p < rootN ; p++ ) |
---|
1237 | { |
---|
1238 | unsigned int index = offset_x + p; |
---|
1239 | unsigned int c_id = index / (points_per_cluster); |
---|
1240 | unsigned int c_offset = index % (points_per_cluster); |
---|
1241 | printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] ); |
---|
1242 | } |
---|
1243 | printf("\n"); |
---|
1244 | #endif |
---|
1245 | |
---|
1246 | // This implements the multi-stages, in place Butterfly network |
---|
1247 | for (q = 1; q <= M/2 ; q++) // loop on stages |
---|
1248 | { |
---|
1249 | L = 1 << q; // number of points per subset for current stage |
---|
1250 | r = rootN / L; // number of subsets |
---|
1251 | Lstar = L / 2; |
---|
1252 | u1 = &u[2 * (Lstar - 1)]; |
---|
1253 | for (k = 0; k < r; k++) // loop on the subsets |
---|
1254 | { |
---|
1255 | offset_x1 = offset_x + (k * L); // index first point |
---|
1256 | offset_x2 = offset_x + (k * L + Lstar); // index second point |
---|
1257 | |
---|
1258 | #if (DEBUG_ROW & 1) |
---|
1259 | printf("\n ### q = %d / k = %d / x1 = %d / x2 = %d\n", q , k , offset_x1 , offset_x2 ); |
---|
1260 | #endif |
---|
1261 | // makes all in-place butterfly(s) for subset |
---|
1262 | for (j = 0; j < Lstar; j++) |
---|
1263 | { |
---|
1264 | // get coef |
---|
1265 | omega_r = u1[2*j]; |
---|
1266 | omega_c = direction * u1[2*j+1]; |
---|
1267 | |
---|
1268 | // get d[x1] address and value |
---|
1269 | c_id_1 = (offset_x1 + j) / (points_per_cluster); |
---|
1270 | c_offset_1 = (offset_x1 + j) % (points_per_cluster); |
---|
1271 | d1_r = x[c_id_1][2*c_offset_1]; |
---|
1272 | d1_c = x[c_id_1][2*c_offset_1+1]; |
---|
1273 | |
---|
1274 | // get d[x2] address and value |
---|
1275 | c_id_2 = (offset_x2 + j) / (points_per_cluster); |
---|
1276 | c_offset_2 = (offset_x2 + j) % (points_per_cluster); |
---|
1277 | d2_r = x[c_id_2][2*c_offset_2]; |
---|
1278 | d2_c = x[c_id_2][2*c_offset_2+1]; |
---|
1279 | |
---|
1280 | #if (DEBUG_ROW & 1) |
---|
1281 | printf("\n ### d1_in = (%f , %f) / d2_in = (%f , %f) / coef = (%f , %f)\n", |
---|
1282 | d1_r , d1_c , d2_r , d2_c , omega_r , omega_c); |
---|
1283 | #endif |
---|
1284 | // tau = omega * d[x2] |
---|
1285 | tau_r = omega_r * d2_r - omega_c * d2_c; |
---|
1286 | tau_c = omega_r * d2_c + omega_c * d2_r; |
---|
1287 | |
---|
1288 | // set new value for d[x1] = d[x1] + omega * d[x2] |
---|
1289 | x[c_id_1][2*c_offset_1] = d1_r + tau_r; |
---|
1290 | x[c_id_1][2*c_offset_1+1] = d1_c + tau_c; |
---|
1291 | |
---|
1292 | // set new value for d[x2] = d[x1] - omega * d[x2] |
---|
1293 | x[c_id_2][2*c_offset_2] = d1_r - tau_r; |
---|
1294 | x[c_id_2][2*c_offset_2+1] = d1_c - tau_c; |
---|
1295 | |
---|
1296 | #if (DEBUG_ROW & 1) |
---|
1297 | printf("\n ### d1_out = (%f , %f) / d2_out = (%f , %f)\n", |
---|
1298 | d1_r + tau_r , d1_c + tau_c , d2_r - tau_r , d2_c - tau_c ); |
---|
1299 | #endif |
---|
1300 | } |
---|
1301 | } |
---|
1302 | } |
---|
1303 | |
---|
1304 | #if DEBUG_ROW |
---|
1305 | printf("\n[fft] ROW data out / %d points / offset = %d\n", rootN , offset_x ); |
---|
1306 | for ( p = 0 ; p < rootN ; p++ ) |
---|
1307 | { |
---|
1308 | unsigned int index = offset_x + p; |
---|
1309 | unsigned int c_id = index / (points_per_cluster); |
---|
1310 | unsigned int c_offset = index % (points_per_cluster); |
---|
1311 | printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] ); |
---|
1312 | } |
---|
1313 | printf("\n"); |
---|
1314 | #endif |
---|
1315 | |
---|
1316 | } // end FFTRow() |
---|
1317 | |
---|
1318 | /////////////////////////////////////// |
---|
1319 | void PrintArray( double ** array, |
---|
1320 | unsigned int size ) |
---|
1321 | { |
---|
1322 | unsigned int i; |
---|
1323 | unsigned int c_id; |
---|
1324 | unsigned int c_offset; |
---|
1325 | |
---|
1326 | // float display |
---|
1327 | for (i = 0; i < size ; i++) |
---|
1328 | { |
---|
1329 | c_id = i / (points_per_cluster); |
---|
1330 | c_offset = i % (points_per_cluster); |
---|
1331 | |
---|
1332 | printf(" %f %f |", array[c_id][2*c_offset], array[c_id][2*c_offset+1]); |
---|
1333 | |
---|
1334 | if ( (i+1) % 4 == 0) printf("\n"); |
---|
1335 | } |
---|
1336 | printf("\n"); |
---|
1337 | } |
---|
1338 | |
---|
1339 | |
---|
1340 | // Local Variables: |
---|
1341 | // tab-width: 4 |
---|
1342 | // c-basic-offset: 4 |
---|
1343 | // c-file-offsets:((innamespace . 0)(inline-open . 0)) |
---|
1344 | // indent-tabs-mode: nil |
---|
1345 | // End: |
---|
1346 | |
---|
1347 | // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 |
---|
1348 | |
---|