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

Last change on this file since 640 was 640, checked in by alain, 4 years ago

Remove all RPCs in page-fault handling.

File size: 47.4 KB
Line 
1/*************************************************************************/
2/*                                                                       */
3/*  Copyright (c) 1994 Stanford University                               */
4/*                                                                       */
5/*  All rights reserved.                                                 */
6/*                                                                       */
7/*  Permission is given to use, copy, and modify this software for any   */
8/*  non-commercial purpose as long as this copyright notice is not       */
9/*  removed.  All other uses, including redistribution in whole or in    */
10/*  part, are forbidden without prior written permission.                */
11/*                                                                       */
12/*  This software is provided with absolutely no warranty and no         */
13/*  support.                                                             */
14/*                                                                       */
15/*************************************************************************/
16
17///////////////////////////////////////////////////////////////////////////
18// This port of the SPLASH FFT benchmark on the ALMOS-MKH OS has been
19// done by Alain Greiner (august 2018).
20//
21// This application performs the 1D fast Fourier transfom for an array
22// of N complex points, using the Cooley-Tuckey FFT method.
23// The N data points are seen as a 2D array (rootN rows * rootN columns).
24// Each thread handle (rootN / nthreads) rows.
25// The N input data points can 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// The main parameters for this generic application are the following:     
31//  - M : N = 2**M = number of data points / M must be an even number.
32//  - T : nthreads = ncores defined by the hardware / must be power of 2.
33// The number of threads cannot be larger than the number of rows.
34//
35// This application uses 3 shared data arrays, that are dynamically
36// allocated and distributed in clusters, with one sub-buffer per cluster:
37// - data[N] contains N input data points,
38// - trans[N] contains N intermediate data points,
39// - twid[N] contains N coefs : exp(2*pi*i*j/N) / i and j in [0,rootN-1]
40// Each sub-buffer contains (N/nclusters) entries, with 2 double per entry.
41// These distributed buffers are allocated and initialised in parallel
42// by the working threads running on core 0 in each cluster.
43//
44// Each working thread allocates also a private coefs[rootN-1] buffer,
45// that contains all coefs required for a rootN points FFT.
46//
47// There is one working thread per core.
48// The actual number of cores and cluster in a given hardware architecture
49// is obtained by the get_config() syscall (x_size, y_size, ncores).
50// The max number of clusters is bounded by (X_MAX * Y_MAX).
51// The max number of cores per cluster is bounded by CORES_MAX.
52//
53// Several configuration parameters can be defined below:
54//  - PRINT_ARRAY : Print out complex data points arrays.
55//  - CHECK       : Perform both FFT and inverse FFT to check output/input.
56//  - DEBUG_MAIN  : Display intermediate results in main()
57//  - DEBUG_FFT1D : Display intermediate results in FFT1D()
58//  - DEBUG_ROW   : Display intermedite results in FFTrow()
59//
60// Regarding final instrumentation:
61// - the sequencial initialisation time (init_time) is computed
62//   by the main thread in the main() function.
63// - The parallel execution time (parallel_time[i]) is computed by each
64//   working thread(i) in the work() function.
65// - The synchronisation time related to the barriers (sync_time[i])
66//   is computed by each thread(i) in the work() function.
67// The results are displayed on the TXT terminal, and registered on disk.
68///////////////////////////////////////////////////////////////////////////
69
70#include <math.h>
71#include <stdio.h>
72#include <stdlib.h>
73#include <fcntl.h>
74#include <unistd.h>
75#include <pthread.h>
76#include <almosmkh.h>
77#include <hal_macros.h>
78
79// constants
80
81#define PI                      3.14159265359
82#define PAGE_SIZE               4096
83#define X_MAX                   16              // max number of clusters in a row
84#define Y_MAX                   16              // max number of clusters in a column
85#define CORES_MAX               4               // max number of cores in a cluster
86#define CLUSTERS_MAX            X_MAX * Y_MAX
87#define THREADS_MAX             CLUSTERS_MAX * CORES_MAX
88#define RANDOM                  0
89#define COSIN                   1
90#define CONSTANT                2
91
92// parameters
93
94#define DEFAULT_M               14              // 16384 data points
95#define USE_DQT_BARRIER         1               // use DDT barrier if non zero
96#define MODE                    COSIN           // DATA array initialisation mode
97#define CHECK                   0               
98#define DEBUG_MAIN              0               // trace main() function (detailed if odd)
99#define DEBUG_WORK              0               // trace work() function (detailed if odd)
100#define DEBUG_FFT1D             0               // trace FFT1D() function (detailed if odd)
101#define DEBUG_ROW               0               // trace FFTRow() function (detailed if odd)
102#define PRINT_ARRAY             0
103#define DISPLAY_SCHED_AND_VMM   0               // display final VMM state in all clusters
104
105// macro to swap two variables
106#define SWAP(a,b) { double tmp; tmp = a; a = b; b = tmp; }
107
108/////////////////////////////////////////////////////////////////////////////////////
109//             FFT specific global variables
110/////////////////////////////////////////////////////////////////////////////////////
111
112// work function arguments
113typedef struct work_args_s
114{
115    unsigned int        tid;               // thread continuous index
116    unsigned int        lid;               // core local index
117    unsigned int        cid;               // cluster continuous index
118    pthread_barrier_t * parent_barrier;    // parent barrier to signal completion
119}
120work_args_t;
121
122unsigned int   nthreads;                   // total number of threads (one thread per core)
123unsigned int   nclusters;                  // total number of clusters
124unsigned int   M = DEFAULT_M;              // log2(number of points)
125unsigned int   N;                          // number of points (N = 2^M)         
126unsigned int   rootN;                      // rootN = 2^M/2   
127unsigned int   rows_per_thread;            // number of data "rows" handled by a single thread
128unsigned int   points_per_cluster;         // number of data points per cluster
129
130// arrays of pointers on distributed buffers (one sub-buffer per cluster)
131double *       data[CLUSTERS_MAX];         // original time-domain data
132double *       trans[CLUSTERS_MAX];        // used as auxiliary space for transpose
133double *       twid[CLUSTERS_MAX];         // twiddle factor : exp(-2iPI*k*n/N)
134double *       bloup[CLUSTERS_MAX];        // used as auxiliary space for DFT
135
136// instrumentation counters
137unsigned int   parallel_time[THREADS_MAX]; // total computation time (per thread)
138unsigned int   sync_time[THREADS_MAX];     // cumulated waiting time in barriers (per thread)
139unsigned int   init_time;                  // initialisation time (in main)
140
141// synchronisation barrier (all threads)
142pthread_barrier_t      barrier;
143pthread_barrierattr_t  barrier_attr;
144
145/////////////////////////////////////////////////////////////////////////////////////
146//             Global variables required by parallel_pthread_create()
147/////////////////////////////////////////////////////////////////////////////////////
148
149// 2D arrays of input arguments for the <work> threads
150// These arrays are initialised by the application main thread
151
152work_args_t       work_args[CLUSTERS_MAX][CORES_MAX];  // work function arguments
153work_args_t     * work_ptrs[CLUSTERS_MAX][CORES_MAX];  // pointers on arguments
154
155// 1D array of barriers to allow the <work> threads to signal termination
156// this array is initialised in each cluster by the <build[cxy][0]> thread
157 
158pthread_barrier_t parent_barriers[CLUSTERS_MAX];        // termination barrier
159
160/////////////////////////////////////////////////////////////////////////////////////
161//           functions declaration
162/////////////////////////////////////////////////////////////////////////////////////
163
164void work( work_args_t * args );
165
166double CheckSum( void );
167
168void InitD( double    ** data , 
169            unsigned int mode,
170            unsigned int tid );
171
172void InitT( double    ** twid,
173            unsigned int tid );
174
175void InitU( double * coefs );
176
177unsigned int BitReverse( unsigned int k );
178
179void FFT1D( int          direction,
180            double    ** x,
181            double    ** tmp,
182            double     * upriv, 
183            double    ** twid,
184            unsigned int tid,
185            unsigned int MyFirst,
186            unsigned int MyLast );
187
188void TwiddleOneCol( int          direction,
189                    unsigned int j,
190                    double    ** u,
191                    double    ** x,
192                    unsigned int offset_x );
193
194void Scale( double    ** x,
195            unsigned int offset_x );
196
197void Transpose( double    ** src, 
198                double    ** dest,
199                unsigned int MyFirst,
200                unsigned int MyLast );
201
202void Copy( double    ** src,
203           double    ** dest,
204           unsigned int MyFirst,
205           unsigned int MyLast );
206
207void Reverse( double    ** x, 
208              unsigned int offset_x );
209
210void FFTRow( int          direction,
211                double     * u,
212                double    ** x,
213                unsigned int offset_x );
214
215void PrintArray( double ** x,
216                 unsigned int size );
217
218void SimpleDft( int          direction,
219                unsigned int size,
220                double    ** src,
221                unsigned int src_offset,
222                double    ** dst,
223                unsigned int dst_offset );
224
225///////////////////////////////////////////////////////////////////
226// This main() function execute the sequencial initialisation
227// launch the parallel execution, and makes the instrumentation.
228///////////////////////////////////////////////////////////////////
229void main ( void )
230{
231    int                 error;
232
233    unsigned int        x_size;            // number of clusters per row
234    unsigned int        y_size;            // number of clusters per column
235    unsigned int        ncores;            // max number of cores per cluster
236
237
238    unsigned int        x;                 // current index for cluster X coordinate
239    unsigned int        y;                 // current index for cluster Y coordinate
240    unsigned int        lid;               // current index for core in a cluster
241    unsigned int        tid;               // continuous thread index
242    unsigned int        cid;               // cluster continuous index
243    unsigned int        cxy;               // hardware specific cluster identifier
244
245    char                name[64];          // instrumentation file name
246    char                path[128];         // instrumentation path name
247    char                string[256];
248    int                 ret;
249
250    unsigned long long  start_init_cycle; 
251    unsigned long long  end_init_cycle;
252
253#if DEBUG_MAIN
254    unsigned long long  debug_cycle;
255#endif
256
257#if CHECK
258    double              ck1;               // for input/output checking
259    double              ck3;               // for input/output checking
260#endif
261   
262    int                 pid = getpid();
263
264    // get FFT application start cycle
265    get_cycle( &start_init_cycle );
266
267    // get platform parameters
268    if( get_config( &x_size , &y_size , &ncores ) )
269    {
270        printf("\n[fft error] cannot get hardware configuration\n");
271        exit( 0 );
272    }
273
274    // check ncores
275    if( (ncores != 1) && (ncores != 2) && (ncores != 4) )
276    {
277        printf("\n[fft error] number of cores per cluster must be 1/2/4\n");
278        exit( 0 );
279    }
280
281    // check x_size
282    if( (x_size != 1) && (x_size != 2) && (x_size != 4) && (x_size != 8) && (x_size != 16) )
283    {
284        printf("\n[fft error] x_size must be 1/2/4/8/16\n");
285        exit( 0 );
286    }
287
288    // check y_size
289    if( (y_size != 1) && (y_size != 2) && (y_size != 4) && (y_size != 8) && (y_size != 16) )
290    {
291        printf("\n[fft error] y_size must be 1/2/4/8/16\n");
292        exit( 0 );
293    }
294
295    // compute nthreads and nclusters
296    nthreads  = x_size * y_size * ncores;
297    nclusters = x_size * y_size;
298
299    // compute covering DQT size an level
300    unsigned int z = (x_size > y_size) ? x_size : y_size;
301    unsigned int root_level = (z == 1) ? 0 : (z == 2) ? 1 : (z == 4) ? 2 : (z == 8) ? 3 : 4;
302
303    // compute various constants depending on N and T
304    N                  = 1 << M;
305    rootN              = 1 << (M / 2);
306    rows_per_thread    = rootN / nthreads;
307    points_per_cluster = N / nclusters;
308 
309    // check N versus T
310    if( rootN < nthreads )
311    {
312        printf("\n[fft error] sqrt(N) must be larger than T\n");
313        exit( 0 );
314    }
315
316    printf("\n[fft] starts / %d points / %d thread(s) / PID %x / cycle %d\n",
317    N, nthreads, pid, (unsigned int)start_init_cycle );
318
319    // build instrumentation file name
320    if( USE_DQT_BARRIER )
321    snprintf( name , 64 , "p_fft_dqt_%d_%d_%d", N , x_size * y_size , ncores );
322    else
323    snprintf( name , 64 , "p_fft_smp_%d_%d_%d", N , x_size * y_size , ncores );
324
325    // build pathname
326    snprintf( path , 128 , "/home/%s", name );
327
328    // open instrumentation file
329    FILE * f = fopen( path , NULL );
330    if ( f == NULL ) 
331    { 
332        printf("\n[fft error] cannot open instrumentation file <%s>\n", path );
333        exit( 0 );
334    }
335
336#if DEBUG_MAIN
337get_cycle( &debug_cycle );
338printf("\n[fft] main open file <%s> at cycle %d\n",
339path, (unsigned int)debug_cycle );
340#endif
341
342#if CHECK
343ck1 = CheckSum();
344#endif
345
346#if PRINT_ARRAY
347printf("\nData values / base = %x\n", &data[0][0] );
348PrintArray( data , N );
349
350printf("\nTwiddle values / base = %x\n", &twid[0][0] );
351PrintArray( twid , N );
352
353SimpleDft( 1 , N , data , 0 , bloup , 0 );
354
355printf("\nExpected results / base = %x\n", &bloup[0][0] );
356PrintArray( bloup , N );
357#endif
358
359    // initialise barrier synchronizing all <work> threads
360    if( USE_DQT_BARRIER )
361    {
362        barrier_attr.x_size   = x_size;
363        barrier_attr.y_size   = y_size;
364        barrier_attr.nthreads = ncores;
365        error = pthread_barrier_init( &barrier, &barrier_attr , nthreads );
366    }
367    else
368    {
369        error = pthread_barrier_init( &barrier, NULL , nthreads );
370    }
371
372    if( error )
373    {
374        printf("\n[fft error] cannot initialize barrier\n");
375        exit( 0 );
376    }
377
378#if DEBUG_MAIN
379get_cycle( &debug_cycle );
380printf("\n[fft] main completes barrier init at cycle %d\n",
381(unsigned int)debug_cycle );
382#endif
383
384    // build array of arguments for the <work> threads
385    for (x = 0 ; x < x_size ; x++)
386    {
387        for (y = 0 ; y < y_size ; y++)
388        {
389            // compute cluster identifier
390            cxy = HAL_CXY_FROM_XY( x , y );
391
392            for ( lid = 0 ; lid < ncores ; lid++ )
393            {
394                // compute cluster continuous index
395                cid = (x * y_size) + y;
396
397                // compute work thread continuous index
398                tid = (cid * ncores) + lid;
399               
400                // initialize 2D array of arguments
401                work_args[cxy][lid].tid            = tid;
402                work_args[cxy][lid].lid            = lid;
403                work_args[cxy][lid].cid            = cid;
404                work_args[cxy][lid].parent_barrier = &parent_barriers[cxy];
405
406                // initialize 2D array of pointers
407                work_ptrs[cxy][lid] = &work_args[cxy][lid];
408            }
409        }
410    }
411
412    // register sequencial time
413    get_cycle( &end_init_cycle );
414    init_time = (unsigned int)(end_init_cycle - start_init_cycle);
415
416#if DEBUG_MAIN
417printf("\n[fft] main completes <work> threads arguments at cycle %d\n",
418(unsigned int)end_init_cycle );
419#endif
420
421    // create and execute the working threads
422    if( pthread_parallel_create( root_level,
423                                 &work,
424                                 &work_ptrs[0][0],
425                                 &parent_barriers[0] ) )
426    {
427        printf("\n[fft error] creating threads\n");
428        exit( 0 );
429    }
430
431#if DEBUG_MAIN
432get_cycle( &debug_cycle );
433printf("\n[fft] main resume for instrumentation at cycle %d\n",
434(unsigned int)debug_cycle) ;
435#endif
436
437#if PRINT_ARRAY
438printf("\nData values after FFT:\n");
439PrintArray( data , N );
440#endif
441
442#if CHECK
443ck3 = CheckSum();
444printf("\n*** Results ***\n");
445printf("Checksum difference is %f (%f, %f)\n", ck1 - ck3, ck1, ck3);
446if (fabs(ck1 - ck3) < 0.001)  printf("Results OK\n");
447else                          printf("Results KO\n");
448#endif
449
450    // display header on terminal, and save to file
451    printf("\n----- %s -----\n", name );
452
453    ret = fprintf( f , "\n----- %s -----\n", name );
454    if( ret < 0 )
455    {
456        printf("\n[fft error] cannot write header to file <%s>\n", path );
457        exit(0);
458    }
459
460    // get instrumentation results for each thread
461    for (tid = 0 ; tid < nthreads ; tid++) 
462    {
463        snprintf( string , 256 , "- tid %d : Sequencial %d / Parallel %d / Barrier %d\n",
464        tid, init_time, parallel_time[tid], sync_time[tid] );
465
466        // save  to instrumentation file
467        fprintf( f , "%s" , string );
468        if( ret < 0 )
469        {
470            printf("\n[fft error] cannot write thread %d to file <%s>\n", tid, path );
471            printf("%s", string );
472            exit(0);
473        }
474    }
475
476    // compute min/max values
477    unsigned int min_para = parallel_time[0];
478    unsigned int max_para = parallel_time[0];
479    unsigned int min_sync = sync_time[0];
480    unsigned int max_sync = sync_time[0];
481
482    for (tid = 0 ; tid < nthreads ; tid++) 
483    {
484        if (parallel_time[tid] > max_para)  max_para = parallel_time[tid];
485        if (parallel_time[tid] < min_para)  min_para = parallel_time[tid];
486        if (sync_time[tid]     > max_sync)  max_sync = sync_time[tid];
487        if (sync_time[tid]     < min_sync)  min_sync = sync_time[tid];
488    }
489
490    // display MIN/MAX values on terminal and save to file
491    snprintf( string , 256 , "\n      Sequencial  Parallel       Barrier\n"
492                             "MIN : %d\t | %d\t | %d\t   (cycles)\n" 
493                             "MAX : %d\t | %d\t | %d\t   (cycles)\n",
494                             (int)init_time, (int)min_para, (int)min_sync,
495                             (int)init_time, (int)max_para, (int)max_sync );
496    printf("%s", string );
497    ret = fprintf( f , "%s", string );
498    if( ret < 0 )
499    {
500        printf("\n[fft error] cannot write MIN/MAX to file <%s>\n", path );
501        exit(0);
502    }
503
504#if DEBUG_MAIN
505get_cycle( &debug_cycle );
506printf("\n[fft] main close file <%s> at cycle %d\n",
507path, (unsigned int)debug_cycle );
508#endif
509
510    exit( 0 );
511
512} // end main()
513
514/////////////////////////////////////////////////////////////////
515// This function is executed in parallel by all <work> threads.
516/////////////////////////////////////////////////////////////////
517void work( work_args_t * args ) 
518{
519    unsigned int        tid;              // this thread continuous index
520    unsigned int        lid;              // core local index
521    unsigned int        cid;              // cluster continuous index
522    pthread_barrier_t * parent_barrier;   // pointer on parent barrier
523
524    unsigned int        MyFirst;          // index first row allocated to thread
525    unsigned int        MyLast;           // index last row allocated to thread
526    double            * upriv;            // private array of FFT coefs
527
528    unsigned long long  parallel_start;
529    unsigned long long  parallel_stop;
530    unsigned long long  barrier_start;
531    unsigned long long  barrier_stop;
532
533    // get thread arguments
534    tid            = args->tid; 
535    lid            = args->lid;             
536    cid            = args->cid;             
537    parent_barrier = args->parent_barrier;
538
539    get_cycle( &parallel_start );
540
541#if DEBUG_WORK
542printf("\n[fft] %s : thread %d enter / cycle %d\n",
543__FUNCTION__, tid, (unsigned int)parallel_start );
544#endif
545
546    // thread on core 0 allocates memory from the local cluster
547    // for the distributed data[], trans[], twid[] buffers
548    if( lid == 0 )
549    {
550        unsigned int data_size = (N / nclusters) * 2 * sizeof(double);
551
552        data[cid] = (double *)malloc( data_size ); 
553        if( data[cid] == NULL )
554        {
555            printf("\n[fft_error] in work : cannot allocate data[%d] buffer\n", cid );
556            pthread_barrier_wait( parent_barrier );
557            pthread_exit( NULL );
558        }
559       
560        trans[cid] = (double *)malloc( data_size ); 
561        if( trans[cid] == NULL )
562        {
563            printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid );
564            pthread_barrier_wait( parent_barrier );
565            pthread_exit( NULL );
566        }
567       
568        twid[cid] = (double *)malloc( data_size ); 
569        if( twid[cid] == NULL )
570        {
571            printf("\n[fft_error] in work : cannot allocate twid[%d] buffer\n", cid );
572            pthread_barrier_wait( parent_barrier );
573            pthread_exit( NULL );
574        }
575    }
576
577    // BARRIER to wait distributed buffers allocation
578    get_cycle( &barrier_start );
579    pthread_barrier_wait( &barrier );
580    get_cycle( &barrier_stop );
581    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
582
583#if DEBUG_WORK
584printf("\n[fft] %s : thread %d exit barrier for buffer allocation / cycle %d\n",
585__FUNCTION__, tid, (unsigned int)barrier_stop );
586#endif
587
588#if DISPLAY_SCHED_AND_VMM
589    unsigned int x_size;
590    unsigned int y_size;
591    unsigned int ncores;
592    get_config( &x_size , &y_size , &ncores );
593    unsigned int x   = cid / y_size;
594    unsigned int y   = cid % y_size;
595    unsigned int cxy = HAL_CXY_FROM_XY( x , y );
596display_sched( cxy , lid );
597if( lid == 0 ) display_vmm( cxy , getpid() , 0 );
598#endif
599
600    // all threads contribute to data[] local array initialisation
601    InitD( data , MODE , tid ); 
602
603    // all threads contribute to data[] local array initialisation
604    InitT( twid , tid );
605   
606    // BARRIER to wait distributed buffers initialisation
607    get_cycle( &barrier_start );
608    pthread_barrier_wait( &barrier );
609    get_cycle( &barrier_stop );
610    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
611
612#if DEBUG_WORK
613printf("\n[fft] %s : thread %d exit barrier for buffer initialisation / cycle %d\n",
614__FUNCTION__, tid, (unsigned int)barrier_stop );
615#endif
616
617    // all threads allocate memory from the local cluster
618    // for the private upriv[] buffer
619    upriv = (double *)malloc( (rootN - 1) * 2 * sizeof(double) );
620    if( upriv == NULL )
621    {
622        printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid );
623        pthread_barrier_wait( parent_barrier );
624        pthread_exit( NULL );
625    }
626
627    // all threads initialise the private upriv[] array
628    InitU( upriv );
629
630    // all threads compute first and last rows handled by the thread
631    MyFirst = rootN * tid / nthreads;
632    MyLast  = rootN * (tid + 1) / nthreads;
633
634    // all threads perform forward FFT
635    FFT1D( 1 , data , trans , upriv , twid , tid , MyFirst , MyLast );
636
637#if CHECK
638get_cycle( &barrier_start );
639pthread_barrier_wait( &barrier );
640get_cycle( &barrier_stop );
641sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
642FFT1D( -1 , data , trans , upriv , twid , tid , MyFirst , MyLast );
643#endif
644
645    get_cycle( &parallel_stop );
646
647    // register parallel time
648    parallel_time[tid] = (unsigned int)(parallel_stop - parallel_start);
649
650#if DEBUG_WORK
651printf("\n[fft] %s : thread %d completes fft / p_start %d / p_stop %d\n", 
652__FUNCTION__, tid, (unsigned int)parallel_start, (unsigned int)parallel_stop );
653#endif
654
655    //  work thread signals completion to main
656    pthread_barrier_wait( parent_barrier );
657
658#if DEBUG_WORK
659printf("\n[fft] %s : thread %d exit\n", 
660__FUNCTION__, tid );
661#endif
662
663#if DISPLAY_SCHED_AND_VMM
664printf("\n[fft] %s : thread %d exit\n", __FUNCTION__, tid );
665if( lid == 0 ) display_vmm( cxy , getpid() , 0 );
666#endif
667
668    //  work thread exit
669    pthread_exit( NULL );
670
671}  // end work()
672
673////////////////////////////////////////////////////////////////////////////////////////
674// This function makes the DFT from the src[nclusters][points_per_cluster] distributed
675// buffer, to the dst[nclusters][points_per_cluster] distributed buffer.
676////////////////////////////////////////////////////////////////////////////////////////
677void SimpleDft( int             direction,      // 1 direct / -1 reverse
678                unsigned int    size,           // number of points
679                double       ** src,            // source distributed buffer
680                unsigned int    src_offset,     // offset in source array
681                double       ** dst,            // destination distributed buffer
682                unsigned int    dst_offset )    // offset in destination array
683{
684    unsigned int  n , k;
685    double        phi;            // 2*PI*n*k/N
686    double        u_r;            // cos( phi )
687    double        u_c;            // sin( phi )
688    double        d_r;            // Re(data[n])
689    double        d_c;            // Im(data[n])
690    double        accu_r;         // Re(accu)
691    double        accu_c;         // Im(accu)
692    unsigned int  c_id;           // distributed buffer cluster index
693    unsigned int  c_offset;       // offset in distributed buffer
694
695    for ( k = 0 ; k < size ; k++ )       // loop on the output data points
696    {
697        // initialise accu
698        accu_r = 0;
699        accu_c = 0;
700
701        for ( n = 0 ; n < size ; n++ )   // loop on the input data points
702        {
703            // compute coef
704            phi = (double)(2*PI*n*k) / size;
705            u_r =  cos( phi );
706            u_c = -sin( phi ) * direction;
707
708            // get input data point
709            c_id     = (src_offset + n) / (points_per_cluster);
710            c_offset = (src_offset + n) % (points_per_cluster);
711            d_r      = src[c_id][2*c_offset];
712            d_c      = src[c_id][2*c_offset+1];
713
714            // increment accu
715            accu_r += ((u_r*d_r) - (u_c*d_c));
716            accu_c += ((u_r*d_c) + (u_c*d_r));
717        }
718
719        // scale for inverse DFT
720        if ( direction == -1 )
721        {
722            accu_r /= size;
723            accu_c /= size;
724        }
725
726        // set output data point
727        c_id     = (dst_offset + k) / (points_per_cluster);
728        c_offset = (dst_offset + k) % (points_per_cluster);
729        dst[c_id][2*c_offset]   = accu_r;
730        dst[c_id][2*c_offset+1] = accu_c;
731    }
732
733}  // end SimpleDft()
734
735///////////////////////
736double CheckSum( void )
737{
738    unsigned int         i , j;
739    unsigned int         c_id;
740    unsigned int         c_offset;
741    double               cks;
742
743    cks = 0.0;
744    for (j = 0; j < rootN ; j++) 
745    {
746        for (i = 0; i < rootN ; i++) 
747        {
748            c_id      = (rootN * j + i) / (points_per_cluster);
749            c_offset  = (rootN * j + i) % (points_per_cluster);
750
751            cks += data[c_id][2*c_offset] + data[c_id][2*c_offset+1];
752        }
753    }
754    return(cks);
755}
756
757//////////////////////////////////////////////////////////////////////////////////////
758// Each working thread <tid> contributes to initialize (rootN / nthreads) rows,
759// in the shared - and distributed - <data> array.
760//////////////////////////////////////////////////////////////////////////////////////
761void InitD(double      ** data,
762           unsigned int   mode,
763           unsigned int   tid ) 
764{
765    unsigned int    i , j;
766    unsigned int    c_id;
767    unsigned int    c_offset;
768    unsigned int    index;
769
770    // compute row_min and row_max
771    unsigned int    row_min = tid * rows_per_thread;
772    unsigned int    row_max = row_min + rows_per_thread;
773
774    for ( j = row_min ; j < row_max ; j++ )      // loop on rows
775    { 
776        for ( i = 0 ; i < rootN ; i++ )          // loop on points in a row
777        { 
778            index     = j * rootN + i;
779            c_id      = index / (points_per_cluster);
780            c_offset  = index % (points_per_cluster);
781
782            // complex input signal is random
783            if ( mode == RANDOM )               
784            {
785                data[c_id][2*c_offset]   = ( (double)rand() ) / 65536;
786                data[c_id][2*c_offset+1] = ( (double)rand() ) / 65536;
787            }
788           
789
790            // complex input signal is cos(n/N) / sin(n/N)
791            if ( mode == COSIN )               
792            {
793                double phi = (double)( 2 * PI * index) / N;
794                data[c_id][2*c_offset]   = cos( phi );
795                data[c_id][2*c_offset+1] = sin( phi );
796            }
797
798            // complex input signal is constant
799            if ( mode == CONSTANT )               
800            {
801                data[c_id][2*c_offset]   = 1.0;
802                data[c_id][2*c_offset+1] = 0.0;
803            }
804        }
805    }
806}
807
808///////////////////////////////////////////////////////////////////////////////////////
809// Each working thread <tid> contributes to initialize (rootN / nthreads) rows,
810// in the shared - and distributed - <twiddle> array.
811///////////////////////////////////////////////////////////////////////////////////////
812void InitT( double      ** twid,
813            unsigned int   tid )
814{
815    unsigned int    i, j;
816    unsigned int    index;
817    unsigned int    c_id;
818    unsigned int    c_offset;
819    double  phi;
820
821    // compute row_min and row_max
822    unsigned int    row_min = tid * rows_per_thread;
823    unsigned int    row_max = row_min + rows_per_thread;
824
825    for ( j = row_min ; j < row_max ; j++ )      // loop on rows
826    { 
827        for ( i = 0 ; i < rootN ; i++ )          // loop on points in a row
828        { 
829            index     = j * rootN + i;
830            c_id      = index / (points_per_cluster);
831            c_offset  = index % (points_per_cluster);
832
833            phi = (double)(2.0 * PI * i * j) / N;
834            twid[c_id][2*c_offset]   = cos( phi );
835            twid[c_id][2*c_offset+1] = -sin( phi );
836        }
837    }
838}
839
840///////////////////////////////////////////////////////////////////////////////////////
841// Each working thread initialize the private <upriv> array / (rootN - 1) entries.
842///////////////////////////////////////////////////////////////////////////////////////
843void InitU( double * upriv ) 
844{
845    unsigned int    q; 
846    unsigned int    j; 
847    unsigned int    base; 
848    unsigned int    n1;
849    double          phi;
850
851    for (q = 0 ; ((unsigned int)(1 << q) < N) ; q++) 
852    { 
853        n1 = 1 << q;    // n1 == 2**q
854        base = n1 - 1;
855        for (j = 0; (j < n1) ; j++) 
856        {
857            if (base + j > rootN - 1) return;
858
859            phi = (double)(2.0 * PI * j) / (2 * n1);
860            upriv[2*(base+j)]   = cos( phi );
861            upriv[2*(base+j)+1] = -sin( phi );
862        }
863    }
864}
865
866////////////////////////////////////////////////////////////////////////////////////////
867// This function returns an index value that is the bit reverse of the input value.
868////////////////////////////////////////////////////////////////////////////////////////
869unsigned int BitReverse( unsigned int k ) 
870{
871    unsigned int i; 
872    unsigned int j; 
873    unsigned int tmp;
874
875    j = 0;
876    tmp = k;
877    for (i = 0; i < M/2 ; i++) 
878    {
879        j = 2 * j + (tmp & 0x1);
880        tmp = tmp >> 1;
881    }
882    return j;
883}
884
885////////////////////////////////////////////////////////////////////////////////////////
886// This function perform the in place (direct or inverse) FFT on the N data points
887// contained in the distributed buffers x[nclusters][points_per_cluster].
888// It handles the (N) points 1D array as a (rootN*rootN) points 2D array.
889// 1) it transpose (rootN/nthreads ) rows from x to tmp.
890// 2) it make (rootN/nthreads) FFT on the tmp rows and apply the twiddle factor.
891// 3) it transpose (rootN/nthreads) columns from tmp to x.
892// 4) it make (rootN/nthreads) FFT on the x rows.
893// It calls the FFTRow() 2*(rootN/nthreads) times to perform the in place FFT
894// on the rootN points contained in a row.
895////////////////////////////////////////////////////////////////////////////////////////
896void FFT1D( int              direction,       // direct 1 / inverse -1
897            double       **  x,               // input & output distributed data points array
898            double       **  tmp,             // auxiliary distributed data points array
899            double        *  upriv,           // local array containing coefs for rootN FFT
900            double       **  twid,            // distributed arrays containing N twiddle factors
901            unsigned int     tid,             // thread continuous index
902            unsigned int     MyFirst, 
903            unsigned int     MyLast )
904{
905    unsigned int j;
906    unsigned long long barrier_start;
907    unsigned long long barrier_stop;
908
909#if DEBUG_FFT1D
910unsigned long long cycle;
911get_cycle( &cycle );
912printf("\n[fft] %s : thread %d enter / first %d / last %d / cycle %d\n",
913__FUNCTION__, tid, MyFirst, MyLast, (unsigned int)cycle );
914#endif
915
916    // transpose (rootN/nthreads) rows from x to tmp
917    Transpose( x , tmp , MyFirst , MyLast );
918
919#if( DEBUG_FFT1D & 1 )
920get_cycle( &cycle );
921printf("\n[fft] %s : thread %d after first transpose / cycle %d\n",
922__FUNCTION__, tid, (unsigned int)cycle );
923if( PRINT_ARRAY ) PrintArray( tmp , N );
924#endif
925
926    // BARRIER
927    get_cycle( &barrier_start );
928    pthread_barrier_wait( &barrier );
929    get_cycle( &barrier_stop );
930    sync_time[tid] = (unsigned int)(barrier_stop - barrier_start);
931
932#if( DEBUG_FFT1D & 1 )
933get_cycle( &cycle );
934printf("\n[fft] %s : thread %d exit barrier after first transpose / cycle %d\n",
935__FUNCTION__, tid, (unsigned int)cycle );
936#endif
937
938    // do FFTs on rows of tmp (i.e. columns of x) and apply twiddle factor
939    for (j = MyFirst; j < MyLast; j++) 
940    {
941        FFTRow( direction , upriv , tmp , j * rootN );
942
943        TwiddleOneCol( direction , j , twid , tmp , j * rootN );
944    } 
945
946#if( DEBUG_FFT1D & 1 )
947printf("\n[fft] %s : thread %d after first twiddle\n", __FUNCTION__, tid);
948if( PRINT_ARRAY ) PrintArray( tmp , N );
949#endif
950
951    // BARRIER
952    get_cycle( &barrier_start );
953    pthread_barrier_wait( &barrier );
954    get_cycle( &barrier_stop );
955
956#if( DEBUG_FFT1D & 1 )
957printf("\n[fft] %s : thread %d exit barrier after first twiddle\n", __FUNCTION__, tid);
958#endif
959
960    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
961
962    // transpose tmp to x
963    Transpose( tmp , x , MyFirst , MyLast );
964
965#if( DEBUG_FFT1D & 1 )
966printf("\n[fft] %s : thread %d after second transpose\n", __FUNCTION__, tid);
967if( PRINT_ARRAY ) PrintArray( x , N );
968#endif
969
970    // BARRIER
971    get_cycle( &barrier_start );
972    pthread_barrier_wait( &barrier );
973    get_cycle( &barrier_stop );
974
975#if( DEBUG_FFT1D & 1 )
976printf("\n[fft] %s : thread %d exit barrier after second transpose\n", __FUNCTION__, tid);
977#endif
978
979    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
980
981    // do FFTs on rows of x and apply the scaling factor
982    for (j = MyFirst; j < MyLast; j++) 
983    {
984        FFTRow( direction , upriv , x , j * rootN );
985        if (direction == -1) Scale( x , j * rootN );
986    }
987
988#if( DEBUG_FFT1D & 1 )
989printf("\n[fft] %s : thread %d after FFT on rows\n", __FUNCTION__, tid);
990if( PRINT_ARRAY ) PrintArray( x , N );
991#endif
992
993    // BARRIER
994    get_cycle( &barrier_start );
995    pthread_barrier_wait( &barrier );
996    get_cycle( &barrier_stop );
997
998#if( DEBUG_FFT1D & 1 )
999printf("\n[fft] %s : thread %d exit barrier after FFT on rows\n", __FUNCTION__, tid);
1000#endif
1001    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1002
1003    // transpose x to tmp
1004    Transpose( x , tmp , MyFirst , MyLast );
1005
1006#if( DEBUG_FFT1D & 1 )
1007printf("\n[fft] %s : thread %x after third transpose\n", __FUNCTION__, tid);
1008if( PRINT_ARRAY ) PrintArray( x , N );
1009#endif
1010
1011    // BARRIER
1012    get_cycle( &barrier_start );
1013    pthread_barrier_wait( &barrier );
1014    get_cycle( &barrier_stop );
1015
1016#if( DEBUG_FFT1D & 1 )
1017printf("\n[fft] %s : thread %d exit barrier after third transpose\n", __FUNCTION__, tid);
1018#endif
1019
1020    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1021    sync_time[tid] += (long)(barrier_stop - barrier_start);
1022
1023    // copy tmp to x
1024    Copy( tmp , x , MyFirst , MyLast );
1025
1026#if DEBUG_FFT1D
1027printf("\n[fft] %s : thread %d completed\n", __FUNCTION__, tid);
1028if( PRINT_ARRAY ) PrintArray( x , N );
1029#endif
1030
1031}  // end FFT1D()
1032
1033/////////////////////////////////////////////////////////////////////////////////////
1034// This function multiply all points contained in a row (rootN points) of the
1035// x[] array by the corresponding twiddle factor, contained in the u[] array.
1036/////////////////////////////////////////////////////////////////////////////////////
1037void TwiddleOneCol( int             direction, 
1038                    unsigned int    j,              // y coordinate in 2D view of coef array
1039                    double       ** u,              // coef array base address
1040                    double       ** x,              // data array base address
1041                    unsigned int    offset_x )      // first point in N points data array
1042{
1043    unsigned int i;
1044    double       omega_r; 
1045    double       omega_c; 
1046    double       x_r; 
1047    double       x_c;
1048    unsigned int c_id;
1049    unsigned int c_offset;
1050
1051    for (i = 0; i < rootN ; i++)  // loop on the rootN points
1052    {
1053        // get coef
1054        c_id      = (j * rootN + i) / (points_per_cluster);
1055        c_offset  = (j * rootN + i) % (points_per_cluster);
1056        omega_r = u[c_id][2*c_offset];
1057        omega_c = direction * u[c_id][2*c_offset+1];
1058
1059        // access data
1060        c_id      = (offset_x + i) / (points_per_cluster);
1061        c_offset  = (offset_x + i) % (points_per_cluster);   
1062        x_r = x[c_id][2*c_offset]; 
1063        x_c = x[c_id][2*c_offset+1];
1064
1065        x[c_id][2*c_offset]   = omega_r*x_r - omega_c * x_c;
1066        x[c_id][2*c_offset+1] = omega_r*x_c + omega_c * x_r;
1067    }
1068}  // end TwiddleOneCol()
1069
1070////////////////////////////
1071void Scale( double      ** x,           // data array base address
1072            unsigned int   offset_x )   // first point of the row to be scaled
1073{
1074    unsigned int i;
1075    unsigned int c_id;
1076    unsigned int c_offset;
1077
1078    for (i = 0; i < rootN ; i++) 
1079    {
1080        c_id      = (offset_x + i) / (points_per_cluster);
1081        c_offset  = (offset_x + i) % (points_per_cluster);
1082        x[c_id][2*c_offset]     /= N;
1083        x[c_id][2*c_offset + 1] /= N;
1084    }
1085}
1086
1087///////////////////////////////////
1088void Transpose( double      ** src,      // source buffer (array of pointers)
1089                double      ** dest,     // destination buffer (array of pointers)
1090                unsigned int   MyFirst,  // first row allocated to the thread
1091                unsigned int   MyLast )  // last row allocated to the thread
1092{
1093    unsigned int row;               // row index
1094    unsigned int point;             // data point index in a row
1095
1096    unsigned int index_src;         // absolute index in the source N points array
1097    unsigned int c_id_src;          // cluster for the source buffer
1098    unsigned int c_offset_src;      // offset in the source buffer
1099
1100    unsigned int index_dst;         // absolute index in the dest N points array
1101    unsigned int c_id_dst;          // cluster for the dest buffer
1102    unsigned int c_offset_dst;      // offset in the dest buffer
1103
1104   
1105    // scan all data points allocated to the thread
1106    // (between MyFirst row and MyLast row) from the source buffer
1107    // and write these points to the destination buffer
1108    for ( row = MyFirst ; row < MyLast ; row++ )       // loop on the rows
1109    {
1110        for ( point = 0 ; point < rootN ; point++ )    // loop on points in row
1111        {
1112            index_src    = row * rootN + point;
1113            c_id_src     = index_src / (points_per_cluster);
1114            c_offset_src = index_src % (points_per_cluster);
1115
1116            index_dst    = point * rootN + row;
1117            c_id_dst     = index_dst / (points_per_cluster);
1118            c_offset_dst = index_dst % (points_per_cluster);
1119
1120            dest[c_id_dst][2*c_offset_dst]   = src[c_id_src][2*c_offset_src];
1121            dest[c_id_dst][2*c_offset_dst+1] = src[c_id_src][2*c_offset_src+1];
1122        }
1123    }
1124}  // end Transpose()
1125
1126//////////////////////////////
1127void Copy( double      ** src,      // source buffer (array of pointers)
1128           double      ** dest,     // destination buffer (array of pointers)
1129           unsigned int   MyFirst,  // first row allocated to the thread
1130           unsigned int   MyLast )  // last row allocated to the thread
1131{
1132    unsigned int row;                  // row index
1133    unsigned int point;                // data point index in a row
1134
1135    unsigned int index;                // absolute index in the N points array
1136    unsigned int c_id;                 // cluster index
1137    unsigned int c_offset;             // offset in local buffer
1138
1139    // scan all data points allocated to the thread
1140    for ( row = MyFirst ; row < MyLast ; row++ )       // loop on the rows
1141    {
1142        for ( point = 0 ; point < rootN ; point++ )    // loop on points in row
1143        {
1144            index    = row * rootN + point;
1145            c_id     = index / (points_per_cluster);
1146            c_offset = index % (points_per_cluster);
1147
1148            dest[c_id][2*c_offset]   = src[c_id][2*c_offset];
1149            dest[c_id][2*c_offset+1] = src[c_id][2*c_offset+1];
1150        }
1151    }
1152}  // end Copy()
1153
1154///////////////////////////////
1155void Reverse( double      ** x, 
1156              unsigned int   offset_x )
1157{
1158    unsigned int j, k;
1159    unsigned int c_id_j;
1160    unsigned int c_offset_j;
1161    unsigned int c_id_k;
1162    unsigned int c_offset_k;
1163
1164    for (k = 0 ; k < rootN ; k++) 
1165    {
1166        j = BitReverse( k );
1167        if (j > k) 
1168        {
1169            c_id_j      = (offset_x + j) / (points_per_cluster);
1170            c_offset_j  = (offset_x + j) % (points_per_cluster);
1171            c_id_k      = (offset_x + k) / (points_per_cluster);
1172            c_offset_k  = (offset_x + k) % (points_per_cluster);
1173
1174            SWAP(x[c_id_j][2*c_offset_j]  , x[c_id_k][2*c_offset_k]);
1175            SWAP(x[c_id_j][2*c_offset_j+1], x[c_id_k][2*c_offset_k+1]);
1176        }
1177    }
1178}
1179
1180/////////////////////////////////////////////////////////////////////////////
1181// This function makes the in-place FFT on all points contained in a row
1182// (i.e. rootN points) of the x[nclusters][points_per_cluster] array.
1183/////////////////////////////////////////////////////////////////////////////
1184void FFTRow( int            direction,  // 1 direct / -1 inverse
1185                double       * u,          // private coefs array
1186                double      ** x,          // array of pointers on distributed buffers
1187                unsigned int   offset_x )  // absolute offset in the x array
1188{
1189    unsigned int     j;
1190    unsigned int     k;
1191    unsigned int     q;
1192    unsigned int     L;
1193    unsigned int     r;
1194    unsigned int     Lstar;
1195    double * u1; 
1196
1197    unsigned int     offset_x1;     // index first butterfly input
1198    unsigned int     offset_x2;     // index second butterfly output
1199
1200    double           omega_r;       // real part butterfy coef
1201    double           omega_c;       // complex part butterfly coef
1202
1203    double           tau_r;
1204    double           tau_c;
1205
1206    double           d1_r;          // real part first butterfly input
1207    double           d1_c;          // imag part first butterfly input
1208    double           d2_r;          // real part second butterfly input
1209    double           d2_c;          // imag part second butterfly input
1210
1211    unsigned int     c_id_1;        // cluster index for first butterfly input
1212    unsigned int     c_offset_1;    // offset for first butterfly input
1213    unsigned int     c_id_2;        // cluster index for second butterfly input
1214    unsigned int     c_offset_2;    // offset for second butterfly input
1215
1216#if DEBUG_ROW
1217unsigned int p;
1218printf("\n[fft] ROW data in / %d points / offset = %d\n", rootN , offset_x );
1219
1220for ( p = 0 ; p < rootN ; p++ )
1221{
1222    unsigned int index    = offset_x + p;
1223    unsigned int c_id     = index / (points_per_cluster);
1224    unsigned int c_offset = index % (points_per_cluster);
1225    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1226}
1227printf("\n");
1228#endif
1229
1230    // This makes the rootN input points reordering
1231    Reverse( x , offset_x ); 
1232
1233#if DEBUG_ROW
1234printf("\n[fft] ROW data after reverse / %d points / offset = %d\n", rootN , offset_x );
1235
1236for ( p = 0 ; p < rootN ; p++ )
1237{
1238    unsigned int index    = offset_x + p;
1239    unsigned int c_id     = index / (points_per_cluster);
1240    unsigned int c_offset = index % (points_per_cluster);
1241    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1242}
1243printf("\n");
1244#endif
1245
1246    // This implements the multi-stages, in place Butterfly network
1247    for (q = 1; q <= M/2 ; q++)     // loop on stages
1248    {
1249        L = 1 << q;       // number of points per subset for current stage
1250        r = rootN / L;    // number of subsets
1251        Lstar = L / 2;
1252        u1 = &u[2 * (Lstar - 1)];
1253        for (k = 0; k < r; k++)     // loop on the subsets
1254        {
1255            offset_x1  = offset_x + (k * L);            // index first point
1256            offset_x2  = offset_x + (k * L + Lstar);    // index second point
1257
1258#if (DEBUG_ROW & 1)
1259printf("\n ### q = %d / k = %d / x1 = %d / x2 = %d\n", q , k , offset_x1 , offset_x2 );
1260#endif
1261            // makes all in-place butterfly(s) for subset
1262            for (j = 0; j < Lstar; j++) 
1263            {
1264                // get coef
1265                omega_r = u1[2*j];
1266                omega_c = direction * u1[2*j+1];
1267
1268                // get d[x1] address and value
1269                c_id_1      = (offset_x1 + j) / (points_per_cluster);
1270                c_offset_1  = (offset_x1 + j) % (points_per_cluster);
1271                d1_r        = x[c_id_1][2*c_offset_1];
1272                d1_c        = x[c_id_1][2*c_offset_1+1];
1273
1274                // get d[x2] address and value
1275                c_id_2      = (offset_x2 + j) / (points_per_cluster);
1276                c_offset_2  = (offset_x2 + j) % (points_per_cluster);
1277                d2_r        = x[c_id_2][2*c_offset_2];
1278                d2_c        = x[c_id_2][2*c_offset_2+1];
1279
1280#if (DEBUG_ROW & 1)
1281printf("\n ### d1_in = (%f , %f) / d2_in = (%f , %f) / coef = (%f , %f)\n", 
1282                d1_r , d1_c , d2_r , d2_c , omega_r , omega_c);
1283#endif
1284                // tau = omega * d[x2]
1285                tau_r = omega_r * d2_r - omega_c * d2_c;
1286                tau_c = omega_r * d2_c + omega_c * d2_r;
1287
1288                // set new value for d[x1] = d[x1] + omega * d[x2]
1289                x[c_id_1][2*c_offset_1]   = d1_r + tau_r;
1290                x[c_id_1][2*c_offset_1+1] = d1_c + tau_c;
1291
1292                // set new value for d[x2] = d[x1] - omega * d[x2]
1293                x[c_id_2][2*c_offset_2]   = d1_r - tau_r;
1294                x[c_id_2][2*c_offset_2+1] = d1_c - tau_c;
1295
1296#if (DEBUG_ROW & 1)
1297printf("\n ### d1_out = (%f , %f) / d2_out = (%f , %f)\n", 
1298                d1_r + tau_r , d1_c + tau_c , d2_r - tau_r , d2_c - tau_c );
1299#endif
1300            }
1301        }
1302    }
1303
1304#if DEBUG_ROW
1305printf("\n[fft] ROW data out / %d points / offset = %d\n", rootN , offset_x );
1306for ( p = 0 ; p < rootN ; p++ )
1307{
1308    unsigned int index    = offset_x + p;
1309    unsigned int c_id     = index / (points_per_cluster);
1310    unsigned int c_offset = index % (points_per_cluster);
1311    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1312}
1313printf("\n");
1314#endif
1315
1316}  // end FFTRow()
1317
1318///////////////////////////////////////
1319void PrintArray( double       ** array,
1320                 unsigned int    size ) 
1321{
1322    unsigned int  i;
1323    unsigned int  c_id;
1324    unsigned int  c_offset;
1325
1326    // float display
1327    for (i = 0; i < size ; i++) 
1328    {
1329        c_id      = i / (points_per_cluster);
1330        c_offset  = i % (points_per_cluster);
1331
1332        printf(" %f  %f |", array[c_id][2*c_offset], array[c_id][2*c_offset+1]);
1333
1334        if ( (i+1) % 4 == 0)  printf("\n");
1335    }
1336    printf("\n");
1337}
1338
1339
1340// Local Variables:
1341// tab-width: 4
1342// c-basic-offset: 4
1343// c-file-offsets:((innamespace . 0)(inline-open . 0))
1344// indent-tabs-mode: nil
1345// End:
1346
1347// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4
1348
Note: See TracBrowser for help on using the repository browser.