source: trunk/user/fft/fft.c

Last change on this file was 659, checked in by alain, 6 years ago

euh...
File size: 54.3 KB
Line 
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// The actual number of cores and cluster in a given hardware architecture
48// is obtained by the get_config() syscall (x_size, y_size, ncores).
49// The max number of clusters is bounded by (X_MAX * Y_MAX).
50// The max number of cores per cluster is bounded by CORES_MAX.
51//
52// The number N of working threads is always defined by the number of cores availables
53// in the architecture, but this application supports three placement modes.
54// In all modes, the working threads are identified by the [tid] continuous index
55// in range [0, NTHREADS-1], and defines how the lines are shared amongst the threads.
56// This continuous index can always be decomposed in two continuous sub-indexes:
57// tid == cid * ncores + lid, where cid is in [0,NCLUSTERS-1] and lid in [0,NCORES-1].
58//
59// - NO_PLACEMENT: the main thread is itsef a working thread. The (N_1) other working
60// threads are created by the main thread, but the placement is done by the OS, using
61// the DQDT for load balancing, and two working threads can be placed on the same core.
62// The [cid,lid] are only abstract identifiers, and cannot be associated to a physical
63// cluster or a physical core. In this mode, the main thread run on any cluster,
64// but has tid = 0 (i.e. cid = 0 & tid = 0).
65//
66// - EXPLICIT_PLACEMENT: the main thread is again a working thread, but the placement of
67// of the threads on the cores is explicitely controled by the main thread to have
68// exactly one working thread per core, and the [cxy][lpid] core coordinates for a given
69// thread[tid] can be directly derived from the [tid] value: [cid] is an alias for the
70// physical cluster identifier, and [lid] is the local core index.
71//
72// - PARALLEL_PLACEMENT: the main thread is not anymore a working thread, and uses the
73// non standard pthread_parallel_create() function to avoid the costly sequencial
74// loops for pthread_create() and pthread_join(). It garanty one working thread
75// per core, and the same relation between the thread[tid] and the core[cxy][lpid].
76//
77// Several others configuration parameters can be defined below:
78// - USE_DQT_BARRIER : use a hierarchical barrier for working threads synchro
79// - PRINT_ARRAY : Print out complex data points arrays.
80// - CHECK : Perform both FFT and inverse FFT to check output/input.
81// - DEBUG_MAIN : Display intermediate results in main()
82// - DEBUG_FFT1D : Display intermediate results in FFT1D()
83// - DEBUG_ROW : Display intermedite results in FFTrow()
84//
85// Regarding final instrumentation:
86// - the sequencial initialisation time (init_time) is computed
87// by the main thread in the main() function.
88// - The parallel execution time (parallel_time[i]) is computed by each
89// working thread(i) in the work() function.
90// - The synchronisation time related to the barriers (sync_time[i])
91// is computed by each thread(i) in the work() function.
92// The results are displayed on the TXT terminal, and registered on disk.
93///////////////////////////////////////////////////////////////////////////////////////
94
95#include <math.h>
96#include <stdio.h>
97#include <stdlib.h>
98#include <fcntl.h>
99#include <unistd.h>
100#include <pthread.h>
101#include <almosmkh.h>
102#include <hal_macros.h>
103
104// constants
105
106#define PI 3.14159265359
107#define PAGE_SIZE 4096
108#define X_MAX 16 // max number of clusters in a row
109#define Y_MAX 16 // max number of clusters in a column
110#define CORES_MAX 4 // max number of cores in a cluster
111#define CLUSTERS_MAX X_MAX * Y_MAX
112#define THREADS_MAX CLUSTERS_MAX * CORES_MAX
113#define RANDOM 0
114#define COSIN 1
115#define CONSTANT 2
116
117// parameters
118
119#define NO_PLACEMENT 1
120#define EXPLICIT_PLACEMENT 0
121#define PARALLEL_PLACEMENT 0
122
123#define DEFAULT_M 18 // 256 K complex points
124#define USE_DQT_BARRIER 1 // use DDT barrier if non zero
125#define MODE COSIN // DATA array initialisation mode
126#define CHECK 0
127#define DEBUG_MAIN 1 // trace main() function (detailed if odd)
128#define DEBUG_WORK 0 // trace work() function (detailed if odd)
129#define DEBUG_FFT1D 0 // trace FFT1D() function (detailed if odd)
130#define DEBUG_ROW 0 // trace FFTRow() function (detailed if odd)
131#define PRINT_ARRAY 0
132#define DISPLAY_SCHED_AND_VMM 0 // display final VMM state in all clusters
133
134// macro to swap two variables
135#define SWAP(a,b) { double tmp; tmp = a; a = b; b = tmp; }
136
137/////////////////////////////////////////////////////////////////////////////////////
138// FFT global variables
139/////////////////////////////////////////////////////////////////////////////////////
140
141unsigned int x_size; // platform global parameter
142unsigned int y_size; // platform global parameter
143unsigned int ncores; // platform global parameter
144
145unsigned int nthreads; // total number of threads (one thread per core)
146unsigned int nclusters; // total number of clusters
147unsigned int M = DEFAULT_M; // log2(number of points)
148unsigned int N; // number of points (N = 2^M)
149unsigned int rootN; // rootN = 2^M/2
150unsigned int rows_per_thread; // number of data "rows" handled by a single thread
151unsigned int points_per_cluster; // number of data points per cluster
152
153// arrays of pointers on distributed buffers (one sub-buffer per cluster)
154double * data[CLUSTERS_MAX]; // original time-domain data
155double * trans[CLUSTERS_MAX]; // used as auxiliary space for fft
156double * twid[CLUSTERS_MAX]; // twiddle factor : exp(-2iPI*k*n/N)
157double * bloup[CLUSTERS_MAX]; // used as auxiliary space for DFT
158
159// instrumentation counters
160unsigned int pgfault_nr[THREADS_MAX]; // total number of page faults (per thread)
161unsigned int pgfault_cost[THREADS_MAX]; // total page faults cost (per thread)
162unsigned int pgfault_max[THREADS_MAX]; // max page faults cost (per thread)
163unsigned int parallel_time[THREADS_MAX]; // total computation time (per thread)
164unsigned int sync_time[THREADS_MAX]; // cumulated waiting time in barriers (per thread)
165unsigned int init_time; // initialisation time (in main)
166
167// synchronisation barrier (all threads)
168pthread_barrier_t barrier;
169pthread_barrierattr_t barrier_attr;
170
171//return values at thread exit
172unsigned int THREAD_EXIT_SUCCESS = 0;
173unsigned int THREAD_EXIT_FAILURE = 1;
174
175// main thread continuous index
176unsigned int tid_main;
177
178// array of kernel thread identifiers / indexed by [tid]
179pthread_t work_trdid[CLUSTERS_MAX * CORES_MAX];
180
181// array of thread attributes / indexed by [tid]
182pthread_attr_t work_attr[CLUSTERS_MAX * CORES_MAX];
183
184// array of work function arguments / indexed by [tid]
185pthread_parallel_work_args_t work_args[CLUSTERS_MAX * CORES_MAX];
186
187/////////////////////////////////////////////////////////////////////////////////////
188// functions declaration
189/////////////////////////////////////////////////////////////////////////////////////
190
191void * work( void * args );
192
193double CheckSum( void );
194
195void InitD( double ** data ,
196 unsigned int mode,
197 unsigned int tid );
198
199void InitT( double ** twid,
200 unsigned int tid );
201
202void InitU( double * coefs );
203
204unsigned int BitReverse( unsigned int k );
205
206void FFT1D( int direction,
207 double ** x,
208 double ** tmp,
209 double * upriv,
210 double ** twid,
211 unsigned int tid,
212 unsigned int MyFirst,
213 unsigned int MyLast );
214
215void TwiddleOneCol( int direction,
216 unsigned int j,
217 double ** u,
218 double ** x,
219 unsigned int offset_x );
220
221void Scale( double ** x,
222 unsigned int offset_x );
223
224void Transpose( double ** src,
225 double ** dest,
226 unsigned int MyFirst,
227 unsigned int MyLast );
228
229void Copy( double ** src,
230 double ** dest,
231 unsigned int MyFirst,
232 unsigned int MyLast );
233
234void Reverse( double ** x,
235 unsigned int offset_x );
236
237void FFTRow( int direction,
238 double * u,
239 double ** x,
240 unsigned int offset_x );
241
242void PrintArray( double ** x,
243 unsigned int size );
244
245void SimpleDft( int direction,
246 unsigned int size,
247 double ** src,
248 unsigned int src_offset,
249 double ** dst,
250 unsigned int dst_offset );
251
252///////////////////////////////////////////////////////////////////
253// This main() function execute the sequencial initialisation
254// launch the parallel execution, and makes the instrumentation.
255///////////////////////////////////////////////////////////////////
256int main ( void )
257{
258 int error;
259
260 unsigned int tid; // continuous thread index
261
262 char name[64]; // instrumentation file name
263 char path[128]; // instrumentation path name
264 char string[256];
265 int ret;
266
267 unsigned long long start_init_cycle;
268 unsigned long long end_init_cycle;
269
270#if DEBUG_MAIN
271 unsigned long long debug_cycle;
272#endif
273
274#if CHECK
275 double ck1; // for input/output checking
276 double ck3; // for input/output checking
277#endif
278
279 int pid = getpid();
280
281 // check placement mode
282 if( (NO_PLACEMENT + EXPLICIT_PLACEMENT + PARALLEL_PLACEMENT) != 1 )
283 {
284 printf("\n[fft error] illegal placement mode\n");
285 exit( 0 );
286 }
287
288 // get FFT application start cycle
289 get_cycle( &start_init_cycle );
290
291 // get platform parameters
292 hard_config_t config;
293 if( get_config( &config ) )
294 {
295 printf("\n[fft error] cannot get hardware configuration\n");
296 exit( 0 );
297 }
298
299 x_size = config.x_size;
300 y_size = config.y_size;
301 ncores = config.ncores;
302
303 // check ncores
304 if( (ncores != 1) && (ncores != 2) && (ncores != 4) )
305 {
306 printf("\n[fft error] number of cores per cluster must be 1/2/4\n");
307 exit( 0 );
308 }
309
310 // check x_size
311 if( (x_size != 1) && (x_size != 2) && (x_size != 4) && (x_size != 8) && (x_size != 16) )
312 {
313 printf("\n[fft error] x_size must be 1/2/4/8/16\n");
314 exit( 0 );
315 }
316
317 // check y_size
318 if( (y_size != 1) && (y_size != 2) && (y_size != 4) && (y_size != 8) && (y_size != 16) )
319 {
320 printf("\n[fft error] y_size must be 1/2/4/8/16\n");
321 exit( 0 );
322 }
323
324 // get identifiers for core executing main
325 unsigned int cxy_main;
326 unsigned int lid_main;
327 get_core_id( &cxy_main , &lid_main );
328
329 // compute nthreads and nclusters
330 nthreads = x_size * y_size * ncores;
331 nclusters = x_size * y_size;
332
333 // compute covering DQT size an level
334 unsigned int z = (x_size > y_size) ? x_size : y_size;
335 unsigned int root_level = (z == 1) ? 0 : (z == 2) ? 1 : (z == 4) ? 2 : (z == 8) ? 3 : 4;
336
337 // compute various constants depending on N and T
338 N = 1 << M;
339 rootN = 1 << (M / 2);
340 rows_per_thread = rootN / nthreads;
341 points_per_cluster = N / nclusters;
342
343 // check N versus T
344 if( rootN < nthreads )
345 {
346 printf("\n[fft error] sqrt(N) must be larger than T\n");
347 exit( 0 );
348 }
349
350 // define instrumentation file name
351 if( NO_PLACEMENT )
352 {
353 printf("\n[fft] starts / %d points / %d thread(s) / PID %x / NO_PLACE\n",
354 N, nthreads, pid );
355
356 // build instrumentation file name
357 if( USE_DQT_BARRIER )
358 snprintf( name , 64 , "fft_dqt_no_place_%d_%d_%d", M , x_size * y_size , ncores );
359 else
360 snprintf( name , 64 , "fft_smp_no_place_%d_%d_%d", M , x_size * y_size , ncores );
361 }
362
363 if( EXPLICIT_PLACEMENT )
364 {
365 printf("\n[fft] starts / %d points / %d thread(s) / PID %x / EXPLICIT\n",
366 N, nthreads, pid );
367
368 // build instrumentation file name
369 if( USE_DQT_BARRIER )
370 snprintf( name , 64 , "fft_dqt_explicit_%d_%d_%d", M , x_size * y_size , ncores );
371 else
372 snprintf( name , 64 , "fft_smp_explicit_%d_%d_%d", M , x_size * y_size , ncores );
373 }
374
375 if( PARALLEL_PLACEMENT )
376 {
377 printf("\n[fft] starts / %d points / %d thread(s) / PID %x / PARALLEL\n",
378 N, nthreads, pid );
379
380 // build instrumentation file name
381 if( USE_DQT_BARRIER )
382 snprintf( name , 64 , "fft_dqt_parallel_%d_%d_%d", M , x_size * y_size , ncores );
383 else
384 snprintf( name , 64 , "fft_smp_parallel_%d_%d_%d", M , x_size * y_size , ncores );
385 }
386
387 // build instrumentation file pathname
388 snprintf( path , 128 , "/home/%s", name );
389
390 // open instrumentation file
391 FILE * f = fopen( path , NULL );
392 if ( f == NULL )
393 {
394 printf("\n[fft error] cannot open instrumentation file <%s>\n", path );
395 exit( 0 );
396 }
397
398#if DEBUG_MAIN
399get_cycle( &debug_cycle );
400printf("\n[fft] main open instrumentation file <%s> at cycle %d\n",
401path, (unsigned int)debug_cycle );
402#endif
403
404#if CHECK
405ck1 = CheckSum();
406#endif
407
408#if PRINT_ARRAY
409printf("\nData values / base = %x\n", &data[0][0] );
410PrintArray( data , N );
411
412printf("\nTwiddle values / base = %x\n", &twid[0][0] );
413PrintArray( twid , N );
414
415SimpleDft( 1 , N , data , 0 , bloup , 0 );
416
417printf("\nExpected results / base = %x\n", &bloup[0][0] );
418PrintArray( bloup , N );
419#endif
420
421 // initialise barrier synchronizing all <work> threads
422 if( USE_DQT_BARRIER )
423 {
424 barrier_attr.x_size = x_size;
425 barrier_attr.y_size = y_size;
426 barrier_attr.nthreads = ncores;
427 error = pthread_barrier_init( &barrier, &barrier_attr , nthreads );
428 }
429 else
430 {
431 error = pthread_barrier_init( &barrier, NULL , nthreads );
432 }
433
434 if( error )
435 {
436 printf("\n[fft error] cannot initialize barrier\n");
437 exit( 0 );
438 }
439
440#if DEBUG_MAIN
441get_cycle( &debug_cycle );
442printf("\n[fft] main completes sequencial initialisation at cycle %d\n",
443(unsigned int)debug_cycle );
444#endif
445
446 // register sequencial time
447 get_cycle( &end_init_cycle );
448 init_time = (unsigned int)(end_init_cycle - start_init_cycle);
449
450 //////////////////
451 if( NO_PLACEMENT )
452 {
453 // the tid value for the main thread is always 0
454 // main thread creates new threads with tid in [1,nthreads-1]
455 unsigned int tid;
456 for ( tid = 0 ; tid < nthreads ; tid++ )
457 {
458 // register tid value in work_args[tid] array
459 work_args[tid].tid = tid;
460
461 // create other threads
462 if( tid > 0 )
463 {
464 if ( pthread_create( &work_trdid[tid],
465 NULL, // no attribute
466 &work,
467 &work_args[tid] ) )
468 {
469 printf("\n[fft error] cannot create thread %d\n", tid );
470 exit( 0 );
471 }
472
473#if DEBUG_MAIN
474printf("\n[fft] main created thread %d\n", tid );
475#endif
476
477 }
478 else
479 {
480 tid_main = 0;
481 }
482 } // end for tid
483
484 // main thread calls itself the execute() function
485 work( &work_args[0] );
486
487 // main thread wait other threads completion
488 for ( tid = 1 ; tid < nthreads ; tid++ )
489 {
490 unsigned int * status;
491
492 // main wait thread[tid] status
493 if ( pthread_join( work_trdid[tid], (void*)(&status)) )
494 {
495 printf("\n[fft error] main cannot join thread %d\n", tid );
496 exit( 0 );
497 }
498
499 // check status
500 if( *status != THREAD_EXIT_SUCCESS )
501 {
502 printf("\n[fft error] thread %x returned failure\n", tid );
503 exit( 0 );
504 }
505
506#if DEBUG_MAIN
507printf("\n[fft] main successfully joined thread %x\n", tid );
508#endif
509
510 } // end for tid
511
512 } // end if no_placement
513
514 ////////////////////////
515 if( EXPLICIT_PLACEMENT )
516 {
517 // main thread places each thread[tid] on a specific core[cxy][lid]
518 // but the actual thread creation is sequencial
519 unsigned int x;
520 unsigned int y;
521 unsigned int l;
522 unsigned int cxy; // cluster identifier
523 unsigned int tid; // thread continuous index
524
525 for( x = 0 ; x < x_size ; x++ )
526 {
527 for( y = 0 ; y < y_size ; y++ )
528 {
529 cxy = HAL_CXY_FROM_XY( x , y );
530 for( l = 0 ; l < ncores ; l++ )
531 {
532 // compute thread continuous index
533 tid = (((x * y_size) + y) * ncores) + l;
534
535 // register tid value in work_args[tid] array
536 work_args[tid].tid = tid;
537
538 // no thread created on the core running the main
539 if( (cxy != cxy_main) || (l != lid_main) )
540 {
541 // define thread attributes
542 work_attr[tid].attributes = PT_ATTR_CLUSTER_DEFINED |
543 PT_ATTR_CORE_DEFINED;
544 work_attr[tid].cxy = cxy;
545 work_attr[tid].lid = l;
546
547 // create thread[tid] on core[cxy][l]
548 if ( pthread_create( &work_trdid[tid],
549 &work_attr[tid],
550 &work,
551 &work_args[tid] ) )
552 {
553 printf("\n[fft error] cannot create thread %d\n", tid );
554 exit( 0 );
555 }
556#if DEBUG_MAIN
557printf("\n[fft] main created thread[%d] on core[%x,%d]\n", tid, cxy, l );
558#endif
559 }
560 else
561 {
562 tid_main = tid;
563 }
564 }
565 }
566 }
567
568 // main thread calls itself the execute() function
569 work( &work_args[tid_main] );
570
571 // main thread wait other threads completion
572 for( tid = 0 ; tid < nthreads ; tid++ )
573 {
574 // no other thread on the core running the main
575 if( tid != tid_main )
576 {
577 unsigned int * status;
578
579 // wait thread[tid]
580 if( pthread_join( work_trdid[tid] , (void*)(&status) ) )
581 {
582 printf("\n[fft error] main cannot join thread %d\n", tid );
583 exit( 0 );
584 }
585
586 // check status
587 if( *status != THREAD_EXIT_SUCCESS )
588 {
589 printf("\n[fft error] thread %d returned failure\n", tid );
590 exit( 0 );
591 }
592#if DEBUG_MAIN
593printf("\n[fft] main joined thread %d on core[%x,%d]\n", tid , cxy , l );
594#endif
595 }
596 }
597 } // end if explicit_placement
598
599 ////////////////////////
600 if( PARALLEL_PLACEMENT )
601 {
602 // create and execute the working threads
603 if( pthread_parallel_create( root_level , &work ) )
604 {
605 printf("\n[fft error] cannot create threads\n");
606 exit( 0 );
607 }
608 }
609
610#if DEBUG_MAIN
611get_cycle( &debug_cycle );
612printf("\n[fft] main resume for instrumentation at cycle %d\n",
613(unsigned int)debug_cycle) ;
614#endif
615
616#if PRINT_ARRAY
617printf("\nData values after FFT:\n");
618PrintArray( data , N );
619#endif
620
621#if CHECK
622ck3 = CheckSum();
623printf("\n*** Results ***\n");
624printf("Checksum difference is %f (%f, %f)\n", ck1 - ck3, ck1, ck3);
625if (fabs(ck1 - ck3) < 0.001) printf("Results OK\n");
626else printf("Results KO\n");
627#endif
628
629 // display header on terminal, and save to file
630 printf("\n----- %s -----\n", name );
631
632 ret = fprintf( f , "\n----- %s -----\n", name );
633 if( ret < 0 )
634 {
635 printf("\n[fft error] cannot write header to file <%s>\n", path );
636 exit(0);
637 }
638
639 // initializes global (all threads) instrumentation values
640 unsigned int time_para = parallel_time[0];
641 unsigned int time_sync = sync_time[0];
642 unsigned int pgfaults_nr = 0;
643 unsigned int pgfaults_cost = 0;
644 unsigned int pgfaults_max = pgfault_max[0];
645
646 // loop on threads to compute global instrumentation results
647 for (tid = 0 ; tid < nthreads ; tid++)
648 {
649 snprintf( string , 256 ,
650 "- tid %d : Seq %d / Para %d / Sync %d / Pgfaults %d ( cost %d / max %d )\n",
651 tid, init_time, parallel_time[tid], sync_time[tid],
652 pgfault_nr[tid], (pgfault_cost[tid] / pgfault_nr[tid]) , pgfault_max[tid] );
653
654 // save to instrumentation file
655 fprintf( f , "%s" , string );
656 if( ret < 0 )
657 {
658 printf("\n[fft error] cannot save thread %d results to file <%s>\n", tid, path );
659 printf("%s", string );
660 exit(0);
661 }
662
663 // compute global values
664 if (parallel_time[tid] > time_para) time_para = parallel_time[tid];
665 if (sync_time[tid] > time_sync) time_sync = sync_time[tid];
666
667 pgfaults_nr += pgfault_nr[tid];
668 pgfaults_cost += pgfault_cost[tid];
669
670 if (pgfault_max[tid] > pgfaults_max) pgfaults_max = pgfault_max[tid];
671 }
672
673 // display global values on terminal and save to file
674 snprintf( string , 256 ,
675 "\nSeq %d / Para %d / Sync %d / Pgfaults %d ( cost %d / max %d )\n",
676 init_time, time_para, time_sync, pgfaults_nr, (pgfaults_cost / pgfaults_nr), pgfaults_max );
677
678 printf("%s", string );
679
680 // save global values to file
681 ret = fprintf( f , "%s", string );
682
683 if( ret < 0 )
684 {
685 printf("\n[fft error] cannot save global results to file <%s>\n", path );
686 exit(0);
687 }
688
689 // close instrumentation file
690 ret = fclose( f );
691
692 if( ret < 0 )
693 {
694 printf("\n[fft error] cannot close file <%s>\n", path );
695 exit(0);
696 }
697
698#if DEBUG_MAIN
699get_cycle( &debug_cycle );
700printf("\n[fft] main exit <%s> at cycle %d\n",
701path, (unsigned int)debug_cycle );
702#endif
703
704 exit( 0 );
705
706 return 0;
707
708} // end main()
709
710/////////////////////////////////////////////////////////////////
711// This function is executed in parallel by all <work> threads.
712/////////////////////////////////////////////////////////////////
713void * work( void * arguments )
714{
715 unsigned int tid; // this thread continuous index
716 unsigned int lid; // core local index
717 unsigned int cid; // cluster continuous index
718
719 unsigned int MyFirst; // index first row allocated to thread
720 unsigned int MyLast; // index last row allocated to thread
721 double * upriv; // private array of FFT coefs
722
723 unsigned long long parallel_start;
724 unsigned long long parallel_stop;
725 unsigned long long barrier_start;
726 unsigned long long barrier_stop;
727
728 get_cycle( &parallel_start );
729
730 // get thread arguments
731 pthread_parallel_work_args_t * args = (pthread_parallel_work_args_t *)arguments;
732
733 tid = args->tid;
734 pthread_barrier_t * parent_barrier = args->barrier;
735
736 // compute lid and cid from tid
737 lid = tid % ncores;
738 cid = tid / ncores;
739
740#if DEBUG_WORK
741printf("\n[fft] %s : thread %d enter / cycle %d\n",
742__FUNCTION__, tid, (unsigned int)parallel_start );
743#endif
744
745 // thread on core 0 allocates memory from the local cluster
746 // for the distributed data[], trans[], twid[] buffers
747 if( lid == 0 )
748 {
749 unsigned int data_size = (N / nclusters) * 2 * sizeof(double);
750
751 data[cid] = (double *)malloc( data_size );
752 if( data[cid] == NULL )
753 {
754 printf("\n[fft_error] in work : cannot allocate data[%d] buffer\n", cid );
755 pthread_barrier_wait( parent_barrier );
756 pthread_exit( NULL );
757 }
758
759 trans[cid] = (double *)malloc( data_size );
760 if( trans[cid] == NULL )
761 {
762 printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid );
763 pthread_barrier_wait( parent_barrier );
764 pthread_exit( NULL );
765 }
766
767 twid[cid] = (double *)malloc( data_size );
768 if( twid[cid] == NULL )
769 {
770 printf("\n[fft_error] in work : cannot allocate twid[%d] buffer\n", cid );
771 pthread_barrier_wait( parent_barrier );
772 pthread_exit( NULL );
773 }
774 }
775
776 // BARRIER to wait distributed buffers allocation
777 get_cycle( &barrier_start );
778 pthread_barrier_wait( &barrier );
779 get_cycle( &barrier_stop );
780 sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
781
782#if DEBUG_WORK
783printf("\n[fft] %s : thread %d exit barrier for buffer allocation / cycle %d\n",
784__FUNCTION__, tid, (unsigned int)barrier_stop );
785#endif
786
787 // all threads contribute to data[] local array initialisation
788 InitD( data , MODE , tid );
789
790 // all threads contribute to data[] local array initialisation
791 InitT( twid , tid );
792
793 // BARRIER to wait distributed buffers initialisation
794 get_cycle( &barrier_start );
795 pthread_barrier_wait( &barrier );
796 get_cycle( &barrier_stop );
797 sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
798
799#if DEBUG_WORK
800printf("\n[fft] %s : thread %d exit barrier for buffer initialisation / cycle %d\n",
801__FUNCTION__, tid, (unsigned int)barrier_stop );
802#endif
803
804 // all threads allocate memory from the local cluster
805 // for the private upriv[] buffer
806 upriv = (double *)malloc( (rootN - 1) * 2 * sizeof(double) );
807 if( upriv == NULL )
808 {
809 printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid );
810 pthread_barrier_wait( parent_barrier );
811 pthread_exit( NULL );
812 }
813
814 // all threads initialise the private upriv[] array
815 InitU( upriv );
816
817 // all threads compute first and last rows handled by the thread
818 MyFirst = rootN * tid / nthreads;
819 MyLast = rootN * (tid + 1) / nthreads;
820
821 // all threads perform forward FFT
822 FFT1D( 1 , data , trans , upriv , twid , tid , MyFirst , MyLast );
823
824#if CHECK
825get_cycle( &barrier_start );
826pthread_barrier_wait( &barrier );
827get_cycle( &barrier_stop );
828sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
829FFT1D( -1 , data , trans , upriv , twid , tid , MyFirst , MyLast );
830#endif
831
832 get_cycle( &parallel_stop );
833
834 // register parallel time in instrumentation counters
835 parallel_time[tid] = (unsigned int)(parallel_stop - parallel_start);
836
837 // get work thread info for page faults
838 thread_info_t info;
839 get_thread_info( &info );
840
841 // register page faults in instrumentation counters
842 pgfault_nr[tid] = info.false_pgfault_nr +
843 info.local_pgfault_nr +
844 info.global_pgfault_nr;
845 pgfault_cost[tid] = info.false_pgfault_cost +
846 info.local_pgfault_cost +
847 info.global_pgfault_cost;
848 pgfault_max[tid] = info.false_pgfault_max +
849 info.local_pgfault_max +
850 info.global_pgfault_max;
851#if DEBUG_WORK
852printf("\n[fft] %s : thread %d completes fft / p_start %d / p_stop %d\n",
853__FUNCTION__, tid, (unsigned int)parallel_start, (unsigned int)parallel_stop );
854#endif
855
856 // work thread signals completion to main
857 pthread_barrier_wait( parent_barrier );
858
859#if DEBUG_WORK
860printf("\n[fft] %s : thread %d exit\n",
861__FUNCTION__, tid );
862#endif
863
864#if DISPLAY_SCHED_AND_VMM
865printf("\n[fft] %s : thread %d exit\n", __FUNCTION__, tid );
866if( lid == 0 ) display_vmm( cxy , getpid() , 0 );
867#endif
868
869 // work thread exit
870 pthread_exit( NULL );
871
872 return NULL;
873
874} // end work()
875
876////////////////////////////////////////////////////////////////////////////////////////
877// This function makes the DFT from the src[nclusters][points_per_cluster] distributed
878// buffer, to the dst[nclusters][points_per_cluster] distributed buffer.
879////////////////////////////////////////////////////////////////////////////////////////
880void SimpleDft( int direction, // 1 direct / -1 reverse
881 unsigned int size, // number of points
882 double ** src, // source distributed buffer
883 unsigned int src_offset, // offset in source array
884 double ** dst, // destination distributed buffer
885 unsigned int dst_offset ) // offset in destination array
886{
887 unsigned int n , k;
888 double phi; // 2*PI*n*k/N
889 double u_r; // cos( phi )
890 double u_c; // sin( phi )
891 double d_r; // Re(data[n])
892 double d_c; // Im(data[n])
893 double accu_r; // Re(accu)
894 double accu_c; // Im(accu)
895 unsigned int c_id; // distributed buffer cluster index
896 unsigned int c_offset; // offset in distributed buffer
897
898 for ( k = 0 ; k < size ; k++ ) // loop on the output data points
899 {
900 // initialise accu
901 accu_r = 0;
902 accu_c = 0;
903
904 for ( n = 0 ; n < size ; n++ ) // loop on the input data points
905 {
906 // compute coef
907 phi = (double)(2*PI*n*k) / size;
908 u_r = cos( phi );
909 u_c = -sin( phi ) * direction;
910
911 // get input data point
912 c_id = (src_offset + n) / (points_per_cluster);
913 c_offset = (src_offset + n) % (points_per_cluster);
914 d_r = src[c_id][2*c_offset];
915 d_c = src[c_id][2*c_offset+1];
916
917 // increment accu
918 accu_r += ((u_r*d_r) - (u_c*d_c));
919 accu_c += ((u_r*d_c) + (u_c*d_r));
920 }
921
922 // scale for inverse DFT
923 if ( direction == -1 )
924 {
925 accu_r /= size;
926 accu_c /= size;
927 }
928
929 // set output data point
930 c_id = (dst_offset + k) / (points_per_cluster);
931 c_offset = (dst_offset + k) % (points_per_cluster);
932 dst[c_id][2*c_offset] = accu_r;
933 dst[c_id][2*c_offset+1] = accu_c;
934 }
935
936} // end SimpleDft()
937
938///////////////////////
939double CheckSum( void )
940{
941 unsigned int i , j;
942 unsigned int c_id;
943 unsigned int c_offset;
944 double cks;
945
946 cks = 0.0;
947 for (j = 0; j < rootN ; j++)
948 {
949 for (i = 0; i < rootN ; i++)
950 {
951 c_id = (rootN * j + i) / (points_per_cluster);
952 c_offset = (rootN * j + i) % (points_per_cluster);
953
954 cks += data[c_id][2*c_offset] + data[c_id][2*c_offset+1];
955 }
956 }
957 return(cks);
958}
959
960//////////////////////////////////////////////////////////////////////////////////////
961// Each working thread <tid> contributes to initialize (rootN / nthreads) rows,
962// in the shared - and distributed - <data> array.
963//////////////////////////////////////////////////////////////////////////////////////
964void InitD(double ** data,
965 unsigned int mode,
966 unsigned int tid )
967{
968 unsigned int i , j;
969 unsigned int c_id;
970 unsigned int c_offset;
971 unsigned int index;
972
973 // compute row_min and row_max
974 unsigned int row_min = tid * rows_per_thread;
975 unsigned int row_max = row_min + rows_per_thread;
976
977 for ( j = row_min ; j < row_max ; j++ ) // loop on rows
978 {
979 for ( i = 0 ; i < rootN ; i++ ) // loop on points in a row
980 {
981 index = j * rootN + i;
982 c_id = index / (points_per_cluster);
983 c_offset = index % (points_per_cluster);
984
985 // complex input signal is random
986 if ( mode == RANDOM )
987 {
988 data[c_id][2*c_offset] = ( (double)rand() ) / 65536;
989 data[c_id][2*c_offset+1] = ( (double)rand() ) / 65536;
990 }
991
992
993 // complex input signal is cos(n/N) / sin(n/N)
994 if ( mode == COSIN )
995 {
996 double phi = (double)( 2 * PI * index) / N;
997 data[c_id][2*c_offset] = cos( phi );
998 data[c_id][2*c_offset+1] = sin( phi );
999 }
1000
1001 // complex input signal is constant
1002 if ( mode == CONSTANT )
1003 {
1004 data[c_id][2*c_offset] = 1.0;
1005 data[c_id][2*c_offset+1] = 0.0;
1006 }
1007 }
1008 }
1009}
1010
1011///////////////////////////////////////////////////////////////////////////////////////
1012// Each working thread <tid> contributes to initialize (rootN / nthreads) rows,
1013// in the shared - and distributed - <twiddle> array.
1014///////////////////////////////////////////////////////////////////////////////////////
1015void InitT( double ** twid,
1016 unsigned int tid )
1017{
1018 unsigned int i, j;
1019 unsigned int index;
1020 unsigned int c_id;
1021 unsigned int c_offset;
1022 double phi;
1023
1024 // compute row_min and row_max
1025 unsigned int row_min = tid * rows_per_thread;
1026 unsigned int row_max = row_min + rows_per_thread;
1027
1028 for ( j = row_min ; j < row_max ; j++ ) // loop on rows
1029 {
1030 for ( i = 0 ; i < rootN ; i++ ) // loop on points in a row
1031 {
1032 index = j * rootN + i;
1033 c_id = index / (points_per_cluster);
1034 c_offset = index % (points_per_cluster);
1035
1036 phi = (double)(2.0 * PI * i * j) / N;
1037 twid[c_id][2*c_offset] = cos( phi );
1038 twid[c_id][2*c_offset+1] = -sin( phi );
1039 }
1040 }
1041}
1042
1043///////////////////////////////////////////////////////////////////////////////////////
1044// Each working thread initialize the private <upriv> array / (rootN - 1) entries.
1045///////////////////////////////////////////////////////////////////////////////////////
1046void InitU( double * upriv )
1047{
1048 unsigned int q;
1049 unsigned int j;
1050 unsigned int base;
1051 unsigned int n1;
1052 double phi;
1053
1054 for (q = 0 ; ((unsigned int)(1 << q) < N) ; q++)
1055 {
1056 n1 = 1 << q; // n1 == 2**q
1057 base = n1 - 1;
1058 for (j = 0; (j < n1) ; j++)
1059 {
1060 if (base + j > rootN - 1) return;
1061
1062 phi = (double)(2.0 * PI * j) / (2 * n1);
1063 upriv[2*(base+j)] = cos( phi );
1064 upriv[2*(base+j)+1] = -sin( phi );
1065 }
1066 }
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////////////
1070// This function returns an index value that is the bit reverse of the input value.
1071////////////////////////////////////////////////////////////////////////////////////////
1072unsigned int BitReverse( unsigned int k )
1073{
1074 unsigned int i;
1075 unsigned int j;
1076 unsigned int tmp;
1077
1078 j = 0;
1079 tmp = k;
1080 for (i = 0; i < M/2 ; i++)
1081 {
1082 j = 2 * j + (tmp & 0x1);
1083 tmp = tmp >> 1;
1084 }
1085 return j;
1086}
1087
1088////////////////////////////////////////////////////////////////////////////////////////
1089// This function perform the in place (direct or inverse) FFT on the N data points
1090// contained in the distributed buffers x[nclusters][points_per_cluster].
1091// It handles the (N) points 1D array as a (rootN*rootN) points 2D array.
1092// 1) it fft (rootN/nthreads ) rows from x to tmp.
1093// 2) it make (rootN/nthreads) FFT on the tmp rows and apply the twiddle factor.
1094// 3) it fft (rootN/nthreads) columns from tmp to x.
1095// 4) it make (rootN/nthreads) FFT on the x rows.
1096// It calls the FFTRow() 2*(rootN/nthreads) times to perform the in place FFT
1097// on the rootN points contained in a row.
1098////////////////////////////////////////////////////////////////////////////////////////
1099void FFT1D( int direction, // direct 1 / inverse -1
1100 double ** x, // input & output distributed data points array
1101 double ** tmp, // auxiliary distributed data points array
1102 double * upriv, // local array containing coefs for rootN FFT
1103 double ** twid, // distributed arrays containing N twiddle factors
1104 unsigned int tid, // thread continuous index
1105 unsigned int MyFirst,
1106 unsigned int MyLast )
1107{
1108 unsigned int j;
1109 unsigned long long barrier_start;
1110 unsigned long long barrier_stop;
1111
1112#if DEBUG_FFT1D
1113unsigned long long cycle;
1114get_cycle( &cycle );
1115printf("\n[fft] %s : thread %d enter / first %d / last %d / cycle %d\n",
1116__FUNCTION__, tid, MyFirst, MyLast, (unsigned int)cycle );
1117#endif
1118
1119 // fft (rootN/nthreads) rows from x to tmp
1120 Transpose( x , tmp , MyFirst , MyLast );
1121
1122#if( DEBUG_FFT1D & 1 )
1123get_cycle( &cycle );
1124printf("\n[fft] %s : thread %d after first fft / cycle %d\n",
1125__FUNCTION__, tid, (unsigned int)cycle );
1126if( PRINT_ARRAY ) PrintArray( tmp , N );
1127#endif
1128
1129 // BARRIER
1130 get_cycle( &barrier_start );
1131 pthread_barrier_wait( &barrier );
1132 get_cycle( &barrier_stop );
1133 sync_time[tid] = (unsigned int)(barrier_stop - barrier_start);
1134
1135#if( DEBUG_FFT1D & 1 )
1136get_cycle( &cycle );
1137printf("\n[fft] %s : thread %d exit barrier after first fft / cycle %d\n",
1138__FUNCTION__, tid, (unsigned int)cycle );
1139#endif
1140
1141 // do FFTs on rows of tmp (i.e. columns of x) and apply twiddle factor
1142 for (j = MyFirst; j < MyLast; j++)
1143 {
1144 FFTRow( direction , upriv , tmp , j * rootN );
1145
1146 TwiddleOneCol( direction , j , twid , tmp , j * rootN );
1147 }
1148
1149#if( DEBUG_FFT1D & 1 )
1150printf("\n[fft] %s : thread %d after first twiddle\n", __FUNCTION__, tid);
1151if( PRINT_ARRAY ) PrintArray( tmp , N );
1152#endif
1153
1154 // BARRIER
1155 get_cycle( &barrier_start );
1156 pthread_barrier_wait( &barrier );
1157 get_cycle( &barrier_stop );
1158
1159#if( DEBUG_FFT1D & 1 )
1160printf("\n[fft] %s : thread %d exit barrier after first twiddle\n", __FUNCTION__, tid);
1161#endif
1162
1163 sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1164
1165 // fft tmp to x
1166 Transpose( tmp , x , MyFirst , MyLast );
1167
1168#if( DEBUG_FFT1D & 1 )
1169printf("\n[fft] %s : thread %d after second fft\n", __FUNCTION__, tid);
1170if( PRINT_ARRAY ) PrintArray( x , N );
1171#endif
1172
1173 // BARRIER
1174 get_cycle( &barrier_start );
1175 pthread_barrier_wait( &barrier );
1176 get_cycle( &barrier_stop );
1177
1178#if( DEBUG_FFT1D & 1 )
1179printf("\n[fft] %s : thread %d exit barrier after second fft\n", __FUNCTION__, tid);
1180#endif
1181
1182 sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1183
1184 // do FFTs on rows of x and apply the scaling factor
1185 for (j = MyFirst; j < MyLast; j++)
1186 {
1187 FFTRow( direction , upriv , x , j * rootN );
1188 if (direction == -1) Scale( x , j * rootN );
1189 }
1190
1191#if( DEBUG_FFT1D & 1 )
1192printf("\n[fft] %s : thread %d after FFT on rows\n", __FUNCTION__, tid);
1193if( PRINT_ARRAY ) PrintArray( x , N );
1194#endif
1195
1196 // BARRIER
1197 get_cycle( &barrier_start );
1198 pthread_barrier_wait( &barrier );
1199 get_cycle( &barrier_stop );
1200
1201#if( DEBUG_FFT1D & 1 )
1202printf("\n[fft] %s : thread %d exit barrier after FFT on rows\n", __FUNCTION__, tid);
1203#endif
1204 sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1205
1206 // fft x to tmp
1207 Transpose( x , tmp , MyFirst , MyLast );
1208
1209#if( DEBUG_FFT1D & 1 )
1210printf("\n[fft] %s : thread %x after third fft\n", __FUNCTION__, tid);
1211if( PRINT_ARRAY ) PrintArray( x , N );
1212#endif
1213
1214 // BARRIER
1215 get_cycle( &barrier_start );
1216 pthread_barrier_wait( &barrier );
1217 get_cycle( &barrier_stop );
1218
1219#if( DEBUG_FFT1D & 1 )
1220printf("\n[fft] %s : thread %d exit barrier after third fft\n", __FUNCTION__, tid);
1221#endif
1222
1223 sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1224 sync_time[tid] += (long)(barrier_stop - barrier_start);
1225
1226 // copy tmp to x
1227 Copy( tmp , x , MyFirst , MyLast );
1228
1229#if DEBUG_FFT1D
1230printf("\n[fft] %s : thread %d completed\n", __FUNCTION__, tid);
1231if( PRINT_ARRAY ) PrintArray( x , N );
1232#endif
1233
1234} // end FFT1D()
1235
1236/////////////////////////////////////////////////////////////////////////////////////
1237// This function multiply all points contained in a row (rootN points) of the
1238// x[] array by the corresponding twiddle factor, contained in the u[] array.
1239/////////////////////////////////////////////////////////////////////////////////////
1240void TwiddleOneCol( int direction,
1241 unsigned int j, // y coordinate in 2D view of coef array
1242 double ** u, // coef array base address
1243 double ** x, // data array base address
1244 unsigned int offset_x ) // first point in N points data array
1245{
1246 unsigned int i;
1247 double omega_r;
1248 double omega_c;
1249 double x_r;
1250 double x_c;
1251 unsigned int c_id;
1252 unsigned int c_offset;
1253
1254 for (i = 0; i < rootN ; i++) // loop on the rootN points
1255 {
1256 // get coef
1257 c_id = (j * rootN + i) / (points_per_cluster);
1258 c_offset = (j * rootN + i) % (points_per_cluster);
1259 omega_r = u[c_id][2*c_offset];
1260 omega_c = direction * u[c_id][2*c_offset+1];
1261
1262 // access data
1263 c_id = (offset_x + i) / (points_per_cluster);
1264 c_offset = (offset_x + i) % (points_per_cluster);
1265 x_r = x[c_id][2*c_offset];
1266 x_c = x[c_id][2*c_offset+1];
1267
1268 x[c_id][2*c_offset] = omega_r*x_r - omega_c * x_c;
1269 x[c_id][2*c_offset+1] = omega_r*x_c + omega_c * x_r;
1270 }
1271} // end TwiddleOneCol()
1272
1273////////////////////////////
1274void Scale( double ** x, // data array base address
1275 unsigned int offset_x ) // first point of the row to be scaled
1276{
1277 unsigned int i;
1278 unsigned int c_id;
1279 unsigned int c_offset;
1280
1281 for (i = 0; i < rootN ; i++)
1282 {
1283 c_id = (offset_x + i) / (points_per_cluster);
1284 c_offset = (offset_x + i) % (points_per_cluster);
1285 x[c_id][2*c_offset] /= N;
1286 x[c_id][2*c_offset + 1] /= N;
1287 }
1288}
1289
1290///////////////////////////////////
1291void Transpose( double ** src, // source buffer (array of pointers)
1292 double ** dest, // destination buffer (array of pointers)
1293 unsigned int MyFirst, // first row allocated to the thread
1294 unsigned int MyLast ) // last row allocated to the thread
1295{
1296 unsigned int row; // row index
1297 unsigned int point; // data point index in a row
1298
1299 unsigned int index_src; // absolute index in the source N points array
1300 unsigned int c_id_src; // cluster for the source buffer
1301 unsigned int c_offset_src; // offset in the source buffer
1302
1303 unsigned int index_dst; // absolute index in the dest N points array
1304 unsigned int c_id_dst; // cluster for the dest buffer
1305 unsigned int c_offset_dst; // offset in the dest buffer
1306
1307
1308 // scan all data points allocated to the thread
1309 // (between MyFirst row and MyLast row) from the source buffer
1310 // and write these points to the destination buffer
1311 for ( row = MyFirst ; row < MyLast ; row++ ) // loop on the rows
1312 {
1313 for ( point = 0 ; point < rootN ; point++ ) // loop on points in row
1314 {
1315 index_src = row * rootN + point;
1316 c_id_src = index_src / (points_per_cluster);
1317 c_offset_src = index_src % (points_per_cluster);
1318
1319 index_dst = point * rootN + row;
1320 c_id_dst = index_dst / (points_per_cluster);
1321 c_offset_dst = index_dst % (points_per_cluster);
1322
1323 dest[c_id_dst][2*c_offset_dst] = src[c_id_src][2*c_offset_src];
1324 dest[c_id_dst][2*c_offset_dst+1] = src[c_id_src][2*c_offset_src+1];
1325 }
1326 }
1327} // end Transpose()
1328
1329//////////////////////////////
1330void Copy( double ** src, // source buffer (array of pointers)
1331 double ** dest, // destination buffer (array of pointers)
1332 unsigned int MyFirst, // first row allocated to the thread
1333 unsigned int MyLast ) // last row allocated to the thread
1334{
1335 unsigned int row; // row index
1336 unsigned int point; // data point index in a row
1337
1338 unsigned int index; // absolute index in the N points array
1339 unsigned int c_id; // cluster index
1340 unsigned int c_offset; // offset in local buffer
1341
1342 // scan all data points allocated to the thread
1343 for ( row = MyFirst ; row < MyLast ; row++ ) // loop on the rows
1344 {
1345 for ( point = 0 ; point < rootN ; point++ ) // loop on points in row
1346 {
1347 index = row * rootN + point;
1348 c_id = index / (points_per_cluster);
1349 c_offset = index % (points_per_cluster);
1350
1351 dest[c_id][2*c_offset] = src[c_id][2*c_offset];
1352 dest[c_id][2*c_offset+1] = src[c_id][2*c_offset+1];
1353 }
1354 }
1355} // end Copy()
1356
1357///////////////////////////////
1358void Reverse( double ** x,
1359 unsigned int offset_x )
1360{
1361 unsigned int j, k;
1362 unsigned int c_id_j;
1363 unsigned int c_offset_j;
1364 unsigned int c_id_k;
1365 unsigned int c_offset_k;
1366
1367 for (k = 0 ; k < rootN ; k++)
1368 {
1369 j = BitReverse( k );
1370 if (j > k)
1371 {
1372 c_id_j = (offset_x + j) / (points_per_cluster);
1373 c_offset_j = (offset_x + j) % (points_per_cluster);
1374 c_id_k = (offset_x + k) / (points_per_cluster);
1375 c_offset_k = (offset_x + k) % (points_per_cluster);
1376
1377 SWAP(x[c_id_j][2*c_offset_j] , x[c_id_k][2*c_offset_k]);
1378 SWAP(x[c_id_j][2*c_offset_j+1], x[c_id_k][2*c_offset_k+1]);
1379 }
1380 }
1381}
1382
1383/////////////////////////////////////////////////////////////////////////////
1384// This function makes the in-place FFT on all points contained in a row
1385// (i.e. rootN points) of the x[nclusters][points_per_cluster] array.
1386/////////////////////////////////////////////////////////////////////////////
1387void FFTRow( int direction, // 1 direct / -1 inverse
1388 double * u, // private coefs array
1389 double ** x, // array of pointers on distributed buffers
1390 unsigned int offset_x ) // absolute offset in the x array
1391{
1392 unsigned int j;
1393 unsigned int k;
1394 unsigned int q;
1395 unsigned int L;
1396 unsigned int r;
1397 unsigned int Lstar;
1398 double * u1;
1399
1400 unsigned int offset_x1; // index first butterfly input
1401 unsigned int offset_x2; // index second butterfly output
1402
1403 double omega_r; // real part butterfy coef
1404 double omega_c; // complex part butterfly coef
1405
1406 double tau_r;
1407 double tau_c;
1408
1409 double d1_r; // real part first butterfly input
1410 double d1_c; // imag part first butterfly input
1411 double d2_r; // real part second butterfly input
1412 double d2_c; // imag part second butterfly input
1413
1414 unsigned int c_id_1; // cluster index for first butterfly input
1415 unsigned int c_offset_1; // offset for first butterfly input
1416 unsigned int c_id_2; // cluster index for second butterfly input
1417 unsigned int c_offset_2; // offset for second butterfly input
1418
1419#if DEBUG_ROW
1420unsigned int p;
1421printf("\n[fft] ROW data in / %d points / offset = %d\n", rootN , offset_x );
1422
1423for ( p = 0 ; p < rootN ; p++ )
1424{
1425 unsigned int index = offset_x + p;
1426 unsigned int c_id = index / (points_per_cluster);
1427 unsigned int c_offset = index % (points_per_cluster);
1428 printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1429}
1430printf("\n");
1431#endif
1432
1433 // This makes the rootN input points reordering
1434 Reverse( x , offset_x );
1435
1436#if DEBUG_ROW
1437printf("\n[fft] ROW data after reverse / %d points / offset = %d\n", rootN , offset_x );
1438
1439for ( p = 0 ; p < rootN ; p++ )
1440{
1441 unsigned int index = offset_x + p;
1442 unsigned int c_id = index / (points_per_cluster);
1443 unsigned int c_offset = index % (points_per_cluster);
1444 printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1445}
1446printf("\n");
1447#endif
1448
1449 // This implements the multi-stages, in place Butterfly network
1450 for (q = 1; q <= M/2 ; q++) // loop on stages
1451 {
1452 L = 1 << q; // number of points per subset for current stage
1453 r = rootN / L; // number of subsets
1454 Lstar = L / 2;
1455 u1 = &u[2 * (Lstar - 1)];
1456 for (k = 0; k < r; k++) // loop on the subsets
1457 {
1458 offset_x1 = offset_x + (k * L); // index first point
1459 offset_x2 = offset_x + (k * L + Lstar); // index second point
1460
1461#if (DEBUG_ROW & 1)
1462printf("\n ### q = %d / k = %d / x1 = %d / x2 = %d\n", q , k , offset_x1 , offset_x2 );
1463#endif
1464 // makes all in-place butterfly(s) for subset
1465 for (j = 0; j < Lstar; j++)
1466 {
1467 // get coef
1468 omega_r = u1[2*j];
1469 omega_c = direction * u1[2*j+1];
1470
1471 // get d[x1] address and value
1472 c_id_1 = (offset_x1 + j) / (points_per_cluster);
1473 c_offset_1 = (offset_x1 + j) % (points_per_cluster);
1474 d1_r = x[c_id_1][2*c_offset_1];
1475 d1_c = x[c_id_1][2*c_offset_1+1];
1476
1477 // get d[x2] address and value
1478 c_id_2 = (offset_x2 + j) / (points_per_cluster);
1479 c_offset_2 = (offset_x2 + j) % (points_per_cluster);
1480 d2_r = x[c_id_2][2*c_offset_2];
1481 d2_c = x[c_id_2][2*c_offset_2+1];
1482
1483#if (DEBUG_ROW & 1)
1484printf("\n ### d1_in = (%f , %f) / d2_in = (%f , %f) / coef = (%f , %f)\n",
1485 d1_r , d1_c , d2_r , d2_c , omega_r , omega_c);
1486#endif
1487 // tau = omega * d[x2]
1488 tau_r = omega_r * d2_r - omega_c * d2_c;
1489 tau_c = omega_r * d2_c + omega_c * d2_r;
1490
1491 // set new value for d[x1] = d[x1] + omega * d[x2]
1492 x[c_id_1][2*c_offset_1] = d1_r + tau_r;
1493 x[c_id_1][2*c_offset_1+1] = d1_c + tau_c;
1494
1495 // set new value for d[x2] = d[x1] - omega * d[x2]
1496 x[c_id_2][2*c_offset_2] = d1_r - tau_r;
1497 x[c_id_2][2*c_offset_2+1] = d1_c - tau_c;
1498
1499#if (DEBUG_ROW & 1)
1500printf("\n ### d1_out = (%f , %f) / d2_out = (%f , %f)\n",
1501 d1_r + tau_r , d1_c + tau_c , d2_r - tau_r , d2_c - tau_c );
1502#endif
1503 }
1504 }
1505 }
1506
1507#if DEBUG_ROW
1508printf("\n[fft] ROW data out / %d points / offset = %d\n", rootN , offset_x );
1509for ( p = 0 ; p < rootN ; p++ )
1510{
1511 unsigned int index = offset_x + p;
1512 unsigned int c_id = index / (points_per_cluster);
1513 unsigned int c_offset = index % (points_per_cluster);
1514 printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1515}
1516printf("\n");
1517#endif
1518
1519} // end FFTRow()
1520
1521///////////////////////////////////////
1522void PrintArray( double ** array,
1523 unsigned int size )
1524{
1525 unsigned int i;
1526 unsigned int c_id;
1527 unsigned int c_offset;
1528
1529 // float display
1530 for (i = 0; i < size ; i++)
1531 {
1532 c_id = i / (points_per_cluster);
1533 c_offset = i % (points_per_cluster);
1534
1535 printf(" %f %f |", array[c_id][2*c_offset], array[c_id][2*c_offset+1]);
1536
1537 if ( (i+1) % 4 == 0) printf("\n");
1538 }
1539 printf("\n");
1540}
1541
1542
1543// Local Variables:
1544// tab-width: 4
1545// c-basic-offset: 4
1546// c-file-offsets:((innamespace . 0)(inline-open . 0))
1547// indent-tabs-mode: nil
1548// End:
1549
1550// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4
1551
Note: See TracBrowser for help on using the repository browser.