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

Last change on this file since 495 was 473, checked in by alain, 6 years ago

Fix several GCC warning related to the -Wextra compilation option.

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