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

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

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