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

Last change on this file since 582 was 582, checked in by alain, 3 years ago

New DQDT implementation supporting missing clusters
thanks to the cluster_info[x][y] array.

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