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

Last change on this file since 639 was 638, checked in by alain, 5 years ago

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