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