Line data Source code
1 : /* Declares a family of functions useful for single threaded sorting of
2 : POD types in high performance contexts. Example usage:
3 :
4 : #define SORT_NAME sort_double_descend
5 : #define SORT_KEY_T double
6 : #define SORT_BEFORE(a,b) (a)>(b)
7 : #include "util/tmpl/fd_sort.c"
8 :
9 : will create the following API for use in the local compile unit:
10 :
11 : // Sort key[i] for i in [0,cnt) stable in place in a best / average
12 : // worst case of O(N) / O(N^2) / O(N^2) operations.
13 :
14 : double *
15 : sort_double_descend_insert( double * key,
16 : ulong cnt );
17 :
18 : // Return 1 if cnt is a sensible value for the sorting APIs
19 : // (i.e. cnt*sizeof(double)+alignof(double)-1 will not overflow,
20 : // 3*cnt*ceil(log_3(cnt)) will not overflow, etc).
21 :
22 : int sort_double_descend_cnt_valid( ulong cnt );
23 :
24 : // Return the alignment and footprint required for a scratch region
25 : // adequate for sorting up to cnt elements. The results will be
26 : // standard allocator and standard declaration friendly. (E.g. a
27 : // declaration "double scratch[ cnt ];" will be fine as a scratch
28 : // region.)
29 :
30 : ulong sort_double_descend_stable_scratch_align ( void );
31 : ulong sort_double_descend_stable_scratch_footprint( ulong cnt );
32 :
33 : // Sort key[i] for i in [0,cnt) stable in best / average / worst
34 : // case of O(N lg N) / O(N lg N) / O(N lg N) operations. Scratch
35 : // is a scratch workspace of suitable alignment and footprint.
36 : // Returns where the sorted values ended up. Will be at either key
37 : // or (double *)scratch.
38 :
39 : double *
40 : sort_double_descend_stable_fast( double * key,
41 : ulong cnt,
42 : void * scratch );
43 :
44 : // Same as above but does additional copying if necessary such that
45 : // the final result ends up in key. Returns key.
46 :
47 : double *
48 : sort_double_descend_stable( double * key,
49 : ulong cnt,
50 : void * scratch );
51 :
52 : // Sort key[i] for i in [0,cnt) inplace in best / average / worst
53 : // case of O(N lg N) / O(N lg N) / O(N^2) operations. Returns key.
54 :
55 : double *
56 : sort_double_descend_inplace( double * key,
57 : ulong cnt );
58 :
59 : // Partially sort key[i] such that key[rnk] holds rnk and return
60 : // key[rnk] in best / average / worst of O(N) / O(N) / O(N^2)
61 : // operations. Returns key. Assumes key is non-NULL, cnt is valid
62 : // and positive and rnk is in [0,cnt)
63 :
64 : double *
65 : sort_double_descend_select( double * key,
66 : ulong cnt,
67 : ulong rnk );
68 :
69 : // Given a sorted array sorted indexed [0,cnt),
70 : // sort_double_descend_split returns an index i in [0,cnt] such all
71 : // entries in [0,i) are BEFORE query and all entries [i,cnt) are
72 : // NOT BEFORE query.
73 :
74 : ulong
75 : sort_double_descend_split( double const * sorted,
76 : ulong cnt,
77 : double query );
78 :
79 : // sort_double_descend_{stable_fast,stable,inplace}_para is the
80 : // same as // above is adaptively parallelized over the caller
81 : // (typically tpool thread t0) and tpool threads (t0,t1). Assumes
82 : // tpool is valid, tpool threads (t0,t1) are idle (or soon to be
83 : // idle). These are only available if SORT_PARALLEL was requested.
84 :
85 : double *
86 : sort_double_descend_stable_fast_para( fd_tpool_t * tpool, ulong t0, ulong t1,
87 : double * key,
88 : ulong cnt,
89 : void * scratch );
90 :
91 : double *
92 : sort_double_descend_stable_para( fd_tpool_t * tpool, ulong t0, ulong t1,
93 : double * key,
94 : ulong cnt,
95 : void * scratch );
96 :
97 : double *
98 : sort_double_descend_inplace_para( fd_tpool_t * tpool, ulong t0, ulong t1,
99 : double * key,
100 : ulong cnt );
101 :
102 : // sort_double_descend_fast_para has better asymptotically
103 : // parallelization on average (~(N lg N)/T + lg T) but requires
104 : // more stack space on the caller (T^2). seed gives a random seed
105 : // to use (e.g. fd_tickcount()). If stable is non-zero, a stable
106 : // method will be used (and the resulting order will always be
107 : // deterministic). Otherwise a faster unstable method will be used
108 : // (the order in which equal keys end up in the array will depend
109 : // on seed). Assumes tpool is valid and tpool threads (t0,t1) are
110 : // idle (or soon to be idle). Only available if SORT_PARALLEL was
111 : // requested and the target has FD_HAS_ALLOCA.
112 :
113 : double *
114 : sort_double_descend_fast_para( fd_tpool_t * tpool, ulong t0, ulong t1,
115 : double * key,
116 : ulong cnt,
117 : void * scratch,
118 : ulong seed,
119 : int stable );
120 :
121 : It is fine to include this template multiple times in a compilation
122 : unit. Just provide the specification before each inclusion. Various
123 : additional options to tune the methods are described below. */
124 :
125 : /* SORT_NAME gives the name of the function to declare (and the base
126 : name of auxiliary and/or variant functions). */
127 :
128 : #ifndef SORT_NAME
129 : #error "SORT_NAME must be defined"
130 : #endif
131 :
132 : /* SORT_KEY_T gives the POD datatype to sort. */
133 :
134 : #ifndef SORT_KEY_T
135 : #error "SORT_KEY_T must be defined"
136 : #endif
137 :
138 : /* SORT_IDX_T gives the data type used to index the arrays */
139 :
140 : #ifndef SORT_IDX_T
141 50920412 : #define SORT_IDX_T ulong
142 : #endif
143 :
144 : /* SORT_BEFORE(a,b) evaluates to 1 if a<b is strictly true. SAFETY TIP:
145 : This is not a 3-way comparison function! */
146 :
147 : #ifndef SORT_BEFORE
148 328 : #define SORT_BEFORE(a,b) (a)<(b)
149 : #endif
150 :
151 : /* SORT_MERGE_THRESH / SORT_QUICK_THRESH give largest merge/quick
152 : partition size where insertion sort will be used. Should be at least
153 : 2 (and probably larger than one might expect to get optimal
154 : performance). */
155 :
156 : #ifndef SORT_MERGE_THRESH
157 0 : #define SORT_MERGE_THRESH 32
158 : #endif
159 :
160 : #ifndef SORT_QUICK_THRESH
161 : #define SORT_QUICK_THRESH 32
162 : #endif
163 :
164 : /* SORT_QUICK_ORDER_STYLE selects the method used for ordering two keys.
165 : Roughly, say 1 for floating point keys (see quick method for
166 : details). */
167 :
168 : #ifndef SORT_QUICK_ORDER_STYLE
169 : #define SORT_QUICK_ORDER_STYLE 0
170 : #endif
171 :
172 : /* SORT_QUICK_SWAP_MINIMIZE indicates that quick sort should eat the
173 : non-deterministic branch cost to avoid no-op swaps. This is useful
174 : if the key_t has a large memory footprint such that the no-op is
175 : actually a larger cost than a branch mispredict. */
176 :
177 : #ifndef SORT_QUICK_SWAP_MINIMIZE
178 : #define SORT_QUICK_SWAP_MINIMIZE 0
179 : #endif
180 :
181 : /* SORT_PARALLEL will generate thread parallel versions of many sorts.
182 : Requires tpool. */
183 :
184 : #ifndef SORT_PARALLEL
185 : #define SORT_PARALLEL 0
186 : #endif
187 :
188 : /* SORT_OVERSAMPLE_RATIO indicates the amount of oversampling fast
189 : parallel sorts should use to fine pivots. Only relevant if
190 : SORT_PARALLEL is true and FD_HAS_ALLOCA. */
191 :
192 : #ifndef SORT_OVERSAMPLE_RATIO
193 0 : #define SORT_OVERSAMPLE_RATIO 5
194 : #endif
195 :
196 : /* SORT_IDX_IF(c,t,f) returns sort_idx t if c is non-zero and sort_idx f o.w. */
197 :
198 : #ifndef SORT_IDX_IF
199 50404951 : #define SORT_IDX_IF(c,t,f) ((SORT_IDX_T)fd_ulong_if( (c), (ulong)(t), (ulong)(f) ))
200 : #endif
201 :
202 : /* SORT_FN_ATTR applies extra function attributes. */
203 :
204 : #ifndef SORT_FN_ATTR
205 : #define SORT_FN_ATTR
206 : #endif
207 :
208 : /* 0 - local use only
209 : 1 - library header declaration
210 : 2 - library implementation */
211 :
212 : #ifndef SORT_IMPL_STYLE
213 : #define SORT_IMPL_STYLE 0
214 : #endif
215 :
216 : /* Implementation *****************************************************/
217 :
218 : #if SORT_IMPL_STYLE==0 /* local use only */
219 : #define SORT_STATIC FD_FN_UNUSED static
220 : #else /* library header and/or implementation */
221 : #define SORT_STATIC
222 : #endif
223 :
224 19202 : #define SORT_(x)FD_EXPAND_THEN_CONCAT3(SORT_NAME,_,x)
225 :
226 : #if SORT_IMPL_STYLE!=2 /* need header */
227 :
228 : #if SORT_PARALLEL
229 : #include "../tpool/fd_tpool.h"
230 : #else
231 : #include "../bits/fd_bits.h"
232 : #endif
233 :
234 : FD_PROTOTYPES_BEGIN
235 :
236 : FD_FN_CONST static inline int
237 0 : SORT_(cnt_valid)( SORT_IDX_T cnt ) {
238 0 : /* Written funny for supporting different signed idx types without
239 0 : generating compiler warnings. */
240 0 : SORT_IDX_T max = (SORT_IDX_T)fd_ulong_min( (-alignof(SORT_KEY_T))/sizeof(SORT_KEY_T), /* ==floor( (2^64-align-1)/sz ) */
241 0 : (ulong)(SORT_IDX_T)((1UL<<57)-1UL) ); /* ==SORT_IDX_MAX as ulong */
242 0 : return (!cnt) | ((((SORT_IDX_T)0)<cnt) & (cnt<max)) | (cnt==max);
243 0 : }
244 :
245 0 : FD_FN_CONST static inline ulong SORT_(stable_scratch_align) ( void ) { return alignof(SORT_KEY_T); }
246 0 : FD_FN_CONST static inline ulong SORT_(stable_scratch_footprint)( SORT_IDX_T cnt ) { return sizeof (SORT_KEY_T)*(ulong)cnt; }
247 :
248 : SORT_FN_ATTR SORT_STATIC SORT_KEY_T *
249 : SORT_(private_merge)( SORT_KEY_T * key,
250 : long cnt,
251 : SORT_KEY_T * tmp );
252 :
253 : SORT_FN_ATTR SORT_STATIC SORT_KEY_T *
254 : SORT_(private_quick)( SORT_KEY_T * key,
255 : SORT_IDX_T cnt );
256 :
257 : SORT_FN_ATTR SORT_STATIC SORT_KEY_T *
258 : SORT_(private_select)( SORT_KEY_T * key,
259 : SORT_IDX_T cnt,
260 : SORT_IDX_T rnk );
261 :
262 : SORT_FN_ATTR SORT_STATIC SORT_KEY_T *
263 : SORT_(insert)( SORT_KEY_T * key,
264 : SORT_IDX_T cnt );
265 :
266 : static inline SORT_KEY_T *
267 : SORT_(stable_fast)( SORT_KEY_T * key,
268 : SORT_IDX_T cnt,
269 0 : void * scratch ) {
270 0 : if( FD_UNLIKELY( cnt<(SORT_IDX_T)2 ) ) return key;
271 0 : SORT_KEY_T * tmp = (SORT_KEY_T *)scratch;
272 0 : return SORT_(private_merge)( key, (long)cnt, tmp );
273 0 : }
274 :
275 : static inline SORT_KEY_T *
276 : SORT_(stable)( SORT_KEY_T * key,
277 : SORT_IDX_T cnt,
278 0 : void * scratch ) {
279 0 : if( FD_UNLIKELY( cnt<(SORT_IDX_T)2 ) ) return key;
280 0 : SORT_KEY_T * tmp = (SORT_KEY_T *)scratch;
281 0 : if( SORT_(private_merge)( key, (long)cnt, tmp )==tmp ) /* 50/50 branch prob */
282 0 : memcpy( key, tmp, sizeof(SORT_KEY_T)*(ulong)cnt );
283 0 : return key;
284 0 : }
285 :
286 : static inline SORT_KEY_T *
287 : SORT_(inplace)( SORT_KEY_T * key,
288 5654 : SORT_IDX_T cnt ) {
289 5654 : return SORT_(private_quick)( key, cnt );
290 5654 : }
291 :
292 : static inline SORT_KEY_T *
293 : SORT_(select)( SORT_KEY_T * key,
294 : SORT_IDX_T cnt,
295 4 : SORT_IDX_T rnk ) {
296 4 : return SORT_(private_select)( key, cnt, rnk );
297 4 : }
298 :
299 : static inline SORT_IDX_T
300 : SORT_(split)( SORT_KEY_T const * sorted,
301 : SORT_IDX_T cnt,
302 0 : SORT_KEY_T query ) {
303 0 : SORT_IDX_T j = (SORT_IDX_T)0;
304 0 : SORT_IDX_T k = (SORT_IDX_T)cnt;
305 0 : for(;;) {
306 0 : SORT_IDX_T n = k-j;
307 0 : if( FD_UNLIKELY( n<(SORT_IDX_T)1 ) ) break;
308 :
309 : /* At this point, entries [0,j) are known "lt" query (i.e. before),
310 : entries [j,k) are unknown (this range is non-empty and sorted),
311 : entries [k,n) are known "geq" query (i.e. not before)
312 : Test the entry in the middle of the unknown range. */
313 :
314 0 : SORT_IDX_T m = j + (n>>1);
315 :
316 0 : int c = SORT_BEFORE( sorted[ m ], query );
317 :
318 : /* At this point:
319 : If c is 1, entry m is before query. As such, entries [0,m] are all known "lt" query and range (m,k) is unknown now.
320 : If c is 0, entry m is not before query. As such, entries [m,n) are all known "geq" query and range [j,m) is unknown now. */
321 :
322 0 : j = c ? (m+(SORT_IDX_T)1) : j; /* cmov */
323 0 : k = c ? k : m; /* cmov */
324 0 : }
325 :
326 : /* At this point, [0,j) are known "lt" query and [k,n) are known
327 : "geq" query and j==k such that [j,k) is an empty range. */
328 :
329 0 : return k;
330 0 : }
331 :
332 : #if SORT_PARALLEL
333 :
334 : SORT_STATIC FD_MAP_REDUCE_PROTO( SORT_(private_merge_para) );
335 : SORT_STATIC FD_FOR_ALL_PROTO ( SORT_(private_memcpy_para) );
336 :
337 : static inline SORT_KEY_T *
338 : SORT_(stable_fast_para)( fd_tpool_t * tpool,
339 : ulong t0,
340 : ulong t1,
341 : SORT_KEY_T * key,
342 : SORT_IDX_T cnt,
343 0 : void * scratch ) {
344 :
345 : /* The wallclock time of the below for N keys and T threads is
346 : roughly:
347 :
348 : alpha N + beta N (ln N - ln T) / T + gamma ln T
349 :
350 : where the first term represents the cost of the final sort rounds
351 : (which are parallelized over increasingly fewer threads than the
352 : initial rounds ... wallclock sums up to something proportional to N
353 : in the limit of large T), the second term represents the cost of
354 : the initial parallel rounds (which are parallelized over T threads)
355 : and the last term represents the wallclock to dispatch and
356 : synchronize the threads. Optimizing T for a given N yields:
357 :
358 : -beta N (ln N - ln T_opt) / T_opt^2 - beta N / T_opt^2 + gamma / T_opt = 0
359 :
360 : Solving for T_opt in the limit N >> T_opt yields:
361 :
362 : T_opt = (beta/gamma) N ln N
363 :
364 : where the ratio beta is roughly the merge pass marginal wallclock
365 : per key and gamma is proportional to the marginal wallclock to
366 : start and stop a thread. Since these are typically used over a
367 : domain of moderate N, we can approximate ln N by ln N_ref,
368 : yielding:
369 :
370 : T_opt ~ N / thresh
371 :
372 : where thresh ~ gamma / (beta ln N_ref). We use an empirical value
373 : such threads will have at least one normal page of work of keys to
374 : process in a block typically. */
375 :
376 0 : if( FD_UNLIKELY( cnt<2L ) ) return key;
377 :
378 0 : static ulong const thresh = (4096UL + sizeof(SORT_KEY_T)-1UL) / sizeof(SORT_KEY_T); /* ceil(page_sz/key_sz) */
379 0 : ulong t_cnt = fd_ulong_min( t1 - t0, ((ulong)cnt) / thresh );
380 0 : if( FD_UNLIKELY( t_cnt<2UL ) ) return SORT_(stable_fast)( key, cnt, scratch );
381 0 : t1 = t0 + t_cnt;
382 :
383 0 : SORT_KEY_T * out[1];
384 0 : FD_MAP_REDUCE( SORT_(private_merge_para), tpool,t0,t1, 0L,(long)cnt, out, key, scratch );
385 0 : return out[0];
386 0 : }
387 :
388 : static inline SORT_KEY_T *
389 : SORT_(stable_para)( fd_tpool_t * tpool,
390 : ulong t0,
391 : ulong t1,
392 : SORT_KEY_T * key,
393 : SORT_IDX_T cnt,
394 0 : void * scratch ) {
395 :
396 : /* This works the same as the above but does a parallel memcpy if the
397 : result doesn't end up in the correct place. */
398 :
399 0 : if( FD_UNLIKELY( cnt<2L ) ) return key;
400 :
401 0 : static ulong const thresh = (4096UL + sizeof(SORT_KEY_T)-1UL) / sizeof(SORT_KEY_T); /* ceil(page_sz/key_sz) */
402 0 : ulong t_cnt = fd_ulong_min( t1 - t0, ((ulong)cnt) / thresh );
403 0 : if( FD_UNLIKELY( t_cnt<2UL ) ) return SORT_(stable)( key, cnt, scratch );
404 0 : t1 = t0 + t_cnt;
405 :
406 0 : SORT_KEY_T * out[1];
407 0 : FD_MAP_REDUCE( SORT_(private_merge_para), tpool,t0,t1, 0L,(long)cnt, out, key, scratch );
408 0 : if( out[0]!=key ) FD_FOR_ALL( SORT_(private_memcpy_para), tpool,t0,t1, 0L,(long)cnt, key, scratch ); /* 50/50 branch prob */
409 0 : return key;
410 0 : }
411 :
412 : SORT_FN_ATTR SORT_STATIC void
413 : SORT_(private_quick_node)( void * _tpool,
414 : ulong t0, ulong t1,
415 : void * _args,
416 : void * _reduce, ulong _stride,
417 : ulong _l0, ulong _l1,
418 : ulong _m0, ulong _m1,
419 : ulong _n0, ulong _n1 );
420 :
421 : static inline SORT_KEY_T *
422 : SORT_(inplace_para)( fd_tpool_t * tpool,
423 : ulong t0,
424 : ulong t1,
425 : SORT_KEY_T * key,
426 0 : SORT_IDX_T cnt ) {
427 0 : if( FD_UNLIKELY( cnt<(SORT_IDX_T)2 ) ) return key;
428 0 : SORT_(private_quick_node)( tpool,t0,t1, (void *)key, (void *)(ulong)cnt, 0UL, 0UL, 0UL, 0UL, 0UL, 0UL, 0UL );
429 0 : return key;
430 0 : }
431 :
432 : #if FD_HAS_ALLOCA
433 :
434 : SORT_FN_ATTR SORT_STATIC SORT_KEY_T *
435 : SORT_(fast_para)( fd_tpool_t * tpool, ulong t0, ulong t1,
436 : SORT_KEY_T * key,
437 : SORT_IDX_T cnt,
438 : void * scratch,
439 : ulong seed,
440 : int stable );
441 :
442 : #endif /* FD_HAS_ALLOCA */
443 :
444 : #endif /* SORT_PARALLEL */
445 :
446 : FD_PROTOTYPES_END
447 :
448 : #endif
449 :
450 : #if SORT_IMPL_STYLE!=1 /* need implementations (assumes header already included) */
451 :
452 : SORT_FN_ATTR SORT_KEY_T *
453 : SORT_(insert)( SORT_KEY_T * key,
454 13543 : SORT_IDX_T cnt ) {
455 149953 : for( SORT_IDX_T i=((SORT_IDX_T)1); i<cnt; i++ ) {
456 136410 : SORT_KEY_T key_i = key[i];
457 136410 : SORT_IDX_T j = i;
458 292649 : while( j ) {
459 284919 : SORT_IDX_T k = j - ((SORT_IDX_T)1);
460 284919 : SORT_KEY_T key_j = key[k]; if( !(SORT_BEFORE( key_i, key_j )) ) break;
461 156239 : key[j] = key_j;
462 156239 : j = k;
463 156239 : }
464 136410 : key[j] = key_i;
465 136410 : }
466 13543 : return key;
467 13543 : }
468 :
469 : SORT_FN_ATTR static void
470 : SORT_(private_merge_pass)( SORT_KEY_T const * key_l, long cnt_l,
471 : SORT_KEY_T const * key_r, long cnt_r,
472 0 : SORT_KEY_T * key_m ) {
473 0 : long i = 0L;
474 0 : long j = 0L;
475 0 : long k = 0L;
476 0 : # if 1 /* Minimal C language operations */
477 0 : for(;;) { /* Note that cnt_left>0 and cnt_right>0 as cnt > SORT_MERGE_THRESH >= 1 at this point */
478 0 : if( SORT_BEFORE( key_r[j], key_l[i] ) ) {
479 0 : key_m[k++] = key_r[j++];
480 0 : if( j>=cnt_r ) {
481 0 : memcpy( key_m + k, key_l + i, sizeof(SORT_KEY_T)*(ulong)(cnt_l-i) ); /* append left stragglers (at least one) */
482 0 : break;
483 0 : }
484 0 : } else {
485 0 : key_m[k++] = key_l[i++];
486 0 : if( i>=cnt_l ) {
487 0 : memcpy( key_m + k, key_r + j, sizeof(SORT_KEY_T)*(ulong)(cnt_r-j) ); /* append right stragglers (at least one) */
488 0 : break;
489 0 : }
490 0 : }
491 0 : }
492 : # else /* Branch free variant (Pyth investigations suggested the above is better) */
493 : while( ((i<cnt_l) & (j<cnt_r)) ) {
494 : SORT_KEY_T ki = key_l[ i ];
495 : SORT_KEY_T kj = key_r[ j ];
496 : int use_left = !(SORT_BEFORE( kj, ki ));
497 : key_m[ k ] = use_left ? ki : kj; /* cmov ideally */
498 : i += (long) use_left;
499 : j += (long)!use_left;
500 : k++;
501 : }
502 : if( i<cnt_l ) memcpy( key_m + k, key_l + i, sizeof(SORT_KEY_T)*(ulong)(cnt_l-i) );
503 : else if( j<cnt_r ) memcpy( key_m + k, key_r + j, sizeof(SORT_KEY_T)*(ulong)(cnt_r-j) );
504 : # endif
505 0 : }
506 :
507 : SORT_FN_ATTR SORT_KEY_T *
508 : SORT_(private_merge)( SORT_KEY_T * key,
509 : long cnt,
510 0 : SORT_KEY_T * tmp ) {
511 :
512 : /* FIXME: USE IMPLICIT RECURSION ALA QUICK BELOW? */
513 :
514 : /* If below threshold, use insertion sort */
515 :
516 0 : if( cnt<=((long)(SORT_MERGE_THRESH)) ) return SORT_(insert)( key, (SORT_IDX_T)cnt );
517 :
518 : /* Otherwise, break input in half and sort the halves */
519 :
520 0 : SORT_KEY_T * key_l = key;
521 0 : SORT_KEY_T * tmp_l = tmp;
522 0 : long cnt_l = cnt >> 1;
523 0 : SORT_KEY_T * in_l = SORT_(private_merge)( key_l, cnt_l, tmp_l );
524 :
525 0 : SORT_KEY_T * key_r = key + cnt_l;
526 0 : SORT_KEY_T * tmp_r = tmp + cnt_l;
527 0 : long cnt_r = cnt - cnt_l;
528 0 : SORT_KEY_T * in_r = SORT_(private_merge)( key_r, cnt_r, tmp_r );
529 :
530 : /* Merge in_l / in_r */
531 :
532 0 : SORT_KEY_T * out = (in_l==key) ? tmp : key; /* If in_l overlaps with key, merge into tmp. Otherwise, merge into key */
533 :
534 : /* Note that in_l does not overlap with out at this point. in_r might
535 : overlap with the right half of out but the merge pass is fine for
536 : that case. */
537 :
538 0 : SORT_(private_merge_pass)( in_l, cnt_l, in_r, cnt_r, out );
539 :
540 0 : return out;
541 0 : }
542 :
543 : /* This uses a dual pivot quick sort for better theoretical and
544 : practical mojo. */
545 :
546 : SORT_FN_ATTR SORT_KEY_T *
547 : SORT_(private_quick)( SORT_KEY_T * key,
548 5654 : SORT_IDX_T cnt ) {
549 5654 : SORT_IDX_T stack[ 4UL*8UL*sizeof(SORT_IDX_T) ]; /* See note below on sizing */
550 5654 : ulong stack_cnt = 0UL;
551 :
552 5654 : stack[ stack_cnt++ ] = (SORT_IDX_T)0;
553 5654 : stack[ stack_cnt++ ] = cnt;
554 24008 : for(;;) {
555 :
556 : /* Pop the next partition to process */
557 :
558 24008 : if( !stack_cnt ) return key; /* All done */
559 18354 : SORT_IDX_T h = stack[ --stack_cnt ];
560 18354 : SORT_IDX_T l = stack[ --stack_cnt ];
561 :
562 : /* [l,h) is the partition we are sorting. If this partition is
563 : small enough, sort it via insertion sort. */
564 :
565 18354 : SORT_IDX_T n = h-l;
566 18354 : if( FD_LIKELY( n <= ((SORT_IDX_T)(SORT_QUICK_THRESH)) ) ) {
567 13540 : SORT_(insert)( key+l, n );
568 13540 : continue;
569 13540 : }
570 :
571 : /* This partition is too large to insertion sort. Pick two pivots,
572 : sort the partition into a left, center and right partition and
573 : then recursively sort those partitions. We initially use a
574 : simple choice of pivots that is ideal for nearly sorted data (in
575 : either direction) with a little extra sampling to improve the
576 : partitioning for randomly ordered data. There is a possibility
577 : that this deterministic choice of pivots will not succeed in
578 : making two or more non-empty partitions; we detect and correct
579 : that below. */
580 :
581 : /* Initial pivot selection */
582 :
583 4814 : SORT_KEY_T p1;
584 4814 : SORT_KEY_T p2;
585 4814 : do {
586 4814 : SORT_IDX_T n3 = n / (SORT_IDX_T)3; /* magic multiply under the hood typically */
587 4814 : SORT_IDX_T h1 = h - (SORT_IDX_T)1;
588 4814 : SORT_KEY_T p0 = key[ l ];
589 4814 : /**/ p1 = key[ l + n3 ];
590 4814 : /**/ p2 = key[ h1 - n3 ];
591 4814 : SORT_KEY_T p3 = key[ h1 ];
592 4814 : # if SORT_QUICK_ORDER_STYLE==0
593 : /* Generates better code for integral types (branchless via stack) */
594 4814 : SORT_KEY_T _k[2]; ulong _c;
595 24070 : # define ORDER(k0,k1) _k[0] = k0; _k[1] = k1; _c = (ulong)(SORT_BEFORE(k1,k0)); k0 = _k[_c]; k1 = _k[_c^1UL];
596 : # else
597 : /* Generates much better code for floating point types (branchless
598 : via min / max floating point instructions) */
599 : SORT_KEY_T _k0; SORT_KEY_T _k1; int _c;
600 : # define ORDER(k0,k1) _c = (SORT_BEFORE(k1,k0)); _k0 = _c ? k1 : k0; _k1 = _c ? k0 : k1; k0 = _k0; k1 = _k1;
601 : # endif
602 :
603 4814 : ORDER(p0,p2); ORDER(p1,p3);
604 4814 : ORDER(p0,p1); ORDER(p2,p3);
605 4814 : ORDER(p1,p2);
606 4814 : # undef ORDER
607 4814 : } while(0);
608 :
609 4881 : retry: /* Three way partitioning (silly language restriction) */;
610 4881 : SORT_IDX_T i = l;
611 4881 : SORT_IDX_T j = l;
612 4881 : SORT_IDX_T k = h-(SORT_IDX_T)1;
613 50400541 : do { /* Note [j,k] is non-empty (look ma ... no branches!) */
614 :
615 : /* At this point, p1 <= p2 and:
616 : - keys [l,i) are before p1 (left partition)
617 : - keys [i,j) are in [p1,p2] (center partition)
618 : - keys [j,k] are not yet partitioned (region is non-empty)
619 : - keys (k,h) are after p2 (right partition)
620 : Decide where key j should be moved. */
621 :
622 50400541 : SORT_KEY_T kj = key[j];
623 :
624 50400541 : int to_left = (SORT_BEFORE( kj, p1 ));
625 50400541 : int to_right = (SORT_BEFORE( p2, kj ));
626 50400541 : SORT_IDX_T m = SORT_IDX_IF( to_right, k, SORT_IDX_IF( to_left, i, j ) );
627 :
628 : # if SORT_QUICK_SWAP_MINIMIZE
629 0 : if( FD_LIKELY( j!=h ) ) { key[j] = key[m]; key[m] = kj; } /* ~2/3 prob */
630 : # else
631 : /**/ key[j] = key[m]; key[m] = kj;
632 : # endif
633 :
634 50400541 : i += (SORT_IDX_T) to_left;
635 50400541 : j += (SORT_IDX_T)(to_right^1); /* to_left_or_to_center <=> !to_right */
636 50400541 : k -= (SORT_IDX_T) to_right;
637 50400541 : } while( j<=k );
638 :
639 : /* Schedule recursion */
640 :
641 4881 : if( FD_UNLIKELY( (j-i)==n ) ) {
642 :
643 : /* All the keys ended up in the center partition. If p1<p2, we
644 : had an unlucky choice of pivots such that p1==min(key[i]) and
645 : p2==max(key[i]) for i in [l,h). Since we picked the keys
646 : deterministically above, we would have the same bad luck the
647 : next time we processed this partition. But, because there are
648 : at least two distinct elements in the partition, there is a
649 : choice of pivots that will make a non-trivial partition, so we
650 : need to redo this partition with different pivots.
651 :
652 : Randomly picking a new pivot pair in [l,h) has probability
653 : between ~50% and ~100%. The worst case is half the keys in
654 : the partition equal to p1 and the other half equal to p2; the
655 : best case exactly one key is equal to p1 / p2 and all other
656 : keys are in (p1,p2] / [p1,p2).
657 :
658 : The optimal handling of the worst case though is just to set
659 : p1=p2 (p2=p1). This will pull the keys equal to p1 (p2) into
660 : the left (right) partition and the remaining keys into the
661 : center partition; the right (left) partition will be empty.
662 : Thus, this is guaranteed to yield a non-trivial partition of
663 : the center partition but may not yield as load balanced set of
664 : partitions as random for the next iteration. We go with the
665 : deterministic option for simplicity below. */
666 :
667 471 : if( FD_UNLIKELY( (SORT_BEFORE( p1, p2 )) ) ) { p1 = p2; goto retry; }
668 :
669 : /* p1==p2 here ... all keys in this partition are the same
670 : We don't need to recurse further in this partition. */
671 :
672 4410 : } else {
673 :
674 : /* Between [1,n-1] keys ended up in the center partition such that
675 : at least 1 key ended up in another partition. We push all
676 : three partitions onto the partition stack for recursion. The
677 : order of the pushes is important to bound the maximum stack
678 : usage. Below, we push the partitions in descending order of
679 : their size, such that we process the smallest next iteration.
680 : For this, the same O(log_2 cnt) max recursion depth applies as
681 : quicksort even though this does a three way partitioning.
682 : That is, let fl/fc/fr be the fraction of keys in the left,
683 : center, right partitions. At this point, fl is in [0,1), fc
684 : is in (0,1) and fr is in [0,1) and fl+fc+fr=1. The smallest
685 : partition is what we process next and, as the smallest, its
686 : fraction of the keys is at most 0.5(1-fc)<=0.5. This means
687 : that at each level, the size of the list we are sorting is
688 : no greater than half the size of the previous level. We do
689 : push up to two partitions onto the stack (four SORT_IDX_T)
690 : each level of recursion though instead of one like normal
691 : quick sort though. */
692 :
693 : /* Create an array representing the lower, center and upper
694 : partitions. Then sort this in order of the size of the
695 : partition, so that we can then push them onto the stack in
696 : order of their size. */
697 :
698 4410 : struct { SORT_IDX_T l; SORT_IDX_T h; SORT_IDX_T n; } part[3];
699 :
700 4410 : part[0].l = l; part[0].h = i; part[0].n = i-l;
701 4410 : part[1].l = i; part[1].h = j; part[1].n = SORT_IDX_IF( SORT_BEFORE( p1, p2 ), j-i, (SORT_IDX_T)0 );
702 4410 : part[2].l = j; part[2].h = h; part[2].n = h-j;
703 :
704 13230 : # define ORDER(i,j) fd_swap_if( part[j].n > part[i].n, part[i], part[j] )
705 :
706 4410 : ORDER(0,2); ORDER(0,1); ORDER(1,2);
707 :
708 4410 : # undef ORDER
709 :
710 : /* Branchlessly push the partitions onto the stack, from largest
711 : to smallest. Skip empty partitions. */
712 :
713 4410 : ulong push0 = (ulong)(!!part[0].n); FD_COMPILER_FORGET( push0 );
714 4410 : ulong push1 = (ulong)(!!part[1].n); FD_COMPILER_FORGET( push1 );
715 4410 : ulong push2 = (ulong)(!!part[2].n); FD_COMPILER_FORGET( push2 );
716 :
717 4410 : stack[ stack_cnt ] = part[0].l; stack_cnt += push0;
718 4410 : stack[ stack_cnt ] = part[0].h; stack_cnt += push0;
719 4410 : stack[ stack_cnt ] = part[1].l; stack_cnt += push1;
720 4410 : stack[ stack_cnt ] = part[1].h; stack_cnt += push1;
721 4410 : stack[ stack_cnt ] = part[2].l; stack_cnt += push2;
722 4410 : stack[ stack_cnt ] = part[2].h; stack_cnt += push2;
723 4410 : }
724 4881 : }
725 : /* never get here */
726 5654 : }
727 :
728 : /* This works identical to sort_private_quick. Only differences
729 : relevant to selection are commented (in short, we only recurse on the
730 : partition that could contain rank). See above for comments on the
731 : algo. */
732 :
733 : SORT_FN_ATTR SORT_KEY_T *
734 : SORT_(private_select)( SORT_KEY_T * key,
735 : SORT_IDX_T cnt,
736 4 : SORT_IDX_T rnk ) {
737 4 : SORT_IDX_T l = (SORT_IDX_T)0;
738 4 : SORT_IDX_T h = cnt;
739 4 : for(;;) {
740 :
741 : /* The partition [l,h) contains rnk. If this partition is small
742 : enough, sort it via insertion sort. FIXME: probably could
743 : truncate the insertion sort early for further optimization (e.g.
744 : if rnk==l / rnk==h-1, we only need to find the min/max for O(n)
745 : operation). */
746 :
747 4 : SORT_IDX_T n = h-l;
748 4 : if( FD_LIKELY( n <= ((SORT_IDX_T)(SORT_QUICK_THRESH)) ) ) { SORT_(insert)( key+l, n ); return key; }
749 :
750 0 : SORT_KEY_T p1;
751 0 : SORT_KEY_T p2;
752 0 : do {
753 0 : SORT_IDX_T n3 = n / (SORT_IDX_T)3;
754 0 : SORT_IDX_T h1 = h - (SORT_IDX_T)1;
755 0 : SORT_KEY_T p0 = key[ l ];
756 0 : /**/ p1 = key[ l + n3 ];
757 0 : /**/ p2 = key[ h1 - n3 ];
758 0 : SORT_KEY_T p3 = key[ h1 ];
759 0 : # if SORT_QUICK_ORDER_STYLE==0
760 0 : SORT_KEY_T _k[2]; ulong _c;
761 0 : # define ORDER(k0,k1) _k[0] = k0; _k[1] = k1; _c = (ulong)(SORT_BEFORE(k1,k0)); k0 = _k[_c]; k1 = _k[_c^1UL];
762 : # else
763 : SORT_KEY_T _k0; SORT_KEY_T _k1; int _c;
764 : # define ORDER(k0,k1) _c = (SORT_BEFORE(k1,k0)); _k0 = _c ? k1 : k0; _k1 = _c ? k0 : k1; k0 = _k0; k1 = _k1;
765 : # endif
766 :
767 0 : ORDER(p0,p2); ORDER(p1,p3);
768 0 : ORDER(p0,p1); ORDER(p2,p3);
769 0 : ORDER(p1,p2);
770 0 : # undef ORDER
771 0 : } while(0);
772 :
773 0 : retry: /* Silly language restriction */;
774 0 : SORT_IDX_T i = l;
775 0 : SORT_IDX_T j = l;
776 0 : SORT_IDX_T k = h-(SORT_IDX_T)1;
777 0 : do {
778 :
779 0 : SORT_KEY_T kj = key[j];
780 :
781 0 : int to_left = (SORT_BEFORE( kj, p1 ));
782 0 : int to_right = (SORT_BEFORE( p2, kj ));
783 0 : SORT_IDX_T m = SORT_IDX_IF( to_right, k, SORT_IDX_IF( to_left, i, j ) );
784 :
785 : # if SORT_QUICK_SWAP_MINIMIZE
786 : if( FD_LIKELY( j!=h ) ) { key[j] = key[m]; key[m] = kj; }
787 : # else
788 : /**/ key[j] = key[m]; key[m] = kj;
789 0 : # endif
790 :
791 0 : i += (SORT_IDX_T) to_left;
792 0 : j += (SORT_IDX_T)(to_right^1);
793 0 : k -= (SORT_IDX_T) to_right;
794 0 : } while( j<=k );
795 :
796 0 : if( FD_UNLIKELY( (j-i)==n ) ) {
797 :
798 0 : if( FD_UNLIKELY( (SORT_BEFORE( p1, p2 )) ) ) { p1 = p2; goto retry; }
799 :
800 : /* p1==p2 here ... all keys in this partition are the same. We
801 : don't need to recurse further in this partition as all keys
802 : would do here. */
803 :
804 0 : return key;
805 0 : }
806 :
807 : /* At this point:
808 : - [l,i) is the left partition,
809 : - [i,j) is the center partition
810 : - [j,h) is the right partition
811 : Between [1,n-1] keys ended up in the center partition such that
812 : at least 1 key ended up in another partition. Recurse on the
813 : partition that contains rnk. As this is typically used in median
814 : finding, we assume that the center partition is the most likely
815 : in the below selection (this can also be done branchlessly). */
816 :
817 0 : SORT_IDX_T l_next = i; SORT_IDX_T h_next = j;
818 0 : if( FD_UNLIKELY( rnk< i ) ) l_next = l, h_next = i;
819 0 : if( FD_UNLIKELY( j <=rnk ) ) l_next = j, h_next = h;
820 0 : l = l_next; h = h_next;
821 0 : }
822 : /* never get here */
823 4 : }
824 :
825 : #if SORT_PARALLEL
826 :
827 : /* This works identical to sort_private_quick. Only differences
828 : relevant to parallelization are commented. See above for comments on
829 : the algo. */
830 :
831 : SORT_FN_ATTR void
832 : SORT_(private_quick_node)( void * _tpool,
833 : ulong t0, ulong t1,
834 : void * _args,
835 : void * _reduce, ulong _stride,
836 : ulong _l0, ulong _l1,
837 : ulong _m0, ulong _m1,
838 0 : ulong _n0, ulong _n1 ) {
839 0 : (void)_stride; (void)_l0; (void)_l1; (void)_m0; (void)_m1; (void)_n0; (void)_n1;
840 :
841 0 : fd_tpool_t * tpool = (fd_tpool_t *) _tpool;
842 0 : SORT_KEY_T * key = (SORT_KEY_T *) _args;
843 0 : SORT_IDX_T cnt = (SORT_IDX_T)(ulong)_reduce;
844 :
845 0 : SORT_IDX_T stack[ 4UL*8UL*sizeof(SORT_IDX_T) ];
846 0 : ulong stack_cnt = 0UL;
847 :
848 0 : ulong wait_stack[ 82 ]; /* See note below for sizing considerations */
849 0 : ulong wait_stack_cnt = 0UL;
850 :
851 0 : stack[ stack_cnt++ ] = (SORT_IDX_T)0;
852 0 : stack[ stack_cnt++ ] = cnt;
853 :
854 0 : while( stack_cnt ) {
855 :
856 0 : SORT_IDX_T h = stack[ --stack_cnt ];
857 0 : SORT_IDX_T l = stack[ --stack_cnt ];
858 :
859 0 : SORT_IDX_T n = h-l;
860 0 : if( FD_LIKELY( n <= ((SORT_IDX_T)(SORT_QUICK_THRESH)) ) ) {
861 0 : SORT_(insert)( key+l, n );
862 0 : continue;
863 0 : }
864 :
865 0 : SORT_KEY_T p1;
866 0 : SORT_KEY_T p2;
867 0 : do {
868 0 : SORT_IDX_T n3 = n / (SORT_IDX_T)3;
869 0 : SORT_IDX_T h1 = h - (SORT_IDX_T)1;
870 0 : SORT_KEY_T p0 = key[ l ];
871 0 : /**/ p1 = key[ l + n3 ];
872 0 : /**/ p2 = key[ h1 - n3 ];
873 0 : SORT_KEY_T p3 = key[ h1 ];
874 0 : # if SORT_QUICK_ORDER_STYLE==0
875 0 : SORT_KEY_T _k[2]; ulong _c;
876 0 : # define ORDER(k0,k1) _k[0] = k0; _k[1] = k1; _c = (ulong)(SORT_BEFORE(k1,k0)); k0 = _k[_c]; k1 = _k[_c^1UL];
877 : # else
878 : SORT_KEY_T _k0; SORT_KEY_T _k1; int _c;
879 : # define ORDER(k0,k1) _c = (SORT_BEFORE(k1,k0)); _k0 = _c ? k1 : k0; _k1 = _c ? k0 : k1; k0 = _k0; k1 = _k1;
880 : # endif
881 :
882 0 : ORDER(p0,p2); ORDER(p1,p3);
883 0 : ORDER(p0,p1); ORDER(p2,p3);
884 0 : ORDER(p1,p2);
885 0 : # undef ORDER
886 0 : } while(0);
887 :
888 0 : retry: /* Silly language restriction */;
889 0 : SORT_IDX_T i = l;
890 0 : SORT_IDX_T j = l;
891 0 : SORT_IDX_T k = h-(SORT_IDX_T)1;
892 0 : do {
893 0 : SORT_KEY_T kj = key[j];
894 :
895 0 : int to_left = (SORT_BEFORE( kj, p1 ));
896 0 : int to_right = (SORT_BEFORE( p2, kj ));
897 0 : SORT_IDX_T m = SORT_IDX_IF( to_right, k, SORT_IDX_IF( to_left, i, j ) );
898 :
899 : # if SORT_QUICK_SWAP_MINIMIZE
900 : if( FD_LIKELY( j!=h ) ) { key[j] = key[m]; key[m] = kj; }
901 : # else
902 : /**/ key[j] = key[m]; key[m] = kj;
903 0 : # endif
904 :
905 0 : i += (SORT_IDX_T) to_left;
906 0 : j += (SORT_IDX_T)(to_right^1);
907 0 : k -= (SORT_IDX_T) to_right;
908 0 : } while( j<=k );
909 :
910 0 : if( FD_UNLIKELY( (j-i)==n ) ) {
911 :
912 0 : if( FD_UNLIKELY( (SORT_BEFORE( p1, p2 )) ) ) { p1 = p2; goto retry; }
913 :
914 0 : } else {
915 :
916 : /* At this point, we have at most 3 partitions to sort and can use
917 : the caller and tpool threads (t0,t1) to sort them. To load
918 : balance this, we sort the partitions by an estimate of how much
919 : work each requires. We estimate this as proportional to:
920 :
921 : work = cnt ceil log_3 cnt if cnt>0 and 0 otherwise
922 :
923 : where cnt is the number of keys in the partition. Since this
924 : is a rough estimate used for load balancing, we don't need a
925 : precise log_3 calculation. With log_3 x ~ 0.633 lg x, we
926 : approximate:
927 :
928 : ceil log_3 x ~ ceil( 5 (lg x) / 8 )
929 : ~ ceil( 5 (floor lg x) / 8)
930 : = floor( (5 (floor lg x) + 7)/8 )
931 : = (5 (msb x) + 7) >> 3
932 :
933 : This approximation is trivial to compute, is monotonic and is
934 : either exact or 1 less for x in (0,2^64). The approximation is
935 : such that cnt==1 / cnt>1 yields zero / positive. As such,
936 : work==0 strictly indicates no additional work is necessary. */
937 :
938 0 : # define WORK_APPROX(cnt) (((ulong)(cnt))*(ulong)((5*fd_ulong_find_msb_w_default( (ulong)(cnt), 0 )+7) >> 3))
939 :
940 0 : struct { SORT_IDX_T l; SORT_IDX_T h; ulong work; } part[3];
941 :
942 0 : part[0].l = l; part[0].h = i; part[0].work = WORK_APPROX( i-l );
943 0 : part[1].l = i; part[1].h = j; part[1].work = fd_ulong_if( SORT_BEFORE( p1, p2 ), WORK_APPROX( j-i ), 0UL );
944 0 : part[2].l = j; part[2].h = h; part[2].work = WORK_APPROX( h-j );
945 :
946 0 : # undef WORK_APPROX
947 :
948 0 : ulong work_remain = part[0].work + part[1].work + part[2].work;
949 :
950 0 : # define ORDER(i,j) fd_swap_if( part[j].work > part[i].work, part[i], part[j] )
951 :
952 0 : ORDER(0,2); ORDER(0,1); ORDER(1,2);
953 :
954 0 : # undef ORDER
955 :
956 0 : for( ulong idx=0UL; idx<3UL; idx++ ) {
957 :
958 : /* At this point, we need to schedule partitions [idx,2] requiring
959 : work_remain total estimated work and we have the caller and
960 : tpool threads (t0,t1) to do this work. Compute which threads
961 : should execute partition idx. */
962 :
963 0 : ulong p_work = part[idx].work;
964 0 : if( FD_UNLIKELY( !p_work ) ) break; /* no work need for this partition and remainin partitions */
965 :
966 0 : SORT_IDX_T p_l = part[idx].l;
967 0 : SORT_IDX_T p_h = part[idx].h;
968 :
969 0 : ulong t_cnt = t1 - t0;
970 0 : if( FD_UNLIKELY( (t_cnt>1UL) & (p_work<work_remain) ) ) {
971 :
972 : /* At this point, we have at least two partitions that need
973 : work to sort and at least two threads remaining for
974 : sorting. Schedule this partition for execution on between
975 : 1 and t_cnt-1 of the remaining threads roughly proportional
976 : to the estimated work to sort this partition versus the
977 : total estimated work remaining.
978 :
979 : For the wait stack, in the worst case, we schedule two
980 : partitions for remote execution covering two thirds of the
981 : keys each round. This implies the wait stack depth is
982 : bounded by:
983 :
984 : 2 ceil log_3 thread_cnt
985 :
986 : Practically, thread_cnt<=FD_TILE_MAX. For
987 : FD_TILE_MAX==1024, 14 is sufficient and 82 is sufficient
988 : for any thread_cnt representable by a ulong. */
989 :
990 0 : ulong ts = t1 - fd_ulong_min( fd_ulong_max( (t_cnt*p_work + (work_remain>>1)) / work_remain, 1UL ), t_cnt-1UL );
991 :
992 0 : fd_tpool_exec( tpool, ts, SORT_(private_quick_node), tpool, ts, t1, key + p_l, (void *)(ulong)(p_h - p_l),
993 0 : 0UL, 0UL, 0UL, 0UL, 0UL, 0UL, 0UL );
994 :
995 0 : wait_stack[ wait_stack_cnt++ ] = ts;
996 :
997 : /* Advance t1 and work_remain for the next iteration */
998 :
999 0 : t1 = ts;
1000 0 : work_remain -= p_work;
1001 :
1002 0 : } else {
1003 :
1004 : /* At this point, we have only one thread available and/or
1005 : this is the only partition remaining that needs sorting.
1006 : Schedule this partition for local execution. */
1007 :
1008 0 : stack[ stack_cnt++ ] = p_l;
1009 0 : stack[ stack_cnt++ ] = p_h;
1010 :
1011 0 : }
1012 0 : }
1013 0 : }
1014 0 : }
1015 :
1016 : /* Wait for any children to finish */
1017 :
1018 0 : while( wait_stack_cnt ) fd_tpool_wait( tpool, wait_stack[ --wait_stack_cnt ] );
1019 0 : }
1020 :
1021 0 : FD_MAP_REDUCE_BEGIN( SORT_(private_merge_para), 1L, 0UL, sizeof(SORT_KEY_T *), 1L ) {
1022 :
1023 0 : SORT_KEY_T ** out = (SORT_KEY_T **)arg[0];
1024 0 : SORT_KEY_T * key = (SORT_KEY_T * )arg[1];
1025 0 : SORT_KEY_T * tmp = (SORT_KEY_T * )arg[2];
1026 :
1027 0 : *out = SORT_(stable_fast)( key + block_i0, (SORT_IDX_T)block_cnt, tmp + block_i0 );
1028 :
1029 0 : } FD_MAP_END {
1030 :
1031 0 : SORT_KEY_T * key = (SORT_KEY_T * )arg[1];
1032 0 : SORT_KEY_T * tmp = (SORT_KEY_T * )arg[2];
1033 0 : SORT_KEY_T ** _in_l = (SORT_KEY_T **)arg[0]; SORT_KEY_T * in_l = *_in_l; long cnt_l = block_is - block_i0;
1034 0 : SORT_KEY_T ** _in_r = (SORT_KEY_T **)_r1; SORT_KEY_T * in_r = *_in_r; long cnt_r = block_i1 - block_is;
1035 :
1036 : /* Merge in_l / in_r (see private_merge above for details) */
1037 :
1038 0 : SORT_KEY_T * out = ((in_l==(key+block_i0)) ? tmp : key) + block_i0; /* cmov */
1039 :
1040 0 : SORT_(private_merge_pass)( in_l, cnt_l, in_r, cnt_r, out );
1041 :
1042 0 : *_in_l = out;
1043 :
1044 0 : } FD_REDUCE_END
1045 :
1046 0 : FD_FOR_ALL_BEGIN( SORT_(private_memcpy_para), 1L ) {
1047 :
1048 0 : SORT_KEY_T * dst = (SORT_KEY_T *)arg[0];
1049 0 : SORT_KEY_T const * src = (SORT_KEY_T const *)arg[1];
1050 :
1051 0 : if( FD_LIKELY( block_cnt ) ) memcpy( dst + block_i0, src + block_i0, sizeof(SORT_KEY_T)*(ulong)block_cnt );
1052 :
1053 0 : } FD_FOR_ALL_END
1054 :
1055 : #if FD_HAS_ALLOCA
1056 :
1057 0 : static FD_FOR_ALL_BEGIN( SORT_(private_cntcpy_para), 1L ) {
1058 0 : ulong tpool_base = (ulong )arg[0];
1059 0 : ulong t_cnt = (ulong )arg[1];
1060 0 : ulong * _key_cnt = (ulong *)arg[2];
1061 0 : SORT_KEY_T const * key = (SORT_KEY_T const *)arg[3];
1062 0 : SORT_KEY_T * tmp = (SORT_KEY_T *)arg[4];
1063 0 : SORT_KEY_T const * pivot = (SORT_KEY_T const *)arg[5];
1064 :
1065 : /* Note keys in [ pivot[t-1], pivot[t] ) for t in [0,t_cnt) are keys
1066 : assigned to thread t for sorting. pivot[-1] / pivot[t_cnt-1] are
1067 : an implied -infinity / +infinity and never explicitly accessed. As
1068 : such, pivot is indexed [0,t_cnt-1). */
1069 :
1070 : /* Allocate and clear a local scratch for counting. We don't count
1071 : directly into key_cnt to avoid false sharing while counting. */
1072 :
1073 0 : ulong * key_cnt = fd_alloca( alignof(ulong), sizeof(ulong)*t_cnt );
1074 :
1075 0 : memset( key_cnt, 0, sizeof(ulong)*t_cnt );
1076 :
1077 0 : for( long i=block_i0; i<block_i1; i++ ) {
1078 0 : SORT_KEY_T ki = key[i];
1079 :
1080 : /* Determine which thread is responsible for subsorting this key */
1081 : /* FIXME: ideally, use a unrolled fixed count iteration for here */
1082 :
1083 0 : ulong l = 0UL;
1084 0 : ulong h = t_cnt;
1085 0 : for(;;) {
1086 0 : ulong n = h - l;
1087 0 : if( n<2UL ) break;
1088 :
1089 : /* At this point, the thread for key i is in [l,h) and this range
1090 : contains at least two threads. Split this range in half. */
1091 :
1092 0 : ulong m = l + (n>>1); /* In [1,t_cnt) */
1093 0 : int c = SORT_BEFORE( ki, pivot[m-1UL] );
1094 :
1095 : /* If ki is before pivot[m], the target thread is in [l,m).
1096 : Otherwise, the target thread is in [m,h). */
1097 :
1098 0 : l = fd_ulong_if( c, l, m );
1099 0 : h = fd_ulong_if( c, m, h );
1100 0 : }
1101 :
1102 : /* At this point, key[i] should be assigned to thread l. Count it
1103 : and copy it into tmp to prepare to scatter. */
1104 :
1105 0 : key_cnt[l]++;
1106 0 : tmp[i] = ki;
1107 0 : }
1108 :
1109 : /* Send the counts back to main thread for partitioning */
1110 :
1111 0 : memcpy( _key_cnt + (tpool_t0-tpool_base)*t_cnt, key_cnt, sizeof(ulong)*t_cnt );
1112 :
1113 0 : } FD_FOR_ALL_END
1114 :
1115 0 : static FD_FOR_ALL_BEGIN( SORT_(private_scatter_para), 1L ) {
1116 0 : ulong tpool_base = (ulong )arg[0];
1117 0 : ulong t_cnt = (ulong )arg[1];
1118 0 : ulong const * _part = (ulong const *)arg[2];
1119 0 : SORT_KEY_T * key = (SORT_KEY_T *)arg[3];
1120 0 : SORT_KEY_T const * tmp = (SORT_KEY_T const *)arg[4];
1121 0 : SORT_KEY_T const * pivot = (SORT_KEY_T const *)arg[5];
1122 :
1123 : /* Receive the key array partitioning from the main thread. Like the
1124 : above, we don't operate on the part array directly to avoid false
1125 : sharing while scattering. */
1126 :
1127 0 : ulong * part = fd_alloca( alignof(ulong), sizeof(ulong)*t_cnt );
1128 :
1129 0 : memcpy( part, _part + (tpool_t0-tpool_base)*t_cnt, sizeof(ulong)*t_cnt );
1130 :
1131 0 : for( long i=block_i0; i<block_i1; i++ ) {
1132 0 : SORT_KEY_T ki = tmp[i];
1133 :
1134 : /* Determine which thread is responsible for subsort this key
1135 : (again). Identical considerations to the above. Note: we can
1136 : eliminate this computation if willing to burn O(N) additional
1137 : scratch memory by saving results computed above for use here. */
1138 :
1139 0 : ulong l = 0UL;
1140 0 : ulong h = t_cnt;
1141 0 : for(;;) {
1142 0 : ulong n = h - l;
1143 0 : if( n<2UL ) break;
1144 0 : ulong m = l + (n>>1);
1145 0 : int c = SORT_BEFORE( ki, pivot[m-1UL] );
1146 0 : l = fd_ulong_if( c, l, m );
1147 0 : h = fd_ulong_if( c, m, h );
1148 0 : }
1149 :
1150 : /* Send this key to target thread */
1151 :
1152 0 : key[ part[l]++ ] = ki;
1153 0 : }
1154 :
1155 0 : } FD_FOR_ALL_END
1156 :
1157 0 : static FD_FOR_ALL_BEGIN( SORT_(private_subsort_para), 1L ) {
1158 0 : ulong const * part = (ulong const *)arg[0];
1159 0 : SORT_KEY_T * key = (SORT_KEY_T *)arg[1];
1160 0 : SORT_KEY_T * tmp = (SORT_KEY_T *)arg[2];
1161 0 : int stable = (int )arg[3];
1162 :
1163 0 : ulong j0 = part[ block_i0 ];
1164 0 : ulong j1 = part[ block_i1 ];
1165 :
1166 0 : if( stable ) SORT_(stable) ( key + j0, j1 - j0, tmp + j0 );
1167 0 : else SORT_(inplace)( key + j0, j1 - j0 );
1168 :
1169 0 : } FD_FOR_ALL_END
1170 :
1171 : SORT_FN_ATTR SORT_KEY_T *
1172 : SORT_(fast_para)( fd_tpool_t * tpool, ulong t0, ulong t1,
1173 : SORT_KEY_T * key,
1174 : SORT_IDX_T cnt,
1175 : void * scratch,
1176 : ulong seed,
1177 0 : int stable ) {
1178 :
1179 0 : if( FD_UNLIKELY( cnt<(SORT_IDX_T)2 ) ) return key; /* nothing to do */
1180 :
1181 : /* For the below (if sampling ops are fully threaded), sampling costs
1182 : O(S), the sample sort costs O(S lg(T S)), the downsampling costs
1183 : O(1), the partitioning costs O(T), the counting and scattering cost
1184 : O(N (lg T) / T), the subsorts cost (N/T) lg(N/T) and thread
1185 : synchronization overhead is O(lg T).
1186 :
1187 : Combining all these, assuming the sampling ratio S is a fixed O(1)
1188 : quantity and assuming N (lg T) / T term in the counting/scattering
1189 : approximately cancals the -N lg T / T in the subsorts yields a
1190 : wallclock that scales as:
1191 :
1192 : alpha N (ln N) / T + beta ln T + gamma T
1193 :
1194 : The first term represents the parallelized sorting cost, the second
1195 : term represents the thread dispatch / sync cost, and the last term
1196 : represents the partitioning costs. Minimizing cost with respect
1197 : to T yields:
1198 :
1199 : -alpha (N ln N) / T_opt^2 + beta / T_opt + gamma = 0
1200 :
1201 : whose (positive) solution is:
1202 :
1203 : T_opt ~ (0.5 beta / gamma) [ sqrt( 1 + (4 alpha gamma N ln N)/beta^2) ) - 1 ]
1204 :
1205 : In the limit 4 alpha gamma N ln N / beta^2 >> 1 (that is,
1206 : partitioning costs are much lower than thread start/stop costs), we
1207 : have:
1208 :
1209 : T_opt -> (alpha N ln N) / beta
1210 :
1211 : In the other limit, we have:
1212 :
1213 : T_opt -> sqrt( alpha N ln N / gamma )
1214 :
1215 : The first limit tends to apply in practice. Thus we use:
1216 :
1217 : T_opt ~ (N lg N) / thresh
1218 : ~ floor( (N msb N + N/2 + thresh/2) / thresh)
1219 :
1220 : where thresh is an empirical minimum amount of sorting work to
1221 : justify starting / stopping a thread. (As written below, the gamma
1222 : term is more like gamma T^2, which makes the other limit more like
1223 : the cube root of N ln N but doesn't change the overall
1224 : conclusion.) */
1225 :
1226 0 : ulong thresh = (4096UL + sizeof(SORT_KEY_T)-1UL) / sizeof(SORT_KEY_T);
1227 0 : ulong t_cnt = fd_ulong_min( t1 - t0, (cnt*(ulong)fd_ulong_find_msb( cnt ) + ((cnt+thresh)>>1)) / thresh );
1228 0 : if( FD_UNLIKELY( t_cnt<2UL ) ) return stable ? SORT_(stable)( key, cnt, scratch ) : SORT_(inplace)( key, cnt );
1229 0 : t1 = t0 + t_cnt;
1230 :
1231 : /* At this point, we have at least 2 threads available and at least 2
1232 : items to sort. Sample the keys to get some idea of their
1233 : distribution, sort the samples and downsample the sorted samples
1234 : into pivots that approximately uniformly partition the samples into
1235 : t_cnt groups (and thus approximately partition the keys uniformly
1236 : too). Notes:
1237 :
1238 : - The sampling below is practically uniform IID but it could be
1239 : made more robust with a rejection method (ala fd_rng_ulong_roll).
1240 :
1241 : - Increasing SORT_OVERSAMPLE_RATIO improves the uniformity of the
1242 : partitioning of key space over the inputs but requires more stack
1243 : scratch memory and more overhead.
1244 :
1245 : - All three steps here could be parallelized but it is probably not
1246 : worth it given the overhead for threaded versus the number of
1247 : samples.
1248 :
1249 : - For a stable sort, the result is completely deterministic even
1250 : if the seed provided is non-deterministic. */
1251 :
1252 0 : ulong sample_cnt = t_cnt*(ulong)SORT_OVERSAMPLE_RATIO;
1253 0 : SORT_KEY_T * pivot = fd_alloca( alignof(SORT_KEY_T), sizeof(SORT_KEY_T)*sample_cnt );
1254 :
1255 0 : for( ulong i=0UL; i<sample_cnt; i++ ) pivot[i] = key[ fd_ulong_hash( seed ^ i ) % cnt ];
1256 :
1257 0 : SORT_(inplace)( pivot, sample_cnt );
1258 :
1259 0 : for( ulong i=1UL; i<t_cnt; i++ ) pivot[i-1UL] = pivot[ i*(ulong)SORT_OVERSAMPLE_RATIO ];
1260 :
1261 : /* At this point, keys in [ pivot[t-1], pivot[t] ) should be assigned
1262 : to thread t for thread t's subsort. pivot[-1] / pivot[t_cnt-1] are
1263 : an implicit -infinity / +infinity. Split the input array equally
1264 : over threads and have each thread copy their block of keys into tmp
1265 : and count the number of keys in its block that could be assigned to
1266 : each thread. */
1267 :
1268 0 : # define part(i,j) part[ (i) + t_cnt*(j) ]
1269 :
1270 0 : ulong * part = fd_alloca( alignof(ulong), sizeof(ulong)*t_cnt*t_cnt );
1271 :
1272 0 : SORT_KEY_T * tmp = (SORT_KEY_T *)scratch;
1273 :
1274 0 : FD_FOR_ALL( SORT_(private_cntcpy_para), tpool,t0,t1, 0L,(long)cnt, t0, t_cnt, part, key, tmp, pivot );
1275 :
1276 : /* At this point, part(i,j) is the count of the number of keys in thread
1277 : j's block that were assigned to thread i's subsort. Convert this
1278 : into a partitioning. In-principle this can be parallelized but
1279 : this is rarely worth it practically. */
1280 :
1281 0 : ulong k = 0UL;
1282 0 : for( ulong i=0UL; i<t_cnt; i++ ) {
1283 0 : for( ulong j=0UL; j<t_cnt; j++ ) {
1284 0 : ulong c_ij = part(i,j);
1285 0 : part(i,j) = k;
1286 0 : k += c_ij;
1287 0 : }
1288 0 : }
1289 :
1290 : /* At this point, the range [ part(i,j), part(i+1,j) ) is where keys
1291 : in thread j's block assigned to thread i's subsort should be
1292 : scattered. part(t_cnt,t_cnt-1) is an implied cnt. Scatter the
1293 : keys from tmp back into key in subsorted order. */
1294 :
1295 0 : FD_FOR_ALL( SORT_(private_scatter_para), tpool,t0,t1, 0L,(long)cnt, t0, t_cnt, part, key, tmp, pivot );
1296 :
1297 : /* At this point, keys [ part(i,0), part(i+1,0) ) are the keys assigned
1298 : to thread i for subsorting. part(t_cnt,0) is an implied cnt.
1299 : Since t_cnt>1 though this location exists. We make that explicit
1300 : such that part[i],part[i+1] give the sets of keys each thread
1301 : should sort. Do the thread parallel subsorts to get the keys into
1302 : final sorted order. */
1303 :
1304 0 : part[ t_cnt ] = cnt;
1305 :
1306 0 : FD_FOR_ALL( SORT_(private_subsort_para), tpool,t0,t1, 0L,(long)t_cnt, part, key, tmp, stable );
1307 :
1308 0 : # undef part
1309 :
1310 0 : return key;
1311 0 : }
1312 :
1313 : #endif /* FD_HAS_ALLOCA */
1314 : #endif /* SORT_PARALLEL */
1315 : #endif /* SORT_IMPL_STYLE!=1 */
1316 :
1317 : #undef SORT_
1318 : #undef SORT_STATIC
1319 :
1320 : #undef SORT_IMPL_STYLE
1321 : #undef SORT_FN_ATTR
1322 : #undef SORT_IDX_IF
1323 : #undef SORT_OVERSAMPLE_RATIO
1324 : #undef SORT_PARALLEL
1325 : #undef SORT_QUICK_SWAP_MINIMIZE
1326 : #undef SORT_QUICK_ORDER_STYLE
1327 : #undef SORT_QUICK_THRESH
1328 : #undef SORT_MERGE_THRESH
1329 : #undef SORT_BEFORE
1330 : #undef SORT_IDX_T
1331 : #undef SORT_KEY_T
1332 : #undef SORT_NAME
|
