source: trunk/user/fft/fft.c @ 677

Last change on this file since 677 was 659, checked in by alain, 4 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 = (((* 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.