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

Last change on this file since 591 was 588, checked in by alain, 6 years ago

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