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

Last change on this file since 471 was 469, checked in by alain, 6 years ago

1) Introduce the libsemaphore library.
2) Introduce a small libmath library, required by the "fft" application..
3) Introduce the multithreaded "fft" application.
4) Fix a bad synchronisation bug in the Copy-On-Write mechanism.

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