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

Last change on this file was 824, checked in by alain, 8 years ago

bloup

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