Skip to content

Commit 7defc2a

Browse files
h-2h-2
committed
[INTERNAL] optimizing the radix sort
1 parent 5912105 commit 7defc2a

2 files changed

Lines changed: 110 additions & 69 deletions

File tree

‎src/lambda_indexer.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -402,8 +402,8 @@ createIndexActual(TIndex & index,
402402
indexCreate(index, text, TFibre(),
403403
[&_lastPercent] (uint64_t curPerc)
404404
{
405-
// SEQAN_OMP_PRAGMA(critical(progressBar))
406-
if (TID == 0)
405+
SEQAN_OMP_PRAGMA(critical(progressBar))
406+
// if (TID == 0)
407407
printProgressBar(_lastPercent, curPerc);
408408
});
409409
printProgressBar(_lastPercent, 100);

‎src/radix_inplace.h‎

Lines changed: 108 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,14 @@
4141
#ifndef CORE_INCLUDE_SEQAN_INDEX_RADIX_INPLACE_H_
4242
#define CORE_INCLUDE_SEQAN_INDEX_RADIX_INPLACE_H_
4343

44+
#if defined(_OPENMP) && defined(__GNUC__) && !defined(__clang__)
45+
#include <parallel/algorithm>
46+
#define SORT __gnu_parallel::sort
47+
#else
48+
#define SORT std::sort
49+
#endif
50+
//TODO(h-2): for clang use std::experimenta::parallel if available
51+
4452
namespace SEQAN_NAMESPACE_MAIN
4553
{
4654

@@ -161,52 +169,55 @@ struct _ZeroBucketComparator<TSAValue, Nothing const>
161169
// struct RadixSortContext_
162170
// ----------------------------------------------------------------------------
163171

164-
template <typename TAccessFunctor, // text accessor
165-
typename TOrderFunctor, // For seperate sort of the 0 bucket.
172+
template <typename TSAValue,
173+
typename TText,
166174
typename TSize, // type of depth and bucketCount a.s.o
167175
unsigned Q> // alph size = ValueSize + 1
168176
struct RadixSortContext_
169177
{
170-
static_assert(Q < 256, "Alphabet size must be smaller 256!"); //TODO really?
171-
typedef typename TAccessFunctor::argument_type TSAValue;
172-
typedef typename TAccessFunctor::result_type TOrdValue; // unsigned
178+
typedef typename StringSetLimits<TText const>::Type TLimitsString; // "Nothing" for Strings
179+
typedef RadixTextAccessor<TSAValue, TText> TAccessFunctor;
180+
typedef _ZeroBucketComparator<TSAValue, TLimitsString> TOrderFunctor;
181+
typedef typename TAccessFunctor::result_type TOrdValue;
173182

183+
static_assert(Q < 256, "Alphabet size must be smaller 256!"); //TODO really?
174184
static const unsigned ORACLESIZE = 256;
175-
TAccessFunctor textAccess;
176-
TOrderFunctor comp;
185+
186+
TText const & text;
187+
TAccessFunctor textAccess;
188+
TOrderFunctor comp;
177189

178190
TSize bucketSize[Q];
179191
std::array<TSAValue*,Q> bucketEnd;
180192

181-
RadixSortContext_(TAccessFunctor const & f, TOrderFunctor const & c) :
182-
textAccess(f), comp(c)
193+
RadixSortContext_(TText const & t) :
194+
text(t), textAccess(t), comp(stringSetLimits(t))
183195
{}
184196
};
185197

186-
template <typename TAccessFunctor,
187-
typename TOrderFunctor,
198+
template <typename TSAValue,
199+
typename TText,
188200
typename TSize,
189201
unsigned Q>
190202
inline void
191-
clear(RadixSortContext_<TAccessFunctor, TOrderFunctor, TSize, Q> & context)
203+
clear(RadixSortContext_<TSAValue, TText, TSize, Q> & context)
192204
{
193205
memset(context.bucketSize, 0, sizeof(TSize)*Q);
194206
}
195207

196208
// ----------------------------------------------------------------------------
197-
// Function radixSort()
209+
// Function _radixSort()
198210
// ----------------------------------------------------------------------------
199211

200212
template <typename TSAValue, typename TSize,
201-
typename TAccessFunctor, typename TOrderFunctor, unsigned Q>
213+
typename TText, unsigned Q>
202214
inline void
203-
radixSort(std::vector<std::tuple<TSAValue*, TSAValue*, TSize> > & stack,
204-
RadixSortContext_<TAccessFunctor, TOrderFunctor, TSize, Q> & context,
205-
std::tuple<TSAValue*, TSAValue*, TSize> const & item)
215+
_radixSort(std::vector<std::tuple<TSAValue*, TSAValue*, TSize> > & stack,
216+
RadixSortContext_<TSAValue, TText, TSize, Q> & context,
217+
std::tuple<TSAValue*, TSAValue*, TSize> const & item)
206218
{
207-
typedef RadixSortContext_<TAccessFunctor, TOrderFunctor, TSize, Q> TContext;
219+
typedef RadixSortContext_<TSAValue, TText, TSize, Q> TContext;
208220
typedef typename TContext::TOrdValue TOrdValue;
209-
static_assert(std::is_same<TSAValue, typename TContext::TSAValue>::value, "TSAValue mismatch!");
210221

211222
clear(context);
212223

@@ -261,85 +272,121 @@ radixSort(std::vector<std::tuple<TSAValue*, TSAValue*, TSize> > & stack,
261272
}
262273

263274
// ----------------------------------------------------------------------------
264-
// Function inplaceFullRadixSort() [default]
275+
// Function _radixSortWrapper()
265276
// ----------------------------------------------------------------------------
277+
// switch to quicksort if the interval is sufficiently small
266278

267279
//TODO: play with this value
280+
#ifndef _RADIX_SORT_SWITCH_TO_QUICKSORT_AT
268281
#define _RADIX_SORT_SWITCH_TO_QUICKSORT_AT 100
282+
#endif
283+
284+
template <typename TSAValue, typename TSize,
285+
typename TText, unsigned Q>
286+
inline void
287+
_radixSortWrapper(std::vector<std::tuple<TSAValue*, TSAValue*, TSize> > & stack,
288+
RadixSortContext_<TSAValue, TText, TSize, Q> & context,
289+
std::tuple<TSAValue*, TSAValue*, TSize> const & i)
290+
{
291+
if (std::get<1>(i) - std::get<0>(i) < _RADIX_SORT_SWITCH_TO_QUICKSORT_AT)
292+
std::sort(std::get<0>(i), std::get<1>(i), SuffixLess_<TSAValue, TText const>(context.text, std::get<2>(i)));
293+
else if (std::get<1>(i) - std::get<0>(i) >= 2)
294+
_radixSort(stack, context, i);
295+
}
296+
297+
// ----------------------------------------------------------------------------
298+
// Function inplaceFullRadixSort() [default]
299+
// ----------------------------------------------------------------------------
269300

270301
#ifdef _OPENMP
271-
#define N_THREADS omp_get_max_threads()
272-
#define I_THREAD omp_get_thread_num()
302+
#define N_THREADS omp_get_max_threads()
303+
#define I_THREAD omp_get_thread_num()
304+
#define MIN_BUCKETS 512
273305
#else
274-
#define N_THREADS 1
275-
#define I_THREAD 0
306+
#define N_THREADS 1
307+
#define I_THREAD 0
308+
#define MIN_BUCKETS 100 // for somewhat decent progress reporting
276309
#endif
277310

278311
// TODO: serial version
279312
// TODO: possibly do multiple runs of "secondStep" if alphabet size to small
280313
// TODO: possibly quicksort directly on buckets in third steps, if buckets have been made small enough
314+
// TODO: double-check the effects of the new "secondStep"
281315

282316
template <typename TSA, typename TText, typename TLambda>
283-
void inPlaceRadixSort(TSA & sa, TText const & str, TLambda const & progressCallback = [] (unsigned) {})
317+
void inPlaceRadixSort(TSA & sa, TText const & text, TLambda const & progressCallback = [] (unsigned) {})
284318
{
285-
typedef typename Value<typename Concatenator<TText>::Type>::Type TAlphabet;
286-
typedef typename Value<TSA>::Type TSAValue;
287-
typedef typename Size<TText>::Type TSize;
288-
typedef std::tuple<TSAValue*, TSAValue*, TSize> TItem;
289-
typedef typename StringSetLimits<TText const>::Type TLimitsString; // "Nothing" for Strings
290-
291-
typedef RadixTextAccessor<TSAValue, TText> TAccessor;
292-
typedef _ZeroBucketComparator<TSAValue,TLimitsString> TZeroComp;
319+
typedef typename Value<typename Concatenator<TText>::Type>::Type TAlphabet;
320+
typedef typename Value<TSA>::Type TSAValue;
321+
typedef typename Size<TText>::Type TSize;
322+
typedef std::tuple<TSAValue*, TSAValue*, TSize> TItem;
293323

294324
static const unsigned SIGMA = static_cast<unsigned>(ValueSize<TAlphabet>::VALUE) + 1;
295325
SEQAN_ASSERT_LT_MSG(SIGMA, 1000u, "Attention: inplace radix sort is not suited for large alphabets");
296326

297-
typedef RadixSortContext_<TAccessor, TZeroComp, TSize, SIGMA> TContext;
327+
typedef RadixSortContext_<TSAValue, TText, TSize, SIGMA> TContext;
298328

299-
if (empty(sa)) return; // otherwise access sa[0] fails
329+
if (empty(sa))
330+
return; // otherwise access sa[0] fails
300331

301332
/* stacks */
302333
std::vector<TItem> firstStack;
303334
firstStack.reserve(SIGMA);
304335
std::vector<TItem> secondStack;
305-
secondStack.reserve(SIGMA*SIGMA);
336+
secondStack.reserve(1000);
306337
std::vector<std::vector<TItem>> lStack(N_THREADS); // one per thread
307338
// reduce memory allocations in threads by reserving space
308339
for (auto & stack : lStack)
309340
stack.reserve(length(sa) / 1000);
310341

311342
/* contexts */
312-
TContext firstSecondContext{TAccessor(str), TZeroComp(stringSetLimits(str))};
313-
std::vector<TContext> lContext(N_THREADS, TContext{TAccessor(str), TZeroComp(stringSetLimits(str))});
343+
TContext firstSecondContext{text};
344+
std::vector<TContext> lContext(N_THREADS, TContext{text});
314345

346+
// FIRST STEP
315347
// sort by the first character
316-
radixSort(firstStack, firstSecondContext, TItem(&sa[0], &sa[0]+length(sa), 0));
317-
348+
_radixSortWrapper(firstStack, firstSecondContext, TItem(&sa[0], &sa[0]+length(sa), 0));
318349
progressCallback(5); // 5% progress guess after first char
319350

320-
// sort by second character
321-
SEQAN_OMP_PRAGMA(parallel for schedule(dynamic))
322-
for (unsigned j = 0; j < length(firstStack); ++j)
351+
// SECOND STEP
352+
// sort by next n characters until the stack has reached a good size for distinct parallelization
353+
// NOTE that for small alphabets in combination with small texts, this step might sort the entire SA
354+
while (!firstStack.empty())
323355
{
324-
TItem & i = firstStack[j];
325-
if (std::get<1>(i) - std::get<0>(i) < _RADIX_SORT_SWITCH_TO_QUICKSORT_AT)
326-
std::sort(std::get<0>(i),
327-
std::get<1>(i),
328-
SuffixLess_<TSAValue, TText const>(str, std::get<2>(i)));
329-
else if (std::get<1>(i) - std::get<0>(i) >= 2)
330-
radixSort(lStack[I_THREAD], lContext[I_THREAD], i);
331-
}
356+
SEQAN_OMP_PRAGMA(parallel for schedule(dynamic))
357+
for (unsigned j = 0; j < length(firstStack); ++j)
358+
_radixSortWrapper(lStack[I_THREAD], lContext[I_THREAD], firstStack[j]);
332359

333-
// merge local stacks and clear for next round
334-
for (auto & stack : lStack)
335-
{
336-
secondStack.insert(secondStack.end(), stack.begin(), stack.end());
337-
stack.clear();
338-
}
360+
// merge local stacks and clear for next round or next step
361+
for (auto & stack : lStack)
362+
{
363+
secondStack.insert(secondStack.end(), stack.begin(), stack.end());
364+
stack.clear();
365+
}
339366

340-
progressCallback(10); // 10% progress guess after second char
367+
// sort the stack by interval size so that large intervals are moved to front
368+
// this improves parallelization of dynamic schedule
369+
SORT(secondStack.begin(), secondStack.end(),
370+
[] (TItem const & l, TItem const & r)
371+
{
372+
return (std::get<1>(l) - std::get<0>(l)) > (std::get<1>(r) - std::get<0>(r));
373+
});
374+
375+
// check if largest interval "fits" in one thread efficiently
376+
// this works independently of alphabet size and just depends on the data
377+
// MIN_BUCKETS check additionally guarantees a degree of granularity
378+
if ((uint64_t(std::get<1>(secondStack.front()) - std::get<0>(secondStack.front())) <= (length(sa) / N_THREADS))
379+
&& (secondStack.size() >= MIN_BUCKETS))
380+
break;
381+
382+
// switch buffers for next round
383+
firstStack.clear();
384+
std::swap(firstStack, secondStack);
385+
}
386+
progressCallback(10); // 10% progress guess after second step
341387

342-
// sort rest
388+
// THIRD STEP
389+
// sort the remaining intervals distinctly; here no locking and syncing is required anymore
343390
SEQAN_OMP_PRAGMA(parallel for schedule(dynamic))
344391
for (unsigned j = 0; j < secondStack.size(); ++j)
345392
{
@@ -349,21 +396,15 @@ void inPlaceRadixSort(TSA & sa, TText const & str, TLambda const & progressCallb
349396
{
350397
TItem i = lStack[I_THREAD].back();
351398
lStack[I_THREAD].pop_back();
352-
353-
if (std::get<1>(i) - std::get<0>(i) < _RADIX_SORT_SWITCH_TO_QUICKSORT_AT)
354-
std::sort(std::get<0>(i),
355-
std::get<1>(i),
356-
SuffixLess_<TSAValue, TText const>(str, std::get<2>(i)));
357-
else if (std::get<1>(i) - std::get<0>(i) >= 2)
358-
radixSort(lStack[I_THREAD], lContext[I_THREAD], i);
399+
_radixSortWrapper(lStack[I_THREAD], lContext[I_THREAD], i);
359400
}
360401

361402
// progressCallback must be thread safe and cope with smaller numbers after big numbers
362403
// remaining characters alloted 90% of total progress
363404
progressCallback(10 + (j * 90) / secondStack.size());
364405
}
365406

366-
progressCallback(100); // done
407+
// progressCallback(100); // done
367408
}
368409

369410
}

0 commit comments

Comments
 (0)