source: soft/giet_vm/applications/fft/fft.c

/*************************************************************************/
/*                                                                       */
/*  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 is the port of the SPLASH FFT benchmark on the GIET_VM OS.
// It has been done by Hao LIU (04/12/2014), and Alain Greiner (14/03/2016).
//
// 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 sub-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 two main parameters are the following:      
//  - M : N = 2**M = number of data points / M must be an even number. 
//  - T : nthreads = nprocs defined by the hardware / must be power of 2. 
//
// Several configuration parameters can be defined below:
//  - VERBOSE : Print out complex data points arrays. 
//  - CHECK : Perform both FFT and inverse FFT to check output/input.
//  - DEBUG : Display intermediate results
///////////////////////////////////////////////////////////////////////////

#include "stdio.h"
#include "stdlib.h"
#include "user_barrier.h"
#include "malloc.h"
#include "math.h"

// constants

#define PI                      3.14159265359
#define PAGE_SIZE               4096
#define CLUSTERS_MAX            256
#define THREADS_MAX             1024
#define RANDOM                  0
#define COSIN                   1
#define CONSTANT                2

// parameters

#define DEFAULT_M               12
#define VERBOSE                 0
#define CHECK                   0
#define DEBUG                   0
#define DETAILED_DEBUG          0
#define MODE                    RANDOM

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

// 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   nprocs;                     // number of processors per cluster
long           nthreads;                   // total number of threads (one thread per proc)
long           nclusters;                  // total number of clusters
long           M = DEFAULT_M;              // log2(number of points)
long           N;                          // number of points (N = 2^M)         
long           rootN;                      // rootN = 2^M/2    
long           rows_per_proc;              // number of data "rows" handled by a single thread
long           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
long           parallel_time[THREADS_MAX]; // total computation time (per thread)
long           sync_time[THREADS_MAX];     // cumulative waiting time in barriers (per thread)
long           init_time;                  // initialisation time (in main)

// lock protecting shared TTY
user_lock_t  tty_lock;

// synchronisation barrier (all threads)
giet_sqt_barrier_t barrier;

// functions declaration
void slave();

double CheckSum(double ** x);

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

void InitU(double ** u);

void InitT(double ** u);

long BitReverse( long k );

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

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

void Scale( double **x, long offset_x );

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

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

void Reverse( double ** x, long offset_x );

void FFT1DOnce( long direction , double * u , double ** x , long offset_x );

void PrintArray( double ** x , long size );

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

///////////////////////////////////////////////////////////////////
// This function execute the sequencial initialisation
// launch the parallel execution, and makes the instrumentation.
///////////////////////////////////////////////////////////////////
__attribute__ ((constructor)) void main()
{
    long    x;             // index for cluster X coordinate
    long    y;             // index for cluster Y coordinate
    long    i;             // cluster index = y * x_size + x
    long    n;             // thread index
    long    start_time;    // starting cycle for main
    double  ck1;           // for input/output checking
    double  ck3;           // for input/output checking

    // FFT application starts
    start_time = giet_proctime();

    // allocate shared TTY 
    giet_tty_alloc( 1 );     
    lock_init( &tty_lock);

    // get platform parameters to compute nthreads & nclusters
    giet_procs_number( &x_size , &y_size , &nprocs );

    giet_pthread_assert( ((nprocs == 1) || (nprocs == 2) || (nprocs == 4)),
                         "[FFT ERROR] number of procs per cluster must be 1, 2 or 4");

    giet_pthread_assert( ((x_size == 1) || (x_size == 2) || (x_size == 4) || 
                          (x_size == 8) || (x_size == 16)),
                         "[FFT ERROR] x_size must be 1,2,4,8,16");

    giet_pthread_assert( ((y_size == 1) || (y_size == 2) || (y_size == 4) || 
                          (y_size == 8) || (y_size == 16)),
                         "[FFT ERROR] y_size must be 1,2,4,8,16");

    nthreads  = x_size * y_size * nprocs;
    nclusters = x_size * y_size;

    // compute various constants depending on N and T 
    N                  = 1 << M;
    rootN              = 1 << (M / 2);
    rows_per_proc      = rootN / nthreads;
    points_per_cluster = N / nclusters;
 
    giet_pthread_assert( (rootN >= nthreads) ,
                         "[FFT ERROR] sqrt(N) must be larger than nthreads");

    giet_tty_printf("\n*** Entering FFT  ***\n"
                    " - %d Complex points\n"
                    " - %d Thread(s)\n",
                    N , nthreads );

    // distributed heap initialisation
    for ( x = 0 ; x < x_size ; x++ ) 
    {
        for ( y = 0 ; y < y_size ; y++ ) 
        {
            heap_init( x , y );
        }
    }

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

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

    if ( CHECK ) 
    {
        ck1 = CheckSum( data );
    }

    if ( VERBOSE ) 
    {
        giet_tty_printf("Data values / base = %x\n", &data[0][0] );
        PrintArray( data , N );

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

        SimpleDft( 1 , N , data , 0 , bloup , 0 );
        giet_tty_printf("Expected results / base = %x\n", &bloup[0][0] );
        PrintArray( bloup , N );
    }

    // allocate thread[] array : thread identidiers defined by the kernel
    pthread_t* thread = malloc( nthreads * sizeof(pthread_t) );

    // allocate tid[] array : continuous thred index defined by the user
    long* tid = malloc( nthreads * sizeof(long) );

    // initialise distributed barrier
    sqt_barrier_init( &barrier, x_size , y_size , nprocs );

    // Initialisation completed
    giet_tty_printf("\n*** Starting Parallel FFT  ***\n");

    // register initialisation time
    init_time = giet_proctime() - start_time;

    // launch other threads to execute the slave() function
    for ( n = 1 ; n < nthreads ; n++ )
    {
        tid[n] = n;
        if ( giet_pthread_create( &thread[n],
                                  NULL,                  // no attribute
                                  &slave,
                                  &tid[n] ) )            // pointer on thread index
        {
            giet_tty_printf("\n[FFT ERROR] creating thread %d\n", tid[n] );
            giet_pthread_exit( NULL );
        }
    }

    // main execute itself the slave() function
    tid[0] = 0;
    slave( &tid[0] );

    // wait other threads completion
    for ( n = 1 ; n < nthreads ; n++ )
    {
        if ( giet_pthread_join( thread[n], NULL ) )
        {
            giet_tty_printf("\n[FFT ERROR] joining thread %d\n", tid[n] );
            giet_pthread_exit( NULL );
        }
    }

    if ( VERBOSE ) 
    {
        giet_tty_printf("Data values after FFT:\n");
        PrintArray( data , N );
    }

    if ( CHECK )
    { 
        ck3 = CheckSum( data );
        giet_tty_printf("\n*** Results ***\n");
        giet_tty_printf("Checksum difference is %f (%f, %f)\n", ck1 - ck3, ck1, ck3);
        if (fabs(ck1 - ck3) < 0.001)  giet_tty_printf("Results OK.\n");
        else                          giet_tty_printf("Results KO.\n");
    }

    // instrumentation (display and save on disk)
    char string[256];

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

    // open instrumentation file
    long fd = giet_fat_open( string , O_CREAT );
    if ( fd < 0 ) 
    { 
        giet_tty_printf("\n[FFT ERROR] cannot open instrumentation file %s\n", string );
        giet_pthread_exit( NULL );
    }

    snprintf( string , 256 , "\n*** FFT instrumentation : (%dx%dx%d) procs / %d data points ---\n",
                             x_size, y_size, nprocs , N ); 

    giet_tty_printf( "%s" , string );
    giet_fat_fprintf( fd , "%s" , string );

    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 (i = 1 ; i < nthreads ; i++) 
    {
        if (parallel_time[i] > max_para)  max_para = parallel_time[i];
        if (parallel_time[i] < min_para)  min_para = parallel_time[i];
        if (sync_time[i]     > max_sync)  max_sync = sync_time[i];
        if (sync_time[i]     < min_sync)  min_sync = sync_time[i];
    }

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

    giet_tty_printf("%s" , string );
    giet_fat_fprintf( fd , "%s" , string );

    // close instrumentation file and exit
    giet_fat_close( fd );

    giet_pthread_exit("main() completed");

    return ;
} // end main()

///////////////////////////////////////////////////////////////
// This function is executed in parallel by all threads.
///////////////////////////////////////////////////////////////
__attribute__ ((constructor)) void slave( long* tid ) 
{
    long     i;
    long     MyNum = *tid;
    long     MyFirst;         // index first row allocated to thread
    long     MyLast;          // index last row allocated to thread
    double * upriv;
    long     c_id;
    long     c_offset;
    long     parallel_start;
    long     barrier_start;

    // initialise instrumentation
    parallel_start   = giet_proctime();
    sync_time[MyNum] = 0;

    // 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
    barrier_start = giet_proctime();
    sqt_barrier_wait( &barrier );
    sync_time[MyNum] = giet_proctime() - barrier_start;

    // perform backward FFT if required
    if ( CHECK ) 
    {
        // BARRIER
        barrier_start = giet_proctime();
        sqt_barrier_wait( &barrier );
        sync_time[MyNum] = giet_proctime() - barrier_start;

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

    // register computation time
    parallel_time[MyNum] = giet_proctime() - parallel_start;

    // exit if MyNum != 0
    if( MyNum ) giet_pthread_exit("slave() completed\n");

}  // 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( long      direction,
                long      size,           // number of points
                double ** src,            // source distributed buffer
                long      src_offset,     // offset in source array
                double ** dst,            // destination distributed buffer
                long      dst_offset )    // offset in destination array
{
    long    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)
    long    c_id;           // distributed buffer cluster index
    long    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      = data[c_id][2*c_offset];
            d_c      = data[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(double ** x) 
{
    long i , j;
    double       cks;
    long c_id;
    long 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 ) 
{
    long    i , j;
    long    c_id;
    long    c_offset;
    long    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 )                
            {
                data[c_id][2*c_offset]   = ( (double)giet_rand() ) / 65536;
                data[c_id][2*c_offset+1] = ( (double)giet_rand() ) / 65536;
            }
            

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

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

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

    for (q = 0 ; ((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 )
{
    long    i, j;
    long    index;
    long    c_id;
    long    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.
////////////////////////////////////////////////////////////////////////////////////////
long BitReverse( long k ) 
{
    long i; 
    long j; 
    long 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 FFT1DOnce() 2*(rootN/nthreads) times to perform the in place FFT
// on the rootN points contained in a row.
////////////////////////////////////////////////////////////////////////////////////////
void FFT1D( long       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
            long       MyNum,
            long       MyFirst, 
            long       MyLast )
{
    long j;
    long barrier_start;

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

#if DEBUG
giet_tty_printf("\n@@@ tmp after first transpose\n");
PrintArray( tmp , N );
#endif

    // BARRIER
    barrier_start = giet_proctime();
    sqt_barrier_wait( &barrier );
    sync_time[MyNum] = giet_proctime() - barrier_start;

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

#if DEBUG
giet_tty_printf("\n@@@ tmp after columns FFT + twiddle \n");
PrintArray( tmp , N );
#endif

    // BARRIER
    barrier_start = giet_proctime();
    sqt_barrier_wait( &barrier );
    sync_time[MyNum] = giet_proctime() - barrier_start;

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

#if DEBUG
giet_tty_printf("\n@@@ x after second transpose \n");
PrintArray( x , N );
#endif

    // BARRIER
    barrier_start = giet_proctime();
    sqt_barrier_wait( &barrier );
    sync_time[MyNum] = giet_proctime() - barrier_start;

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

#if DEBUG
giet_tty_printf("\n@@@ x after rows FFT + scaling \n");
PrintArray( x , N );
#endif

    // BARRIER
    barrier_start = giet_proctime();
    sqt_barrier_wait( &barrier );
    sync_time[MyNum] = giet_proctime() - barrier_start;

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

#if DEBUG
giet_tty_printf("\n@@@ tmp after third transpose \n");
PrintArray( tmp , N );
#endif

    // BARRIER
    barrier_start = giet_proctime();
    sqt_barrier_wait( &barrier );
    sync_time[MyNum] = giet_proctime() - barrier_start;

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

#if DEBUG
giet_tty_printf("\n@@@ x after final copy \n");
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( long      direction, 
                    long      j,              // y coordinate in 2D view of coef array
                    double ** u,              // coef array base address
                    double ** x,              // data array base address
                    long      offset_x )      // first point in N points data array
{
    long i;
    double       omega_r; 
    double       omega_c; 
    double       x_r; 
    double       x_c;
    long         c_id;
    long         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 
            long      offset_x )   // first point of the row to be scaled
{
    long i;
    long c_id;
    long c_offset;

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

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

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

    long index_dst;         // absolute index in the dest N points array
    long c_id_dst;          // cluster for the dest buffer
    long 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)
           long      MyFirst,  // first row allocated to the thread
           long      MyLast )  // last row allocated to the thread
{
    long row;                  // row index
    long point;                // data point index in a row

    long index;                // absolute index in the N points array
    long c_id;                 // cluster index
    long 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, 
              long      offset_x )
{
    long j, k;
    long c_id_j;
    long c_offset_j;
    long c_id_k;
    long 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 FFT1DOnce( long      direction,  // direct / inverse
                double *  u,          // private coefs array 
                double ** x,          // array of pointers on distributed buffers
                long      offset_x )  // absolute offset in the x array
{
    long     j;
    long     k;
    long     q;
    long     L;
    long     r;
    long     Lstar;
    double * u1; 

    long     offset_x1;     // index first butterfly input
    long     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

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

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

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

#if DETAILED_DEBUG
giet_tty_printf("\n@@@ FFT ROW data after reverse\n");
for ( p = 0 ; p < rootN ; p++ )
{
    long index    = offset_x + p;
    long c_id     = index / (points_per_cluster);
    long c_offset = index % (points_per_cluster);
    giet_tty_printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
}
giet_tty_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 DETAILED_DEBUG
giet_tty_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 DETAILED_DEBUG
giet_tty_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 DETAILED_DEBUG
giet_tty_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 DETAILED_DEBUG
giet_tty_printf("\n@@@ FFT ROW data out\n");
for ( p = 0 ; p < rootN ; p++ )
{
    long index    = offset_x + p;
    long c_id     = index / (points_per_cluster);
    long c_offset = index % (points_per_cluster);
    giet_tty_printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
}
giet_tty_printf("\n");
#endif

}  // end FFT1DOnce()

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

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

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

        if ( (i+1) % 4 == 0)  giet_tty_printf("\n");
    }
    giet_tty_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