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

/*************************************************************************/
/*                                                                       */
/*  Copyright (c) 1994 Stanford University                               */
/*                                                                       */
/*  All rights reserved.                                                 */
/*                                                                       */
/*  Permission is given to use, copy, and modify this software for any   */
/*  non-commercial purpose as long as this copyright notice is not       */
/*  removed.  All other uses, including redistribution in whole or in    */
/*  part, are forbidden without prior written permission.                */
/*                                                                       */
/*  This software is provided with absolutely no warranty and no         */
/*  support.                                                             */
/*                                                                       */
/*************************************************************************/

///////////////////////////////////////////////////////////////////////////
// This port of the SPLASH FFT benchmark on the ALMOS-MKH OS has been
// done by Alain Greiner (august 2018).
//
// This application performs the 1D fast Fourier transfom for an array
// of N complex points, using the Cooley-Tuckey FFT method.
// The N data points are seen as a 2D array (rootN rows * rootN columns).
// Each thread handle (rootN / nthreads) rows. The N input data points
// be initialised in three different modes:
// - CONSTANT : all data points have the same [1,0] value
// - COSIN    : data point n has [cos(n/N) , sin(n/N)] values
// - RANDOM   : data points have pseudo random values
//
// This application uses 4 shared data arrays, that are distributed 
// in all clusters (one buffer per cluster):
// - data[N] contains N input data points, with 2 double per point.
// - trans[N] contains N intermediate data points, 2 double per point.
// - umain[rootN] contains rootN coefs required for a rootN points FFT.
// - twid[N] contains N coefs : exp(2*pi*i*j/N) / i and j in [0,rootN-1].
// For data, trans, twid, each sub-buffer contains (N/nclusters) points.
// For umain, each sub-buffer contains (rootN/nclusters) points.
//
// The main parameters for this generic application are the following:      
//  - M : N = 2**M = number of data points / M must be an even number. 
//  - T : nthreads = ncores defined by the hardware / must be power of 2. 
//
// There is one thread per core.
// The max number of clusters is defined by (X_MAX * Y_MAX).
// The max number of cores per cluster is defined by CORES_MAX.
//
// Several configuration parameters can be defined below:
//  - PRINT_ARRAY : Print out complex data points arrays. 
//  - CHECK       : Perform both FFT and inverse FFT to check output/input.
//  - DEBUG_MAIN  : Display intermediate results in main()
//  - DEBUG_FFT1D : Display intermediate results in FFT1D()
//  - DEBUG_ROW  :
//
// Regarding final instrumentation:
// - the sequencial initialisation time (init_time) is computed
//   by the main thread in the main() function.
// - The parallel execution time (parallel_time[i]) is computed by each
//   thread(i) in the slave() function.
// - The synchronisation time related to the barriers (sync_time[i])
//   is computed by each thread(i) in the slave() function.
// The results are displayed on the TXT terminal, and registered on disk.
///////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <unistd.h>
#include <pthread.h>
#include <almosmkh.h>
#include <hal_macros.h>

// constants

#define PI                      3.14159265359
#define PAGE_SIZE               4096
#define X_MAX                   16              // max number of clusters in a row
#define Y_MAX                   16              // max number of clusters in a column
#define CORES_MAX               4               // max number of cores in a cluster
#define CLUSTERS_MAX            X_MAX * Y_MAX
#define THREADS_MAX             CLUSTERS_MAX * CORES_MAX
#define RANDOM                  0
#define COSIN                   1
#define CONSTANT                2

// parameters

#define DEFAULT_M               14              // 16 K data points
#define MODE                    COSIN
#define CHECK                   0
#define DEBUG_MAIN              0               // trace main() function (detailed if odd)
#define DEBUG_SLAVE             0               // trace slave() function (detailed if odd)
#define DEBUG_FFT1D             0               // trace FFT1D() function (detailed if odd)
#define DEBUG_ROW               0               // trace FFTRow() function (detailed if odd)
#define PRINT_ARRAY             0

// macro to swap two variables
#define SWAP(a,b) { double tmp; tmp = a; a = b; b = tmp; }

/////////////////////////////////////////////////////////////////////////////////////
//             structure containing the arguments for the slave() function
/////////////////////////////////////////////////////////////////////////////////////

typedef struct args_s
{
    unsigned int   tid;                    // thread continuous index
    unsigned int   main_tid;               // main thread continuous index
}
args_t;

/////////////////////////////////////////////////////////////////////////////////////
//             global variables
/////////////////////////////////////////////////////////////////////////////////////

unsigned int   x_size;                     // number of clusters per row in the mesh
unsigned int   y_size;                     // number of clusters per column in the mesh
unsigned int   ncores;                     // number of cores per cluster
unsigned int   nthreads;                   // total number of threads (one thread per core)
unsigned int   nclusters;                  // total number of clusters
unsigned int   M = DEFAULT_M;              // log2(number of points)
unsigned int   N;                          // number of points (N = 2^M)         
unsigned int   rootN;                      // rootN = 2^M/2    
unsigned int   rows_per_thread;            // number of data "rows" handled by a single thread
unsigned int   points_per_cluster;         // number of data points per cluster 

// arrays of pointers on distributed buffers (one sub-buffer per cluster) 
double *       data[CLUSTERS_MAX];         // original time-domain data
double *       trans[CLUSTERS_MAX];        // used as auxiliary space for transpose
double *       bloup[CLUSTERS_MAX];        // used as auxiliary space for DFT
double *       umain[CLUSTERS_MAX];        // roots of unity used fo rootN points FFT    
double *       twid[CLUSTERS_MAX];         // twiddle factor : exp(-2iPI*k*n/N) 

// instrumentation counters
unsigned int   parallel_time[THREADS_MAX]; // total computation time (per thread)
unsigned int   sync_time[THREADS_MAX];     // cumulated waiting time in barriers (per thread)
unsigned int   init_time;                  // initialisation time (in main)

// synchronisation barrier (all threads)
pthread_barrier_t      barrier;
pthread_barrierattr_t  barrierattr;

// threads identifiers, attributes, and arguments 
pthread_t       trdid[THREADS_MAX];        // kernel threads identifiers
pthread_attr_t  attr[THREADS_MAX];         // POSIX thread attributes
args_t          args[THREADS_MAX];         // slave function arguments

/////////////////////////////////////////////////////////////////////////////////
//           functions declaration
/////////////////////////////////////////////////////////////////////////////////

void slave( args_t * args );

double CheckSum( void );

void InitX(double ** x , unsigned int mode);

void InitU(double ** u);

void InitT(double ** u);

unsigned int BitReverse( unsigned int k );

void FFT1D( int          direction,
            double    ** x,
            double    ** tmp,
            double     * upriv, 
            double    ** twid,
            unsigned int MyNum,
            unsigned int MyFirst,
            unsigned int MyLast );

void TwiddleOneCol( int          direction,
                    unsigned int j,
                    double    ** u,
                    double    ** x,
                    unsigned int offset_x );

void Scale( double    ** x,
            unsigned int offset_x );

void Transpose( double    ** src, 
                double    ** dest,
                unsigned int MyFirst,
                unsigned int MyLast );

void Copy( double    ** src,
           double    ** dest,
           unsigned int MyFirst,
           unsigned int MyLast );

void Reverse( double    ** x, 
              unsigned int offset_x );

void FFTRow( int          direction,
                double     * u,
                double    ** x,
                unsigned int offset_x );

void PrintArray( double ** x,
                 unsigned int size );

void SimpleDft( int          direction,
                unsigned int size,
                double    ** src,
                unsigned int src_offset,
                double    ** dst,
                unsigned int dst_offset );

///////////////////////////////////////////////////////////////////
// This main() function execute the sequencial initialisation
// launch the parallel execution, and makes the instrumentation.
///////////////////////////////////////////////////////////////////
void main ( void )
{
    unsigned int        main_cxy;          // main thread cluster
    unsigned int        main_x;            // main thread X coordinate
    unsigned int        main_y;            // main thread y coordinate
    unsigned int        main_lid;          // main thread local core index
    unsigned int        main_tid;          // main thread continuous index

    unsigned int        x;                 // current index for cluster X coordinate
    unsigned int        y;                 // current index for cluster Y coordinate
    unsigned int        lid;               // current index for core in a cluster
    unsigned int        ci;                // continuous cluster index (from x,y)
    unsigned int        cxy;               // hardware specific cluster identifier
    unsigned int        tid;               // continuous thread index

    unsigned long long  start_init_cycle; 
    unsigned long long  end_init_cycle;

#if CHECK
double     ck1;           // for input/output checking
double     ck3;           // for input/output checking
#endif
    
    // get FFT application start cycle
    get_cycle( &start_init_cycle );

    // get platform parameters to compute nthreads & nclusters
    if( get_config( &x_size , &y_size , &ncores ) )
    {
        printf("\n[FFT ERROR] cannot get hardware configuration\n");
        exit( 0 );
    }

    // check ncores
    if( (ncores != 1) && (ncores != 2) && (ncores != 4) )
    {
        printf("\n[FFT ERROR] number of cores per cluster must be 1/2/4\n");
        exit( 0 );
    }

    // check x_size
    if( (x_size != 1) && (x_size != 2) && (x_size != 4) && (x_size != 8) && (x_size != 16) )
    {
        printf("\n[FFT ERROR] x_size must be 1/2/4/8/16\n");
        exit( 0 );
    }

    // check y_size
    if( (y_size != 1) && (y_size != 2) && (y_size != 4) && (y_size != 8) && (y_size != 16) )
    {
        printf("\n[FFT ERROR] y_size must be 1/2/4/8/16\n");
        exit( 0 );
    }

    nthreads  = x_size * y_size * ncores;
    nclusters = x_size * y_size;

    // compute various constants depending on N and T 
    N                  = 1 << M;
    rootN              = 1 << (M / 2);
    rows_per_thread    = rootN / nthreads;
    points_per_cluster = N / nclusters;
 
    // check N versus T
    if( rootN < nthreads )
    {
        printf("\n[FFT ERROR] sqrt(N) must be larger than T\n");
        exit( 0 );
    }

    // get main thread coordinates (main_x, main_y, main_lid)
    get_core( &main_cxy , &main_lid );
    main_x   = HAL_X_FROM_CXY( main_cxy );
    main_y   = HAL_Y_FROM_CXY( main_cxy );
    main_tid = (((main_x * y_size) + main_y) * ncores) + main_lid; 

    printf("\n[FFT] starts on core[%x,%d] / %d complex points / %d thread(s) / PID %x\n",
    main_cxy, main_lid, N, nthreads, getpid() );

    // allocate memory for the distributed data[i], trans[i], umain[i], twid[i] buffers
    // the index (i) is a continuous cluster index
    unsigned int data_size   = (N / nclusters) * 2 * sizeof(double);
    unsigned int coefs_size  = (rootN / nclusters) * 2 * sizeof(double);
    for (x = 0 ; x < x_size ; x++)
    {
        for (y = 0 ; y < y_size ; y++)
        {
            ci         = x * y_size + y;
            cxy        = HAL_CXY_FROM_XY( x , y );
            data[ci]   = (double *)remote_malloc( data_size  , cxy ); 
            trans[ci]  = (double *)remote_malloc( data_size  , cxy ); 
            bloup[ci]  = (double *)remote_malloc( data_size  , cxy ); 
            umain[ci]  = (double *)remote_malloc( coefs_size , cxy ); 
            twid[ci]   = (double *)remote_malloc( data_size  , cxy ); 
        }
    }

    printf("\n[FFT] main completes remote_malloc\n");

    // arrays initialisation
    InitX( data , MODE ); 
    InitU( umain ); 
    InitT( twid );

    printf("\n[FFT] main completes arrays init\n");

#if CHECK 
ck1 = CheckSum();
#endif

#if PRINT_ARRAY 
printf("\nData values / base = %x\n", &data[0][0] );
PrintArray( data , N );

printf("\nTwiddle values / base = %x\n", &twid[0][0] );
PrintArray( twid , N );

SimpleDft( 1 , N , data , 0 , bloup , 0 );

printf("\nExpected results / base = %x\n", &bloup[0][0] );
PrintArray( bloup , N );
#endif

    // initialise barrier
    barrierattr.x_size   = x_size;
    barrierattr.y_size   = y_size;
    barrierattr.nthreads = ncores;
    if( pthread_barrier_init( &barrier, &barrierattr , nthreads) )
    {
        printf("\n[FFT ERROR] cannot initialize barrier\n");
        exit( 0 );
    }

    printf("\n[FFT] main completes barrier init\n");

    // launch other threads to execute the slave() function
    // on cores other than the core running the main thread
    for (x = 0 ; x < x_size ; x++)
    {
        for (y = 0 ; y < y_size ; y++)
        {
            // compute cluster identifier
            cxy = HAL_CXY_FROM_XY( x , y );

            for ( lid = 0 ; lid < ncores ; lid++ )
            {
                // compute thread user index (continuous index)
                tid = (((x * y_size) + y) * ncores) + lid;

                // set thread attributes
                attr[tid].attributes = PT_ATTR_CLUSTER_DEFINED | PT_ATTR_CORE_DEFINED;
                attr[tid].cxy        = cxy;
                attr[tid].lid        = lid;

                // set slave function argument
                args[tid].tid      = tid;
                args[tid].main_tid = main_tid;

                // create thread
                if( tid != main_tid )
                {
                    if ( pthread_create( &trdid[tid],  // pointer on kernel identifier
                                         &attr[tid],   // pointer on thread attributes
                                         &slave,       // pointer on function 
                                         &args[tid]) ) // pointer on function arguments 
                    {
                        printf("\n[FFT ERROR] creating thread %x\n", tid );
                        exit( 0 );
                    }
#if DEBUG_MAIN
unsigned long long debug_cycle; 
get_cycle( &debug_cycle ); 
printf("\n[FFT] main created thread %x on core[%x,%d] / cycle %d\n",
tid, cxy, lid, (unsigned int)debug_cycle );
#endif
                }
            }
        }
    }

    // register sequencial initalisation completion cycle
    get_cycle( &end_init_cycle );
    init_time = (unsigned int)(end_init_cycle - start_init_cycle);

    printf("\n[FFT] main enters parallel execution\n");
    
    // main itself executes the slave() function
    slave( &args[main_tid] );

    // wait other threads completion
    for (x = 0 ; x < x_size ; x++)
    {
        for (y = 0 ; y < y_size ; y++)
        {
            for ( lid = 0 ; lid < ncores ; lid++ )
            {
                // compute thread continuous index
                tid = (((x * y_size) + y) * ncores) + lid;

                if( tid != main_tid )
                {
                    if( pthread_join( trdid[tid] , NULL ) )
                    {
                        printf("\n[FFT ERROR] in main thread joining thread %x\n", tid );
                        exit( 0 );
                    }
                    
#if DEBUG_MAIN
printf("\n[FFT] main thread %d joined thread %d\n", main_tid, tid );
#endif

                }
            }
        }
    }

#if PRINT_ARRAY 
printf("\nData values after FFT:\n");
PrintArray( data , N );
#endif

#if CHECK
ck3 = CheckSum();
printf("\n*** Results ***\n");
printf("Checksum difference is %f (%f, %f)\n", ck1 - ck3, ck1, ck3);
if (fabs(ck1 - ck3) < 0.001)  printf("Results OK\n");
else                          printf("Results KO\n");
#endif

    // instrumentation 
    char string[256];

    snprintf( string , 256 , "/home/fft_%d_%d_%d_%d", x_size , y_size , ncores , N );

    // open instrumentation file
//  FILE * f = fopen( string , NULL );
//  if ( f == NULL ) 
//  { 
//      printf("\n[FFT ERROR] cannot open instrumentation file %s\n", string );
//      exit( 0 );
//  }

    snprintf( string , 256 , "\n[FFT] instrumentation : (%dx%dx%d) threads / %d points\n",
    x_size, y_size, ncores , N ); 

    // display on terminal, and save to instrumentation file
    printf( "%s" , string );
//  fprintf( f , string );

    for (tid = 0 ; tid < nthreads ; tid++) 
    {
        snprintf( string , 256 , "\ntid %d : Init %d / Parallel %d / Sync %d\n",
        tid, init_time, parallel_time[tid], sync_time[tid] );

        // display on terminal, and save to instrumentation file
        printf("%s" , string );
//      fprintf( f , string );
    }

    // close instrumentation file and exit
//  fclose( f );
                             

/*
    long min_para = parallel_time[0];
    long max_para = parallel_time[0];
    long min_sync = sync_time[0];
    long max_sync = sync_time[0];

    for (tid = 1 ; tid < nthreads ; tid++) 
    {
        if (parallel_time[tid] > max_para)  max_para = parallel_time[tid];
        if (parallel_time[tid] < min_para)  min_para = parallel_time[tid];
        if (sync_time[tid]     > max_sync)  max_sync = sync_time[tid];
        if (sync_time[tid]     < min_sync)  min_sync = sync_time[tid];
    }

    snprintf( string , 256 , "\n      Init       Parallel   Barrier\n"
                             "MIN : %d\t | %d\t | %d\t   (cycles)\n" 
                             "MAX : %d\t | %d\t | %d\t   (cycles)\n",
                             (int)init_time, (int)min_para, (int)min_sync,
                             (int)init_time, (int)max_para, (int)max_sync );
*/

    pthread_exit( NULL );

} // end main()

///////////////////////////////////////////////////////////////
// This function is executed in parallel by all threads.
///////////////////////////////////////////////////////////////
void slave( args_t * args ) 
{
    unsigned int   i;
    unsigned int   MyNum;           // this thread index
    unsigned int   MainNum;         // main thread index
    unsigned int   MyFirst;         // index first row allocated to thread
    unsigned int   MyLast;          // index last row allocated to thread
    double       * upriv;
    unsigned int   c_id;
    unsigned int   c_offset;

    unsigned long long  parallel_start;
    unsigned long long  parallel_stop;
    unsigned long long  barrier_start;
    unsigned long long  barrier_stop;

    MyNum   = args->tid;
    MainNum = args->main_tid;

    // initialise instrumentation
    get_cycle( &parallel_start );

#if DEBUG_SLAVE
printf("\n[FFT] %s : thread %x enter / cycle %d\n",
__FUNCTION__, MyNum, (unsigned int)parallel_start );
#endif

    // allocate and initialise local array upriv[] 
    // that is a local copy of the rootN coefs defined in umain[]
    upriv = (double *)malloc(2 * (rootN - 1) * sizeof(double));  
    for ( i = 0 ; i < (rootN - 1) ; i++) 
    {
        c_id     = i / (rootN / nclusters);
        c_offset = i % (rootN / nclusters);
        upriv[2*i]   = umain[c_id][2*c_offset];
        upriv[2*i+1] = umain[c_id][2*c_offset+1];
    }

    // compute first and last rows handled by the thread
    MyFirst = rootN * MyNum / nthreads;
    MyLast  = rootN * (MyNum + 1) / nthreads;

    // perform forward FFT 
    FFT1D( 1 , data , trans , upriv , twid , MyNum , MyFirst , MyLast );

    // BARRIER
    get_cycle( &barrier_start );
    pthread_barrier_wait( &barrier );
    get_cycle( &barrier_stop );

    sync_time[MyNum] += (barrier_stop - barrier_start);

#if CHECK 
get_cycle( &barrier_start );
pthread_barrier_wait( &barrier );
get_cycle( &barrier_stop );

sync_time[MyNum] += (long)(barrier_stop - barrier_start);

FFT1D( -1 , data , trans , upriv , twid , MyNum , MyFirst , MyLast );
#endif

    // register computation time
    get_cycle( &parallel_stop );
    parallel_time[MyNum] = (parallel_stop - parallel_start);

#if DEBUG_SLAVE
printf("\n[FFT] %s : thread %x exit / parallel_time %d / sync_time %d / cycle %d\n",
__FUNCTION__, MyNum, parallel_time[MyNum], sync_time[MyNum], (unsigned int)parallel_stop );
#endif

    // exit only if MyNum != MainNum
    if( MyNum != MainNum ) pthread_exit( NULL );

}  // end slave()

////////////////////////////////////////////////////////////////////////////////////////
// This function makes the DFT from the src[nclusters][points_per_cluster] distributed
// buffer, to the dst[nclusters][points_per_cluster] distributed buffer.
////////////////////////////////////////////////////////////////////////////////////////
void SimpleDft( int             direction,      // 1 direct / -1 reverse
                unsigned int    size,           // number of points
                double       ** src,            // source distributed buffer
                unsigned int    src_offset,     // offset in source array
                double       ** dst,            // destination distributed buffer
                unsigned int    dst_offset )    // offset in destination array
{
    unsigned int  n , k;
    double        phi;            // 2*PI*n*k/N
    double        u_r;            // cos( phi )
    double        u_c;            // sin( phi )
    double        d_r;            // Re(data[n])
    double        d_c;            // Im(data[n])
    double        accu_r;         // Re(accu)
    double        accu_c;         // Im(accu)
    unsigned int  c_id;           // distributed buffer cluster index
    unsigned int  c_offset;       // offset in distributed buffer

    for ( k = 0 ; k < size ; k++ )       // loop on the output data points
    {
        // initialise accu
        accu_r = 0;
        accu_c = 0;

        for ( n = 0 ; n < size ; n++ )   // loop on the input data points
        {
            // compute coef
            phi = (double)(2*PI*n*k) / size;
            u_r =  cos( phi );
            u_c = -sin( phi ) * direction;

            // get input data point
            c_id     = (src_offset + n) / (points_per_cluster);
            c_offset = (src_offset + n) % (points_per_cluster);
            d_r      = src[c_id][2*c_offset];
            d_c      = src[c_id][2*c_offset+1];

            // increment accu
            accu_r += ((u_r*d_r) - (u_c*d_c));
            accu_c += ((u_r*d_c) + (u_c*d_r));
        }

        // scale for inverse DFT
        if ( direction == -1 )
        {
            accu_r /= size;
            accu_c /= size;
        }

        // set output data point
        c_id     = (dst_offset + k) / (points_per_cluster);
        c_offset = (dst_offset + k) % (points_per_cluster);
        dst[c_id][2*c_offset]   = accu_r;
        dst[c_id][2*c_offset+1] = accu_c;
    }

}  // end SimpleDft()

///////////////////////
double CheckSum( void )
{
    unsigned int         i , j;
    double       cks;
    unsigned int         c_id;
    unsigned int         c_offset;

    cks = 0.0;
    for (j = 0; j < rootN ; j++) 
    {
        for (i = 0; i < rootN ; i++) 
        {
            c_id      = (rootN * j + i) / (points_per_cluster);
            c_offset  = (rootN * j + i) % (points_per_cluster);

            cks += data[c_id][2*c_offset] + data[c_id][2*c_offset+1];
        }
    }
    return(cks);
}


////////////////////////////
void InitX(double      ** x,
           unsigned int   mode ) 
{
    unsigned int    i , j;
    unsigned int    c_id;
    unsigned int    c_offset;
    unsigned int    index;

    for ( j = 0 ; j < rootN ; j++ )      // loop on row index
    {  
        for ( i = 0 ; i < rootN ; i++ )  // loop on point in a row
        {  
            index     = j * rootN + i;
            c_id      = index / (points_per_cluster);
            c_offset  = index % (points_per_cluster);

            // complex input signal is random
            if ( mode == RANDOM )                
            {
                x[c_id][2*c_offset]   = ( (double)rand() ) / 65536;
                x[c_id][2*c_offset+1] = ( (double)rand() ) / 65536;
            }
            

            // complex input signal is cos(n/N) / sin(n/N) 
            if ( mode == COSIN )                
            {
                double phi = (double)( 2 * PI * index) / N;
                x[c_id][2*c_offset]   = cos( phi );
                x[c_id][2*c_offset+1] = sin( phi );
            }

            // complex input signal is constant 
            if ( mode == CONSTANT )                
            {
                x[c_id][2*c_offset]   = 1.0;
                x[c_id][2*c_offset+1] = 0.0;
            }
        }
    }
}

/////////////////////////
void InitU( double ** u ) 
{
    unsigned int    q; 
    unsigned int    j; 
    unsigned int    base; 
    unsigned int    n1;
    unsigned int    c_id;
    unsigned int    c_offset;
    double  phi;
    unsigned int    stop = 0;

    for (q = 0 ; ((unsigned int)(1 << q) < N) && (stop == 0) ; q++) 
    {  
        n1 = 1 << q;
        base = n1 - 1;
        for (j = 0; (j < n1) && (stop == 0) ; j++) 
        {
            if (base + j > rootN - 1) return;

            c_id      = (base + j) / (rootN / nclusters);
            c_offset  = (base + j) % (rootN / nclusters);
            phi = (double)(2.0 * PI * j) / (2 * n1);
            u[c_id][2*c_offset]   = cos( phi );
            u[c_id][2*c_offset+1] = -sin( phi );
        }
    }
}

//////////////////////////
void InitT( double ** u )
{
    unsigned int    i, j;
    unsigned int    index;
    unsigned int    c_id;
    unsigned int    c_offset;
    double  phi;

    for ( j = 0 ; j < rootN ; j++ )      // loop on row index
    {  
        for ( i = 0 ; i < rootN ; i++ )  // loop on points in a row
        {  
            index     = j * rootN + i;
            c_id      = index / (points_per_cluster);
            c_offset  = index % (points_per_cluster);

            phi = (double)(2.0 * PI * i * j) / N;
            u[c_id][2*c_offset]   = cos( phi );
            u[c_id][2*c_offset+1] = -sin( phi );
        }
    }
}

////////////////////////////////////////////////////////////////////////////////////////
// This function returns an index value that is the bit reverse of the input value.
////////////////////////////////////////////////////////////////////////////////////////
unsigned int BitReverse( unsigned int k ) 
{
    unsigned int i; 
    unsigned int j; 
    unsigned int tmp;

    j = 0;
    tmp = k;
    for (i = 0; i < M/2 ; i++) 
    {
        j = 2 * j + (tmp & 0x1);
        tmp = tmp >> 1;
    }
    return j;
}

////////////////////////////////////////////////////////////////////////////////////////
// This function perform the in place (direct or inverse) FFT on the N data points
// contained in the distributed buffers x[nclusters][points_per_cluster].
// It handles the (N) points 1D array as a (rootN*rootN) points 2D array. 
// 1) it transpose (rootN/nthreads ) rows from x to tmp.
// 2) it make (rootN/nthreads) FFT on the tmp rows and apply the twiddle factor.
// 3) it transpose (rootN/nthreads) columns from tmp to x.
// 4) it make (rootN/nthreads) FFT on the x rows.
// It calls the FFTRow() 2*(rootN/nthreads) times to perform the in place FFT
// on the rootN points contained in a row.
////////////////////////////////////////////////////////////////////////////////////////
void FFT1D( int              direction,       // direct 1 / inverse -1
            double       **  x,               // input & output distributed data points array
            double       **  tmp,             // auxiliary distributed data points array
            double        *  upriv,           // local array containing coefs for rootN FFT
            double       **  twid,            // distributed arrays containing N twiddle factors
            unsigned int     MyNum,           // thread continuous index
            unsigned int     MyFirst, 
            unsigned int     MyLast )
{
    unsigned int j;
    unsigned long long barrier_start;
    unsigned long long barrier_stop;

#if DEBUG_FFT1D
unsigned long long cycle;
get_cycle( &cycle );
printf("\n[FFT] %s : thread %x enter / first %d / last %d / cycle %d\n",
__FUNCTION__, MyNum, MyFirst, MyLast, (unsigned int)cycle );
#endif

    // transpose (rootN/nthreads) rows from x to tmp 
    Transpose( x , tmp , MyFirst , MyLast );

#if( DEBUG_FFT1D & 1 )
get_cycle( &cycle );
printf("\n[FFT] %s : thread %x after first transpose / cycle %d\n",
__FUNCTION__, MyNum, (unsigned int)cycle );
if( PRINT_ARRAY ) PrintArray( tmp , N );
#endif

    // BARRIER
    get_cycle( &barrier_start );
    pthread_barrier_wait( &barrier );
    get_cycle( &barrier_stop );
    sync_time[MyNum] = (long)(barrier_stop - barrier_start);

#if( DEBUG_FFT1D & 1 )
get_cycle( &cycle );
printf("\n[FFT] %s : thread %x exit barrier after first transpose / cycle %d\n",
__FUNCTION__, MyNum, (unsigned int)cycle );
#endif

    // do FFTs on rows of tmp (i.e. columns of x) and apply twiddle factor
    for (j = MyFirst; j < MyLast; j++) 
    {
        FFTRow( direction , upriv , tmp , j * rootN );

        TwiddleOneCol( direction , j , twid , tmp , j * rootN );
    }  

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x after first twiddle\n", __FUNCTION__, MyNum);
if( PRINT_ARRAY ) PrintArray( tmp , N );
#endif

    // BARRIER
    get_cycle( &barrier_start );
    pthread_barrier_wait( &barrier );
    get_cycle( &barrier_stop );

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x exit barrier after first twiddle\n", __FUNCTION__, MyNum);
#endif

    sync_time[MyNum] += (long)(barrier_stop - barrier_start);

    // transpose tmp to x
    Transpose( tmp , x , MyFirst , MyLast );

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x after second transpose\n", __FUNCTION__, MyNum);
if( PRINT_ARRAY ) PrintArray( x , N );
#endif

    // BARRIER
    get_cycle( &barrier_start );
    pthread_barrier_wait( &barrier );
    get_cycle( &barrier_stop );

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x exit barrier after second transpose\n", __FUNCTION__, MyNum);
#endif

    sync_time[MyNum] += (long)(barrier_stop - barrier_start);

    // do FFTs on rows of x and apply the scaling factor 
    for (j = MyFirst; j < MyLast; j++) 
    {
        FFTRow( direction , upriv , x , j * rootN );
        if (direction == -1) Scale( x , j * rootN );
    }

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x after FFT on rows\n", __FUNCTION__, MyNum);
if( PRINT_ARRAY ) PrintArray( x , N );
#endif

    // BARRIER
    get_cycle( &barrier_start );
    pthread_barrier_wait( &barrier );
    get_cycle( &barrier_stop );

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x exit barrier after FFT on rows\n", __FUNCTION__, MyNum);
#endif
    sync_time[MyNum] += (long)(barrier_stop - barrier_start);

    // transpose x to tmp
    Transpose( x , tmp , MyFirst , MyLast );

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x after third transpose\n", __FUNCTION__, MyNum);
if( PRINT_ARRAY ) PrintArray( x , N );
#endif

    // BARRIER
    get_cycle( &barrier_start );
    pthread_barrier_wait( &barrier );
    get_cycle( &barrier_stop );

#if( DEBUG_FFT1D & 1 )
printf("\n[FFT] %s : thread %x exit barrier after third transpose\n", __FUNCTION__, MyNum);
#endif

    sync_time[MyNum] += (long)(barrier_stop - barrier_start);
    sync_time[MyNum] += (long)(barrier_stop - barrier_start);

    // copy tmp to x
    Copy( tmp , x , MyFirst , MyLast );

#if DEBUG_FFT1D
printf("\n[FFT] %s : thread %x completed\n", __FUNCTION__, MyNum);
if( PRINT_ARRAY ) PrintArray( x , N );
#endif

}  // end FFT1D()

/////////////////////////////////////////////////////////////////////////////////////
// This function multiply all points contained in a row (rootN points) of the 
// x[] array by the corresponding twiddle factor, contained in the u[] array.
/////////////////////////////////////////////////////////////////////////////////////
void TwiddleOneCol( int             direction, 
                    unsigned int    j,              // y coordinate in 2D view of coef array
                    double       ** u,              // coef array base address
                    double       ** x,              // data array base address
                    unsigned int    offset_x )      // first point in N points data array
{
    unsigned int i;
    double       omega_r; 
    double       omega_c; 
    double       x_r; 
    double       x_c;
    unsigned int c_id;
    unsigned int c_offset;

    for (i = 0; i < rootN ; i++)  // loop on the rootN points
    {
        // get coef
        c_id      = (j * rootN + i) / (points_per_cluster);
        c_offset  = (j * rootN + i) % (points_per_cluster);
        omega_r = u[c_id][2*c_offset];
        omega_c = direction * u[c_id][2*c_offset+1];

        // access data
        c_id      = (offset_x + i) / (points_per_cluster);
        c_offset  = (offset_x + i) % (points_per_cluster);   
        x_r = x[c_id][2*c_offset]; 
        x_c = x[c_id][2*c_offset+1];

        x[c_id][2*c_offset]   = omega_r*x_r - omega_c * x_c;
        x[c_id][2*c_offset+1] = omega_r*x_c + omega_c * x_r;
    }
}  // end TwiddleOneCol()

////////////////////////////
void Scale( double      ** x,           // data array base address 
            unsigned int   offset_x )   // first point of the row to be scaled
{
    unsigned int i;
    unsigned int c_id;
    unsigned int c_offset;

    for (i = 0; i < rootN ; i++) 
    {
        c_id      = (offset_x + i) / (points_per_cluster);
        c_offset  = (offset_x + i) % (points_per_cluster);
        x[c_id][2*c_offset]     /= N;
        x[c_id][2*c_offset + 1] /= N;
    }
}

///////////////////////////////////
void Transpose( double      ** src,      // source buffer (array of pointers)
                double      ** dest,     // destination buffer (array of pointers)
                unsigned int   MyFirst,  // first row allocated to the thread
                unsigned int   MyLast )  // last row allocated to the thread
{
    unsigned int row;               // row index
    unsigned int point;             // data point index in a row

    unsigned int index_src;         // absolute index in the source N points array
    unsigned int c_id_src;          // cluster for the source buffer
    unsigned int c_offset_src;      // offset in the source buffer

    unsigned int index_dst;         // absolute index in the dest N points array
    unsigned int c_id_dst;          // cluster for the dest buffer
    unsigned int c_offset_dst;      // offset in the dest buffer

    
    // scan all data points allocated to the thread 
    // (between MyFirst row and MyLast row) from the source buffer
    // and write these points to the destination buffer
    for ( row = MyFirst ; row < MyLast ; row++ )       // loop on the rows
    {
        for ( point = 0 ; point < rootN ; point++ )    // loop on points in row
        {
            index_src    = row * rootN + point;
            c_id_src     = index_src / (points_per_cluster);
            c_offset_src = index_src % (points_per_cluster);

            index_dst    = point * rootN + row;
            c_id_dst     = index_dst / (points_per_cluster);
            c_offset_dst = index_dst % (points_per_cluster);

            dest[c_id_dst][2*c_offset_dst]   = src[c_id_src][2*c_offset_src];
            dest[c_id_dst][2*c_offset_dst+1] = src[c_id_src][2*c_offset_src+1];
        }
    }
}  // end Transpose()

//////////////////////////////
void Copy( double      ** src,      // source buffer (array of pointers)
           double      ** dest,     // destination buffer (array of pointers)
           unsigned int   MyFirst,  // first row allocated to the thread
           unsigned int   MyLast )  // last row allocated to the thread
{
    unsigned int row;                  // row index
    unsigned int point;                // data point index in a row

    unsigned int index;                // absolute index in the N points array
    unsigned int c_id;                 // cluster index
    unsigned int c_offset;             // offset in local buffer

    // scan all data points allocated to the thread 
    for ( row = MyFirst ; row < MyLast ; row++ )       // loop on the rows
    {
        for ( point = 0 ; point < rootN ; point++ )    // loop on points in row
        {
            index    = row * rootN + point;
            c_id     = index / (points_per_cluster);
            c_offset = index % (points_per_cluster);

            dest[c_id][2*c_offset]   = src[c_id][2*c_offset];
            dest[c_id][2*c_offset+1] = src[c_id][2*c_offset+1];
        }
    }
}  // end Copy()

///////////////////////////////
void Reverse( double      ** x, 
              unsigned int   offset_x )
{
    unsigned int j, k;
    unsigned int c_id_j;
    unsigned int c_offset_j;
    unsigned int c_id_k;
    unsigned int c_offset_k;

    for (k = 0 ; k < rootN ; k++) 
    {
        j = BitReverse( k );
        if (j > k) 
        {
            c_id_j      = (offset_x + j) / (points_per_cluster);
            c_offset_j  = (offset_x + j) % (points_per_cluster);
            c_id_k      = (offset_x + k) / (points_per_cluster);
            c_offset_k  = (offset_x + k) % (points_per_cluster);

            SWAP(x[c_id_j][2*c_offset_j]  , x[c_id_k][2*c_offset_k]);
            SWAP(x[c_id_j][2*c_offset_j+1], x[c_id_k][2*c_offset_k+1]);
        }
    }
}

/////////////////////////////////////////////////////////////////////////////
// This function makes the in-place FFT on all points contained in a row
// (i.e. rootN points) of the x[nclusters][points_per_cluster] array.
/////////////////////////////////////////////////////////////////////////////
void FFTRow( int            direction,  // 1 direct / -1 inverse
                double       * u,          // private coefs array 
                double      ** x,          // array of pointers on distributed buffers
                unsigned int   offset_x )  // absolute offset in the x array
{
    unsigned int     j;
    unsigned int     k;
    unsigned int     q;
    unsigned int     L;
    unsigned int     r;
    unsigned int     Lstar;
    double * u1; 

    unsigned int     offset_x1;     // index first butterfly input
    unsigned int     offset_x2;     // index second butterfly output

    double           omega_r;       // real part butterfy coef
    double           omega_c;       // complex part butterfly coef

    double           tau_r;
    double           tau_c;

    double           d1_r;          // real part first butterfly input
    double           d1_c;          // imag part first butterfly input
    double           d2_r;          // real part second butterfly input
    double           d2_c;          // imag part second butterfly input

    unsigned int     c_id_1;        // cluster index for first butterfly input
    unsigned int     c_offset_1;    // offset for first butterfly input
    unsigned int     c_id_2;        // cluster index for second butterfly input
    unsigned int     c_offset_2;    // offset for second butterfly input

#if DEBUG_ROW
unsigned int p;
printf("\n[FFT] ROW data in / %d points / offset = %d\n", rootN , offset_x );

for ( p = 0 ; p < rootN ; p++ )
{
    unsigned int index    = offset_x + p;
    unsigned int c_id     = index / (points_per_cluster);
    unsigned int c_offset = index % (points_per_cluster);
    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
}
printf("\n");
#endif

    // This makes the rootN input points reordering
    Reverse( x , offset_x );  

#if DEBUG_ROW
printf("\n[FFT] ROW data after reverse / %d points / offset = %d\n", rootN , offset_x );

for ( p = 0 ; p < rootN ; p++ )
{
    unsigned int index    = offset_x + p;
    unsigned int c_id     = index / (points_per_cluster);
    unsigned int c_offset = index % (points_per_cluster);
    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
}
printf("\n");
#endif

    // This implements the multi-stages, in place Butterfly network
    for (q = 1; q <= M/2 ; q++)     // loop on stages
    {
        L = 1 << q;       // number of points per subset for current stage
        r = rootN / L;    // number of subsets
        Lstar = L / 2;
        u1 = &u[2 * (Lstar - 1)];
        for (k = 0; k < r; k++)     // loop on the subsets
        {
            offset_x1  = offset_x + (k * L);            // index first point
            offset_x2  = offset_x + (k * L + Lstar);    // index second point

#if (DEBUG_ROW & 1)
printf("\n ### q = %d / k = %d / x1 = %d / x2 = %d\n", q , k , offset_x1 , offset_x2 );
#endif
            // makes all in-place butterfly(s) for subset
            for (j = 0; j < Lstar; j++) 
            {
                // get coef
                omega_r = u1[2*j];
                omega_c = direction * u1[2*j+1];

                // get d[x1] address and value
                c_id_1      = (offset_x1 + j) / (points_per_cluster);
                c_offset_1  = (offset_x1 + j) % (points_per_cluster);
                d1_r        = x[c_id_1][2*c_offset_1];
                d1_c        = x[c_id_1][2*c_offset_1+1];

                // get d[x2] address and value
                c_id_2      = (offset_x2 + j) / (points_per_cluster);
                c_offset_2  = (offset_x2 + j) % (points_per_cluster);
                d2_r        = x[c_id_2][2*c_offset_2];
                d2_c        = x[c_id_2][2*c_offset_2+1];

#if (DEBUG_ROW & 1)
printf("\n ### d1_in = (%f , %f) / d2_in = (%f , %f) / coef = (%f , %f)\n", 
                d1_r , d1_c , d2_r , d2_c , omega_r , omega_c);
#endif
                // tau = omega * d[x2]
                tau_r = omega_r * d2_r - omega_c * d2_c;
                tau_c = omega_r * d2_c + omega_c * d2_r;

                // set new value for d[x1] = d[x1] + omega * d[x2]
                x[c_id_1][2*c_offset_1]   = d1_r + tau_r;
                x[c_id_1][2*c_offset_1+1] = d1_c + tau_c;

                // set new value for d[x2] = d[x1] - omega * d[x2]
                x[c_id_2][2*c_offset_2]   = d1_r - tau_r;
                x[c_id_2][2*c_offset_2+1] = d1_c - tau_c;

#if (DEBUG_ROW & 1)
printf("\n ### d1_out = (%f , %f) / d2_out = (%f , %f)\n", 
                d1_r + tau_r , d1_c + tau_c , d2_r - tau_r , d2_c - tau_c );
#endif
            }
        }
    }

#if DEBUG_ROW
printf("\n[FFT] ROW data out / %d points / offset = %d\n", rootN , offset_x );
for ( p = 0 ; p < rootN ; p++ )
{
    unsigned int index    = offset_x + p;
    unsigned int c_id     = index / (points_per_cluster);
    unsigned int c_offset = index % (points_per_cluster);
    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
}
printf("\n");
#endif

}  // end FFTRow()

///////////////////////////////////////
void PrintArray( double       ** array,
                 unsigned int    size ) 
{
    unsigned int  i;
    unsigned int  c_id;
    unsigned int  c_offset;

    // float display
    for (i = 0; i < size ; i++) 
    {
        c_id      = i / (points_per_cluster);
        c_offset  = i % (points_per_cluster);

        printf(" %f  %f |", array[c_id][2*c_offset], array[c_id][2*c_offset+1]);

        if ( (i+1) % 4 == 0)  printf("\n");
    }
    printf("\n");
}


// Local Variables:
// tab-width: 4
// c-basic-offset: 4
// c-file-offsets:((innamespace . 0)(inline-open . 0))
// indent-tabs-mode: nil
// End:

// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4