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

Last change on this file since 635 was 635, checked in by alain, 5 months ago

This version is a major evolution: The physical memory allocators,
defined in the kmem.c, ppm.c, and kcm.c files have been modified
to support remote accesses. The RPCs that were previously user
to allocate physical memory in a remote cluster have been removed.
This has been done to cure a dead-lock in case of concurrent page-faults.

This version 2.2 has been tested on a (4 clusters / 2 cores per cluster)
TSAR architecture, for both the "sort" and the "fft" applications.

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