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

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

update build_disk command in Makefile
File size: 48.5 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
491        pgfaults_nr   += pgfault_nr[tid];
492        pgfaults_cost += pgfault_cost[tid];
493
494        if (pgfault_max[tid]   > pgfaults_max) pgfaults_max   = pgfault_max[tid];
495    }
496
497    // display global values on terminal and save to file
498    snprintf( string , 256 ,
499    "\nSeq %d / Para %d / Sync %d / Pgfaults %d ( cost %d / max %d )\n",
500    init_time, time_para, time_sync, pgfaults_nr, (pgfaults_cost / pgfaults_nr), pgfaults_max );
501
502    printf("%s", string );
503
504    // save global values to file
505    ret = fprintf( f , "%s", string );
506
507    if( ret < 0 )
508    {
509        printf("\n[fft error] cannot save global results to file <%s>\n", path );
510        exit(0);
511    }
512
513    // close instrumentation file
514    ret = fclose( f );
515
516    if( ret < 0 )
517    {
518        printf("\n[fft error] cannot close file <%s>\n", path );
519        exit(0);
520    }
521
522#if DEBUG_MAIN
523get_cycle( &debug_cycle );
524printf("\n[fft] main exit <%s> at cycle %d\n",
525path, (unsigned int)debug_cycle );
526#endif
527
528    exit( 0 );
529
530} // end main()
531
532/////////////////////////////////////////////////////////////////
533// This function is executed in parallel by all <work> threads.
534/////////////////////////////////////////////////////////////////
535void work( work_args_t * args ) 
536{
537    unsigned int        tid;              // this thread continuous index
538    unsigned int        lid;              // core local index
539    unsigned int        cid;              // cluster continuous index
540    pthread_barrier_t * parent_barrier;   // pointer on parent barrier
541
542    unsigned int        MyFirst;          // index first row allocated to thread
543    unsigned int        MyLast;           // index last row allocated to thread
544    double            * upriv;            // private array of FFT coefs
545
546    unsigned long long  parallel_start;
547    unsigned long long  parallel_stop;
548    unsigned long long  barrier_start;
549    unsigned long long  barrier_stop;
550
551    // get thread arguments
552    tid            = args->tid; 
553    lid            = args->lid;             
554    cid            = args->cid;             
555    parent_barrier = args->parent_barrier;
556
557    get_cycle( &parallel_start );
558
559#if DEBUG_WORK
560printf("\n[fft] %s : thread %d enter / cycle %d\n",
561__FUNCTION__, tid, (unsigned int)parallel_start );
562#endif
563
564    // thread on core 0 allocates memory from the local cluster
565    // for the distributed data[], trans[], twid[] buffers
566    if( lid == 0 )
567    {
568        unsigned int data_size = (N / nclusters) * 2 * sizeof(double);
569
570        data[cid] = (double *)malloc( data_size ); 
571        if( data[cid] == NULL )
572        {
573            printf("\n[fft_error] in work : cannot allocate data[%d] buffer\n", cid );
574            pthread_barrier_wait( parent_barrier );
575            pthread_exit( NULL );
576        }
577       
578        trans[cid] = (double *)malloc( data_size ); 
579        if( trans[cid] == NULL )
580        {
581            printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid );
582            pthread_barrier_wait( parent_barrier );
583            pthread_exit( NULL );
584        }
585       
586        twid[cid] = (double *)malloc( data_size ); 
587        if( twid[cid] == NULL )
588        {
589            printf("\n[fft_error] in work : cannot allocate twid[%d] buffer\n", cid );
590            pthread_barrier_wait( parent_barrier );
591            pthread_exit( NULL );
592        }
593    }
594
595    // BARRIER to wait distributed buffers allocation
596    get_cycle( &barrier_start );
597    pthread_barrier_wait( &barrier );
598    get_cycle( &barrier_stop );
599    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
600
601#if DEBUG_WORK
602printf("\n[fft] %s : thread %d exit barrier for buffer allocation / cycle %d\n",
603__FUNCTION__, tid, (unsigned int)barrier_stop );
604#endif
605
606#if DISPLAY_SCHED_AND_VMM
607    unsigned int x_size;
608    unsigned int y_size;
609    unsigned int ncores;
610    get_config( &x_size , &y_size , &ncores );
611    unsigned int x   = cid / y_size;
612    unsigned int y   = cid % y_size;
613    unsigned int cxy = HAL_CXY_FROM_XY( x , y );
614display_sched( cxy , lid );
615if( lid == 0 ) display_vmm( cxy , getpid() , 0 );
616#endif
617
618    // all threads contribute to data[] local array initialisation
619    InitD( data , MODE , tid ); 
620
621    // all threads contribute to data[] local array initialisation
622    InitT( twid , tid );
623   
624    // BARRIER to wait distributed buffers initialisation
625    get_cycle( &barrier_start );
626    pthread_barrier_wait( &barrier );
627    get_cycle( &barrier_stop );
628    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
629
630#if DEBUG_WORK
631printf("\n[fft] %s : thread %d exit barrier for buffer initialisation / cycle %d\n",
632__FUNCTION__, tid, (unsigned int)barrier_stop );
633#endif
634
635    // all threads allocate memory from the local cluster
636    // for the private upriv[] buffer
637    upriv = (double *)malloc( (rootN - 1) * 2 * sizeof(double) );
638    if( upriv == NULL )
639    {
640        printf("\n[fft_error] in work : cannot allocate trans[%d] buffer\n", cid );
641        pthread_barrier_wait( parent_barrier );
642        pthread_exit( NULL );
643    }
644
645    // all threads initialise the private upriv[] array
646    InitU( upriv );
647
648    // all threads compute first and last rows handled by the thread
649    MyFirst = rootN * tid / nthreads;
650    MyLast  = rootN * (tid + 1) / nthreads;
651
652    // all threads perform forward FFT
653    FFT1D( 1 , data , trans , upriv , twid , tid , MyFirst , MyLast );
654
655#if CHECK
656get_cycle( &barrier_start );
657pthread_barrier_wait( &barrier );
658get_cycle( &barrier_stop );
659sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
660FFT1D( -1 , data , trans , upriv , twid , tid , MyFirst , MyLast );
661#endif
662
663    get_cycle( &parallel_stop );
664
665    // register parallel time in instrumentation counters
666    parallel_time[tid] = (unsigned int)(parallel_stop - parallel_start);
667
668    // get work thread info for page faults
669    thread_info_t info;
670    get_thread_info( &info );
671   
672    // register page faults in instrumentation counters
673    pgfault_nr[tid]   = info.false_pgfault_nr + 
674                        info.local_pgfault_nr + 
675                        info.global_pgfault_nr;
676    pgfault_cost[tid] = info.false_pgfault_cost + 
677                        info.local_pgfault_cost + 
678                        info.global_pgfault_cost;
679    pgfault_max[tid]  = info.false_pgfault_max + 
680                        info.local_pgfault_max + 
681                        info.global_pgfault_max;
682#if DEBUG_WORK
683printf("\n[fft] %s : thread %d completes fft / p_start %d / p_stop %d\n", 
684__FUNCTION__, tid, (unsigned int)parallel_start, (unsigned int)parallel_stop );
685#endif
686
687    //  work thread signals completion to main
688    pthread_barrier_wait( parent_barrier );
689
690#if DEBUG_WORK
691printf("\n[fft] %s : thread %d exit\n", 
692__FUNCTION__, tid );
693#endif
694
695#if DISPLAY_SCHED_AND_VMM
696printf("\n[fft] %s : thread %d exit\n", __FUNCTION__, tid );
697if( lid == 0 ) display_vmm( cxy , getpid() , 0 );
698#endif
699
700    //  work thread exit
701    pthread_exit( NULL );
702
703}  // end work()
704
705////////////////////////////////////////////////////////////////////////////////////////
706// This function makes the DFT from the src[nclusters][points_per_cluster] distributed
707// buffer, to the dst[nclusters][points_per_cluster] distributed buffer.
708////////////////////////////////////////////////////////////////////////////////////////
709void SimpleDft( int             direction,      // 1 direct / -1 reverse
710                unsigned int    size,           // number of points
711                double       ** src,            // source distributed buffer
712                unsigned int    src_offset,     // offset in source array
713                double       ** dst,            // destination distributed buffer
714                unsigned int    dst_offset )    // offset in destination array
715{
716    unsigned int  n , k;
717    double        phi;            // 2*PI*n*k/N
718    double        u_r;            // cos( phi )
719    double        u_c;            // sin( phi )
720    double        d_r;            // Re(data[n])
721    double        d_c;            // Im(data[n])
722    double        accu_r;         // Re(accu)
723    double        accu_c;         // Im(accu)
724    unsigned int  c_id;           // distributed buffer cluster index
725    unsigned int  c_offset;       // offset in distributed buffer
726
727    for ( k = 0 ; k < size ; k++ )       // loop on the output data points
728    {
729        // initialise accu
730        accu_r = 0;
731        accu_c = 0;
732
733        for ( n = 0 ; n < size ; n++ )   // loop on the input data points
734        {
735            // compute coef
736            phi = (double)(2*PI*n*k) / size;
737            u_r =  cos( phi );
738            u_c = -sin( phi ) * direction;
739
740            // get input data point
741            c_id     = (src_offset + n) / (points_per_cluster);
742            c_offset = (src_offset + n) % (points_per_cluster);
743            d_r      = src[c_id][2*c_offset];
744            d_c      = src[c_id][2*c_offset+1];
745
746            // increment accu
747            accu_r += ((u_r*d_r) - (u_c*d_c));
748            accu_c += ((u_r*d_c) + (u_c*d_r));
749        }
750
751        // scale for inverse DFT
752        if ( direction == -1 )
753        {
754            accu_r /= size;
755            accu_c /= size;
756        }
757
758        // set output data point
759        c_id     = (dst_offset + k) / (points_per_cluster);
760        c_offset = (dst_offset + k) % (points_per_cluster);
761        dst[c_id][2*c_offset]   = accu_r;
762        dst[c_id][2*c_offset+1] = accu_c;
763    }
764
765}  // end SimpleDft()
766
767///////////////////////
768double CheckSum( void )
769{
770    unsigned int         i , j;
771    unsigned int         c_id;
772    unsigned int         c_offset;
773    double               cks;
774
775    cks = 0.0;
776    for (j = 0; j < rootN ; j++) 
777    {
778        for (i = 0; i < rootN ; i++) 
779        {
780            c_id      = (rootN * j + i) / (points_per_cluster);
781            c_offset  = (rootN * j + i) % (points_per_cluster);
782
783            cks += data[c_id][2*c_offset] + data[c_id][2*c_offset+1];
784        }
785    }
786    return(cks);
787}
788
789//////////////////////////////////////////////////////////////////////////////////////
790// Each working thread <tid> contributes to initialize (rootN / nthreads) rows,
791// in the shared - and distributed - <data> array.
792//////////////////////////////////////////////////////////////////////////////////////
793void InitD(double      ** data,
794           unsigned int   mode,
795           unsigned int   tid ) 
796{
797    unsigned int    i , j;
798    unsigned int    c_id;
799    unsigned int    c_offset;
800    unsigned int    index;
801
802    // compute row_min and row_max
803    unsigned int    row_min = tid * rows_per_thread;
804    unsigned int    row_max = row_min + rows_per_thread;
805
806    for ( j = row_min ; j < row_max ; j++ )      // loop on rows
807    { 
808        for ( i = 0 ; i < rootN ; i++ )          // loop on points in a row
809        { 
810            index     = j * rootN + i;
811            c_id      = index / (points_per_cluster);
812            c_offset  = index % (points_per_cluster);
813
814            // complex input signal is random
815            if ( mode == RANDOM )               
816            {
817                data[c_id][2*c_offset]   = ( (double)rand() ) / 65536;
818                data[c_id][2*c_offset+1] = ( (double)rand() ) / 65536;
819            }
820           
821
822            // complex input signal is cos(n/N) / sin(n/N)
823            if ( mode == COSIN )               
824            {
825                double phi = (double)( 2 * PI * index) / N;
826                data[c_id][2*c_offset]   = cos( phi );
827                data[c_id][2*c_offset+1] = sin( phi );
828            }
829
830            // complex input signal is constant
831            if ( mode == CONSTANT )               
832            {
833                data[c_id][2*c_offset]   = 1.0;
834                data[c_id][2*c_offset+1] = 0.0;
835            }
836        }
837    }
838}
839
840///////////////////////////////////////////////////////////////////////////////////////
841// Each working thread <tid> contributes to initialize (rootN / nthreads) rows,
842// in the shared - and distributed - <twiddle> array.
843///////////////////////////////////////////////////////////////////////////////////////
844void InitT( double      ** twid,
845            unsigned int   tid )
846{
847    unsigned int    i, j;
848    unsigned int    index;
849    unsigned int    c_id;
850    unsigned int    c_offset;
851    double  phi;
852
853    // compute row_min and row_max
854    unsigned int    row_min = tid * rows_per_thread;
855    unsigned int    row_max = row_min + rows_per_thread;
856
857    for ( j = row_min ; j < row_max ; j++ )      // loop on rows
858    { 
859        for ( i = 0 ; i < rootN ; i++ )          // loop on points in a row
860        { 
861            index     = j * rootN + i;
862            c_id      = index / (points_per_cluster);
863            c_offset  = index % (points_per_cluster);
864
865            phi = (double)(2.0 * PI * i * j) / N;
866            twid[c_id][2*c_offset]   = cos( phi );
867            twid[c_id][2*c_offset+1] = -sin( phi );
868        }
869    }
870}
871
872///////////////////////////////////////////////////////////////////////////////////////
873// Each working thread initialize the private <upriv> array / (rootN - 1) entries.
874///////////////////////////////////////////////////////////////////////////////////////
875void InitU( double * upriv ) 
876{
877    unsigned int    q; 
878    unsigned int    j; 
879    unsigned int    base; 
880    unsigned int    n1;
881    double          phi;
882
883    for (q = 0 ; ((unsigned int)(1 << q) < N) ; q++) 
884    { 
885        n1 = 1 << q;    // n1 == 2**q
886        base = n1 - 1;
887        for (j = 0; (j < n1) ; j++) 
888        {
889            if (base + j > rootN - 1) return;
890
891            phi = (double)(2.0 * PI * j) / (2 * n1);
892            upriv[2*(base+j)]   = cos( phi );
893            upriv[2*(base+j)+1] = -sin( phi );
894        }
895    }
896}
897
898////////////////////////////////////////////////////////////////////////////////////////
899// This function returns an index value that is the bit reverse of the input value.
900////////////////////////////////////////////////////////////////////////////////////////
901unsigned int BitReverse( unsigned int k ) 
902{
903    unsigned int i; 
904    unsigned int j; 
905    unsigned int tmp;
906
907    j = 0;
908    tmp = k;
909    for (i = 0; i < M/2 ; i++) 
910    {
911        j = 2 * j + (tmp & 0x1);
912        tmp = tmp >> 1;
913    }
914    return j;
915}
916
917////////////////////////////////////////////////////////////////////////////////////////
918// This function perform the in place (direct or inverse) FFT on the N data points
919// contained in the distributed buffers x[nclusters][points_per_cluster].
920// It handles the (N) points 1D array as a (rootN*rootN) points 2D array.
921// 1) it transpose (rootN/nthreads ) rows from x to tmp.
922// 2) it make (rootN/nthreads) FFT on the tmp rows and apply the twiddle factor.
923// 3) it transpose (rootN/nthreads) columns from tmp to x.
924// 4) it make (rootN/nthreads) FFT on the x rows.
925// It calls the FFTRow() 2*(rootN/nthreads) times to perform the in place FFT
926// on the rootN points contained in a row.
927////////////////////////////////////////////////////////////////////////////////////////
928void FFT1D( int              direction,       // direct 1 / inverse -1
929            double       **  x,               // input & output distributed data points array
930            double       **  tmp,             // auxiliary distributed data points array
931            double        *  upriv,           // local array containing coefs for rootN FFT
932            double       **  twid,            // distributed arrays containing N twiddle factors
933            unsigned int     tid,             // thread continuous index
934            unsigned int     MyFirst, 
935            unsigned int     MyLast )
936{
937    unsigned int j;
938    unsigned long long barrier_start;
939    unsigned long long barrier_stop;
940
941#if DEBUG_FFT1D
942unsigned long long cycle;
943get_cycle( &cycle );
944printf("\n[fft] %s : thread %d enter / first %d / last %d / cycle %d\n",
945__FUNCTION__, tid, MyFirst, MyLast, (unsigned int)cycle );
946#endif
947
948    // transpose (rootN/nthreads) rows from x to tmp
949    Transpose( x , tmp , MyFirst , MyLast );
950
951#if( DEBUG_FFT1D & 1 )
952get_cycle( &cycle );
953printf("\n[fft] %s : thread %d after first transpose / cycle %d\n",
954__FUNCTION__, tid, (unsigned int)cycle );
955if( PRINT_ARRAY ) PrintArray( tmp , N );
956#endif
957
958    // BARRIER
959    get_cycle( &barrier_start );
960    pthread_barrier_wait( &barrier );
961    get_cycle( &barrier_stop );
962    sync_time[tid] = (unsigned int)(barrier_stop - barrier_start);
963
964#if( DEBUG_FFT1D & 1 )
965get_cycle( &cycle );
966printf("\n[fft] %s : thread %d exit barrier after first transpose / cycle %d\n",
967__FUNCTION__, tid, (unsigned int)cycle );
968#endif
969
970    // do FFTs on rows of tmp (i.e. columns of x) and apply twiddle factor
971    for (j = MyFirst; j < MyLast; j++) 
972    {
973        FFTRow( direction , upriv , tmp , j * rootN );
974
975        TwiddleOneCol( direction , j , twid , tmp , j * rootN );
976    } 
977
978#if( DEBUG_FFT1D & 1 )
979printf("\n[fft] %s : thread %d after first twiddle\n", __FUNCTION__, tid);
980if( PRINT_ARRAY ) PrintArray( tmp , N );
981#endif
982
983    // BARRIER
984    get_cycle( &barrier_start );
985    pthread_barrier_wait( &barrier );
986    get_cycle( &barrier_stop );
987
988#if( DEBUG_FFT1D & 1 )
989printf("\n[fft] %s : thread %d exit barrier after first twiddle\n", __FUNCTION__, tid);
990#endif
991
992    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
993
994    // transpose tmp to x
995    Transpose( tmp , x , MyFirst , MyLast );
996
997#if( DEBUG_FFT1D & 1 )
998printf("\n[fft] %s : thread %d after second transpose\n", __FUNCTION__, tid);
999if( PRINT_ARRAY ) PrintArray( x , N );
1000#endif
1001
1002    // BARRIER
1003    get_cycle( &barrier_start );
1004    pthread_barrier_wait( &barrier );
1005    get_cycle( &barrier_stop );
1006
1007#if( DEBUG_FFT1D & 1 )
1008printf("\n[fft] %s : thread %d exit barrier after second transpose\n", __FUNCTION__, tid);
1009#endif
1010
1011    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1012
1013    // do FFTs on rows of x and apply the scaling factor
1014    for (j = MyFirst; j < MyLast; j++) 
1015    {
1016        FFTRow( direction , upriv , x , j * rootN );
1017        if (direction == -1) Scale( x , j * rootN );
1018    }
1019
1020#if( DEBUG_FFT1D & 1 )
1021printf("\n[fft] %s : thread %d after FFT on rows\n", __FUNCTION__, tid);
1022if( PRINT_ARRAY ) PrintArray( x , N );
1023#endif
1024
1025    // BARRIER
1026    get_cycle( &barrier_start );
1027    pthread_barrier_wait( &barrier );
1028    get_cycle( &barrier_stop );
1029
1030#if( DEBUG_FFT1D & 1 )
1031printf("\n[fft] %s : thread %d exit barrier after FFT on rows\n", __FUNCTION__, tid);
1032#endif
1033    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1034
1035    // transpose x to tmp
1036    Transpose( x , tmp , MyFirst , MyLast );
1037
1038#if( DEBUG_FFT1D & 1 )
1039printf("\n[fft] %s : thread %x after third transpose\n", __FUNCTION__, tid);
1040if( PRINT_ARRAY ) PrintArray( x , N );
1041#endif
1042
1043    // BARRIER
1044    get_cycle( &barrier_start );
1045    pthread_barrier_wait( &barrier );
1046    get_cycle( &barrier_stop );
1047
1048#if( DEBUG_FFT1D & 1 )
1049printf("\n[fft] %s : thread %d exit barrier after third transpose\n", __FUNCTION__, tid);
1050#endif
1051
1052    sync_time[tid] += (unsigned int)(barrier_stop - barrier_start);
1053    sync_time[tid] += (long)(barrier_stop - barrier_start);
1054
1055    // copy tmp to x
1056    Copy( tmp , x , MyFirst , MyLast );
1057
1058#if DEBUG_FFT1D
1059printf("\n[fft] %s : thread %d completed\n", __FUNCTION__, tid);
1060if( PRINT_ARRAY ) PrintArray( x , N );
1061#endif
1062
1063}  // end FFT1D()
1064
1065/////////////////////////////////////////////////////////////////////////////////////
1066// This function multiply all points contained in a row (rootN points) of the
1067// x[] array by the corresponding twiddle factor, contained in the u[] array.
1068/////////////////////////////////////////////////////////////////////////////////////
1069void TwiddleOneCol( int             direction, 
1070                    unsigned int    j,              // y coordinate in 2D view of coef array
1071                    double       ** u,              // coef array base address
1072                    double       ** x,              // data array base address
1073                    unsigned int    offset_x )      // first point in N points data array
1074{
1075    unsigned int i;
1076    double       omega_r; 
1077    double       omega_c; 
1078    double       x_r; 
1079    double       x_c;
1080    unsigned int c_id;
1081    unsigned int c_offset;
1082
1083    for (i = 0; i < rootN ; i++)  // loop on the rootN points
1084    {
1085        // get coef
1086        c_id      = (j * rootN + i) / (points_per_cluster);
1087        c_offset  = (j * rootN + i) % (points_per_cluster);
1088        omega_r = u[c_id][2*c_offset];
1089        omega_c = direction * u[c_id][2*c_offset+1];
1090
1091        // access data
1092        c_id      = (offset_x + i) / (points_per_cluster);
1093        c_offset  = (offset_x + i) % (points_per_cluster);   
1094        x_r = x[c_id][2*c_offset]; 
1095        x_c = x[c_id][2*c_offset+1];
1096
1097        x[c_id][2*c_offset]   = omega_r*x_r - omega_c * x_c;
1098        x[c_id][2*c_offset+1] = omega_r*x_c + omega_c * x_r;
1099    }
1100}  // end TwiddleOneCol()
1101
1102////////////////////////////
1103void Scale( double      ** x,           // data array base address
1104            unsigned int   offset_x )   // first point of the row to be scaled
1105{
1106    unsigned int i;
1107    unsigned int c_id;
1108    unsigned int c_offset;
1109
1110    for (i = 0; i < rootN ; i++) 
1111    {
1112        c_id      = (offset_x + i) / (points_per_cluster);
1113        c_offset  = (offset_x + i) % (points_per_cluster);
1114        x[c_id][2*c_offset]     /= N;
1115        x[c_id][2*c_offset + 1] /= N;
1116    }
1117}
1118
1119///////////////////////////////////
1120void Transpose( double      ** src,      // source buffer (array of pointers)
1121                double      ** dest,     // destination buffer (array of pointers)
1122                unsigned int   MyFirst,  // first row allocated to the thread
1123                unsigned int   MyLast )  // last row allocated to the thread
1124{
1125    unsigned int row;               // row index
1126    unsigned int point;             // data point index in a row
1127
1128    unsigned int index_src;         // absolute index in the source N points array
1129    unsigned int c_id_src;          // cluster for the source buffer
1130    unsigned int c_offset_src;      // offset in the source buffer
1131
1132    unsigned int index_dst;         // absolute index in the dest N points array
1133    unsigned int c_id_dst;          // cluster for the dest buffer
1134    unsigned int c_offset_dst;      // offset in the dest buffer
1135
1136   
1137    // scan all data points allocated to the thread
1138    // (between MyFirst row and MyLast row) from the source buffer
1139    // and write these points to the destination buffer
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_src    = row * rootN + point;
1145            c_id_src     = index_src / (points_per_cluster);
1146            c_offset_src = index_src % (points_per_cluster);
1147
1148            index_dst    = point * rootN + row;
1149            c_id_dst     = index_dst / (points_per_cluster);
1150            c_offset_dst = index_dst % (points_per_cluster);
1151
1152            dest[c_id_dst][2*c_offset_dst]   = src[c_id_src][2*c_offset_src];
1153            dest[c_id_dst][2*c_offset_dst+1] = src[c_id_src][2*c_offset_src+1];
1154        }
1155    }
1156}  // end Transpose()
1157
1158//////////////////////////////
1159void Copy( double      ** src,      // source buffer (array of pointers)
1160           double      ** dest,     // destination buffer (array of pointers)
1161           unsigned int   MyFirst,  // first row allocated to the thread
1162           unsigned int   MyLast )  // last row allocated to the thread
1163{
1164    unsigned int row;                  // row index
1165    unsigned int point;                // data point index in a row
1166
1167    unsigned int index;                // absolute index in the N points array
1168    unsigned int c_id;                 // cluster index
1169    unsigned int c_offset;             // offset in local buffer
1170
1171    // scan all data points allocated to the thread
1172    for ( row = MyFirst ; row < MyLast ; row++ )       // loop on the rows
1173    {
1174        for ( point = 0 ; point < rootN ; point++ )    // loop on points in row
1175        {
1176            index    = row * rootN + point;
1177            c_id     = index / (points_per_cluster);
1178            c_offset = index % (points_per_cluster);
1179
1180            dest[c_id][2*c_offset]   = src[c_id][2*c_offset];
1181            dest[c_id][2*c_offset+1] = src[c_id][2*c_offset+1];
1182        }
1183    }
1184}  // end Copy()
1185
1186///////////////////////////////
1187void Reverse( double      ** x, 
1188              unsigned int   offset_x )
1189{
1190    unsigned int j, k;
1191    unsigned int c_id_j;
1192    unsigned int c_offset_j;
1193    unsigned int c_id_k;
1194    unsigned int c_offset_k;
1195
1196    for (k = 0 ; k < rootN ; k++) 
1197    {
1198        j = BitReverse( k );
1199        if (j > k) 
1200        {
1201            c_id_j      = (offset_x + j) / (points_per_cluster);
1202            c_offset_j  = (offset_x + j) % (points_per_cluster);
1203            c_id_k      = (offset_x + k) / (points_per_cluster);
1204            c_offset_k  = (offset_x + k) % (points_per_cluster);
1205
1206            SWAP(x[c_id_j][2*c_offset_j]  , x[c_id_k][2*c_offset_k]);
1207            SWAP(x[c_id_j][2*c_offset_j+1], x[c_id_k][2*c_offset_k+1]);
1208        }
1209    }
1210}
1211
1212/////////////////////////////////////////////////////////////////////////////
1213// This function makes the in-place FFT on all points contained in a row
1214// (i.e. rootN points) of the x[nclusters][points_per_cluster] array.
1215/////////////////////////////////////////////////////////////////////////////
1216void FFTRow( int            direction,  // 1 direct / -1 inverse
1217                double       * u,          // private coefs array
1218                double      ** x,          // array of pointers on distributed buffers
1219                unsigned int   offset_x )  // absolute offset in the x array
1220{
1221    unsigned int     j;
1222    unsigned int     k;
1223    unsigned int     q;
1224    unsigned int     L;
1225    unsigned int     r;
1226    unsigned int     Lstar;
1227    double * u1; 
1228
1229    unsigned int     offset_x1;     // index first butterfly input
1230    unsigned int     offset_x2;     // index second butterfly output
1231
1232    double           omega_r;       // real part butterfy coef
1233    double           omega_c;       // complex part butterfly coef
1234
1235    double           tau_r;
1236    double           tau_c;
1237
1238    double           d1_r;          // real part first butterfly input
1239    double           d1_c;          // imag part first butterfly input
1240    double           d2_r;          // real part second butterfly input
1241    double           d2_c;          // imag part second butterfly input
1242
1243    unsigned int     c_id_1;        // cluster index for first butterfly input
1244    unsigned int     c_offset_1;    // offset for first butterfly input
1245    unsigned int     c_id_2;        // cluster index for second butterfly input
1246    unsigned int     c_offset_2;    // offset for second butterfly input
1247
1248#if DEBUG_ROW
1249unsigned int p;
1250printf("\n[fft] ROW data in / %d points / offset = %d\n", rootN , offset_x );
1251
1252for ( p = 0 ; p < rootN ; p++ )
1253{
1254    unsigned int index    = offset_x + p;
1255    unsigned int c_id     = index / (points_per_cluster);
1256    unsigned int c_offset = index % (points_per_cluster);
1257    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1258}
1259printf("\n");
1260#endif
1261
1262    // This makes the rootN input points reordering
1263    Reverse( x , offset_x ); 
1264
1265#if DEBUG_ROW
1266printf("\n[fft] ROW data after reverse / %d points / offset = %d\n", rootN , offset_x );
1267
1268for ( p = 0 ; p < rootN ; p++ )
1269{
1270    unsigned int index    = offset_x + p;
1271    unsigned int c_id     = index / (points_per_cluster);
1272    unsigned int c_offset = index % (points_per_cluster);
1273    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1274}
1275printf("\n");
1276#endif
1277
1278    // This implements the multi-stages, in place Butterfly network
1279    for (q = 1; q <= M/2 ; q++)     // loop on stages
1280    {
1281        L = 1 << q;       // number of points per subset for current stage
1282        r = rootN / L;    // number of subsets
1283        Lstar = L / 2;
1284        u1 = &u[2 * (Lstar - 1)];
1285        for (k = 0; k < r; k++)     // loop on the subsets
1286        {
1287            offset_x1  = offset_x + (k * L);            // index first point
1288            offset_x2  = offset_x + (k * L + Lstar);    // index second point
1289
1290#if (DEBUG_ROW & 1)
1291printf("\n ### q = %d / k = %d / x1 = %d / x2 = %d\n", q , k , offset_x1 , offset_x2 );
1292#endif
1293            // makes all in-place butterfly(s) for subset
1294            for (j = 0; j < Lstar; j++) 
1295            {
1296                // get coef
1297                omega_r = u1[2*j];
1298                omega_c = direction * u1[2*j+1];
1299
1300                // get d[x1] address and value
1301                c_id_1      = (offset_x1 + j) / (points_per_cluster);
1302                c_offset_1  = (offset_x1 + j) % (points_per_cluster);
1303                d1_r        = x[c_id_1][2*c_offset_1];
1304                d1_c        = x[c_id_1][2*c_offset_1+1];
1305
1306                // get d[x2] address and value
1307                c_id_2      = (offset_x2 + j) / (points_per_cluster);
1308                c_offset_2  = (offset_x2 + j) % (points_per_cluster);
1309                d2_r        = x[c_id_2][2*c_offset_2];
1310                d2_c        = x[c_id_2][2*c_offset_2+1];
1311
1312#if (DEBUG_ROW & 1)
1313printf("\n ### d1_in = (%f , %f) / d2_in = (%f , %f) / coef = (%f , %f)\n", 
1314                d1_r , d1_c , d2_r , d2_c , omega_r , omega_c);
1315#endif
1316                // tau = omega * d[x2]
1317                tau_r = omega_r * d2_r - omega_c * d2_c;
1318                tau_c = omega_r * d2_c + omega_c * d2_r;
1319
1320                // set new value for d[x1] = d[x1] + omega * d[x2]
1321                x[c_id_1][2*c_offset_1]   = d1_r + tau_r;
1322                x[c_id_1][2*c_offset_1+1] = d1_c + tau_c;
1323
1324                // set new value for d[x2] = d[x1] - omega * d[x2]
1325                x[c_id_2][2*c_offset_2]   = d1_r - tau_r;
1326                x[c_id_2][2*c_offset_2+1] = d1_c - tau_c;
1327
1328#if (DEBUG_ROW & 1)
1329printf("\n ### d1_out = (%f , %f) / d2_out = (%f , %f)\n", 
1330                d1_r + tau_r , d1_c + tau_c , d2_r - tau_r , d2_c - tau_c );
1331#endif
1332            }
1333        }
1334    }
1335
1336#if DEBUG_ROW
1337printf("\n[fft] ROW data out / %d points / offset = %d\n", rootN , offset_x );
1338for ( p = 0 ; p < rootN ; p++ )
1339{
1340    unsigned int index    = offset_x + p;
1341    unsigned int c_id     = index / (points_per_cluster);
1342    unsigned int c_offset = index % (points_per_cluster);
1343    printf("%f , %f | ", x[c_id][2*c_offset] , x[c_id][2*c_offset+1] );
1344}
1345printf("\n");
1346#endif
1347
1348}  // end FFTRow()
1349
1350///////////////////////////////////////
1351void PrintArray( double       ** array,
1352                 unsigned int    size ) 
1353{
1354    unsigned int  i;
1355    unsigned int  c_id;
1356    unsigned int  c_offset;
1357
1358    // float display
1359    for (i = 0; i < size ; i++) 
1360    {
1361        c_id      = i / (points_per_cluster);
1362        c_offset  = i % (points_per_cluster);
1363
1364        printf(" %f  %f |", array[c_id][2*c_offset], array[c_id][2*c_offset+1]);
1365
1366        if ( (i+1) % 4 == 0)  printf("\n");
1367    }
1368    printf("\n");
1369}
1370
1371
1372// Local Variables:
1373// tab-width: 4
1374// c-basic-offset: 4
1375// c-file-offsets:((innamespace . 0)(inline-open . 0))
1376// indent-tabs-mode: nil
1377// End:
1378
1379// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4
1380
Note: See TracBrowser for help on using the repository browser.