include/comphelper/parallelsort.hxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  */
   9
  10 #ifndef INCLUDED_COMPHELPER_PARALLELSORT_HXX
  11 #define INCLUDED_COMPHELPER_PARALLELSORT_HXX
  12
  13 #include <comphelper/threadpool.hxx>
  14 #include <tools/cpuid.hxx>
  15
  16 #include <memory>
  17 #include <iterator>
  18 #include <thread>
  19 #include <algorithm>
  20 #include <cmath>
  21 #include <random>
  22 #include <functional>
  23 #include <iostream>
  24 #include <chrono>
  25
  26 namespace comphelper
  27 {
  28 const size_t nThreadCountGlobal = std::thread::hardware_concurrency();
  29 const bool bHyperThreadingActive = cpuid::hasHyperThreading();
  30 static comphelper::ThreadPool& rTPool(comphelper::ThreadPool::getSharedOptimalPool());
  31
  32 static thread_local std::mt19937 aGenerator{ std::random_device{}() };
  33
  34 #define PARALLELSORT_ENABLEPZ 0
  35
  36 namespace
  37 {
  38 class ProfileZone
  39 {
  40 public:
  41 #if PARALLELSORT_ENABLEPZ
  42     ProfileZone(const char* pTag)
  43         : maTag(pTag)
  44         , maStart(std::chrono::steady_clock::now())
  45         , mbFinished(false)
  46     {
  47     }
  48
  49     ~ProfileZone()
  50     {
  51         if (!mbFinished)
  52             showTimeElapsed();
  53     }
  54
  55     void stop()
  56     {
  57         showTimeElapsed();
  58         mbFinished = true;
  59     }
  60 #else
  61     ProfileZone(const char* /*pTag*/)
  62         : mbDummy(true)
  63     {
  64     }
  65
  66     void stop()
  67     {
  68         // Avoid loplugin:staticmethods, loplugin:staticaccess errors
  69         (void)mbDummy;
  70     }
  71 #endif
  72
  73 private:
  74 #if PARALLELSORT_ENABLEPZ
  75
  76     void showTimeElapsed()
  77     {
  78         auto end = std::chrono::steady_clock::now();
  79         size_t elapsed
  80             = std::chrono::duration_cast<std::chrono::milliseconds>(end - maStart).count();
  81         std::cout << maTag << " : " << elapsed << " ms" << std::endl << std::flush;
  82     }
  83
  84     std::string maTag;
  85     std::chrono::steady_clock::time_point maStart;
  86     bool mbFinished;
  87 #else
  88     bool mbDummy;
  89
  90 #endif
  91 };
  92
  93 class ParallelRunner
  94 {
  95     class Executor final : public comphelper::ThreadTask
  96     {
  97     public:
  98         Executor(const std::shared_ptr<comphelper::ThreadTaskTag>& rTag,
  99                  std::function<void()> aFunc)
 100             : comphelper::ThreadTask(rTag)
 101             , maFunc(std::move(aFunc))
 102         {
 103         }
 104
 105         virtual void doWork() override { maFunc(); }
 106
 107     private:
 108         const std::function<void()> maFunc;
 109     };
 110
 111 public:
 112     ParallelRunner() { maTag = comphelper::ThreadPool::createThreadTaskTag(); }
 113
 114     void enqueue(std::function<void()> aFunc)
 115     {
 116         rTPool.pushTask(std::make_unique<Executor>(maTag, aFunc));
 117     }
 118
 119     void wait() { rTPool.waitUntilDone(maTag, false); }
 120
 121 private:
 122     std::shared_ptr<comphelper::ThreadTaskTag> maTag;
 123 };
 124
 125 constexpr size_t nMaxTreeArraySize = 64;
 126
 127 size_t lcl_tree_array_size(size_t nNum)
 128 {
 129     size_t nPow2;
 130     for (nPow2 = 1; nPow2 <= nNum; nPow2 <<= 1)
 131         ;
 132     return std::clamp((nPow2 >> 1), size_t(1), nMaxTreeArraySize);
 133 }
 134
 135 template <class RandItr> struct Sampler
 136 {
 137     using ValueType = typename std::iterator_traits<RandItr>::value_type;
 138
 139     static void sample(RandItr aBegin, RandItr aEnd, ValueType* pSamples, size_t nSamples,
 140                        size_t /*nParallelism*/)
 141     {
 142         ProfileZone aZone("\tsample()");
 143         assert(aBegin <= aEnd);
 144         size_t nLen = static_cast<std::size_t>(aEnd - aBegin);
 145         assert(std::mt19937::max() >= nLen);
 146
 147         for (size_t nIdx = 0; nIdx < nSamples; ++nIdx)
 148         {
 149             size_t nSel = aGenerator() % nLen--;
 150             using namespace std;
 151             swap(*(aBegin + nSel), *(aBegin + nLen));
 152             pSamples[nIdx] = *(aBegin + nLen);
 153         }
 154     }
 155 };
 156
 157 template <class RandItr, class Compare> class Binner
 158 {
 159     using ValueType = typename std::iterator_traits<RandItr>::value_type;
 160
 161     const size_t mnTreeArraySize;
 162     const size_t mnDividers;
 163     constexpr static size_t mnMaxStaticSize = 1024 * 50;
 164     uint8_t maLabels[mnMaxStaticSize];
 165     ValueType maDividers[nMaxTreeArraySize];
 166     std::unique_ptr<uint8_t[]> pLabels;
 167     size_t maSepBinEnds[nMaxTreeArraySize * nMaxTreeArraySize];
 168     bool mbThreaded;
 169
 170 public:
 171     size_t maBinEnds[nMaxTreeArraySize];
 172
 173     Binner(const ValueType* pSamples, size_t nSamples, size_t nBins, bool bThreaded)
 174         : mnTreeArraySize(lcl_tree_array_size(nBins))
 175         , mnDividers(mnTreeArraySize - 1)
 176         , mbThreaded(bThreaded)
 177     {
 178         assert((nSamples % mnTreeArraySize) == 0);
 179         assert(mnTreeArraySize <= nMaxTreeArraySize);
 180         std::fill(maBinEnds, maBinEnds + mnTreeArraySize, 0);
 181         std::fill(maSepBinEnds, maSepBinEnds + mnTreeArraySize * mnTreeArraySize, 0);
 182         fillTreeArray(1, pSamples, pSamples + nSamples);
 183     }
 184
 185     void fillTreeArray(size_t nPos, const ValueType* pLow, const ValueType* pHigh)
 186     {
 187         assert(pLow <= pHigh);
 188         const ValueType* pMid = pLow + (pHigh - pLow) / 2;
 189         maDividers[nPos] = *pMid;
 190
 191         if (2 * nPos < mnDividers) // So that 2*nPos < mnTreeArraySize
 192         {
 193             fillTreeArray(2 * nPos, pLow, pMid);
 194             fillTreeArray(2 * nPos + 1, pMid + 1, pHigh);
 195         }
 196     }
 197
 198     constexpr inline size_t findBin(const ValueType& rVal, Compare& aComp)
 199     {
 200         size_t nIdx = 1;
 201         while (nIdx <= mnDividers)
 202             nIdx = ((nIdx << 1) + aComp(maDividers[nIdx], rVal));
 203         return (nIdx - mnTreeArraySize);
 204     }
 205
 206     void label(const RandItr aBegin, const RandItr aEnd, Compare& aComp)
 207     {
 208         ProfileZone aZoneSetup("\tlabel():setup");
 209         size_t nLen = static_cast<std::size_t>(aEnd - aBegin);
 210         if (nLen > mnMaxStaticSize)
 211             pLabels = std::make_unique<uint8_t[]>(nLen);
 212         uint8_t* pLabelsRaw = (nLen > mnMaxStaticSize) ? pLabels.get() : maLabels;
 213         aZoneSetup.stop();
 214         ProfileZone aZoneFindBins("\tFindBins()");
 215         if (mbThreaded)
 216         {
 217             ParallelRunner aPRunner;
 218             const size_t nBins = mnTreeArraySize;
 219             for (size_t nTIdx = 0; nTIdx < nBins; ++nTIdx)
 220             {
 221                 aPRunner.enqueue([this, nTIdx, nBins, nLen, aBegin, pLabelsRaw, &aComp] {
 222                     ProfileZone aZoneIn("\t\tFindBinsThreaded()");
 223                     size_t nBinEndsStartIdx = nTIdx * mnTreeArraySize;
 224                     size_t* pBinEnds = maSepBinEnds + nBinEndsStartIdx;
 225                     size_t aBinEndsF[nMaxTreeArraySize] = { 0 };
 226                     for (size_t nIdx = nTIdx; nIdx < nLen; nIdx += nBins)
 227                     {
 228                         size_t nBinIdx = findBin(*(aBegin + nIdx), aComp);
 229                         pLabelsRaw[nIdx] = static_cast<uint8_t>(nBinIdx);
 230                         ++aBinEndsF[nBinIdx];
 231                     }
 232
 233                     for (size_t nIdx = 0; nIdx < mnTreeArraySize; ++nIdx)
 234                         pBinEnds[nIdx] = aBinEndsF[nIdx];
 235                 });
 236             }
 237
 238             aPRunner.wait();
 239
 240             // Populate maBinEnds from maSepBinEnds
 241             for (size_t nTIdx = 0; nTIdx < mnTreeArraySize; ++nTIdx)
 242             {
 243                 for (size_t nSepIdx = 0; nSepIdx < mnTreeArraySize; ++nSepIdx)
 244                     maBinEnds[nTIdx] += maSepBinEnds[nSepIdx * mnTreeArraySize + nTIdx];
 245             }
 246         }
 247         else
 248         {
 249             uint8_t* pLabel = pLabelsRaw;
 250             for (RandItr aItr = aBegin; aItr != aEnd; ++aItr)
 251             {
 252                 size_t nBinIdx = findBin(*aItr, aComp);
 253                 *pLabel++ = nBinIdx;
 254                 ++maBinEnds[nBinIdx];
 255             }
 256         }
 257
 258         aZoneFindBins.stop();
 259
 260         size_t nSum = 0;
 261         // Store each bin's starting position in maBinEnds array for now.
 262         for (size_t nIdx = 0; nIdx < mnTreeArraySize; ++nIdx)
 263         {
 264             size_t nSize = maBinEnds[nIdx];
 265             maBinEnds[nIdx] = nSum;
 266             nSum += nSize;
 267         }
 268
 269         // Now maBinEnds has end positions of each bin.
 270     }
 271
 272     void bin(const RandItr aBegin, const RandItr aEnd, ValueType* pOut)
 273     {
 274         ProfileZone aZone("\tbin()");
 275         const size_t nLen = static_cast<std::size_t>(aEnd - aBegin);
 276         uint8_t* pLabelsRaw = (nLen > mnMaxStaticSize) ? pLabels.get() : maLabels;
 277         size_t nIdx;
 278         for (nIdx = 0; nIdx < nLen; ++nIdx)
 279         {
 280             pOut[maBinEnds[pLabelsRaw[nIdx]]++] = *(aBegin + nIdx);
 281         }
 282     }
 283 };
 284
 285 template <class RandItr, class Compare = std::less<>>
 286 void s3sort(const RandItr aBegin, const RandItr aEnd, Compare aComp = Compare(),
 287             bool bThreaded = true)
 288 {
 289     static size_t nThreadCount = nThreadCountGlobal;
 290
 291     constexpr size_t nBaseCaseSize = 1024;
 292     const std::size_t nLen = static_cast<std::size_t>(aEnd - aBegin);
 293     if (nLen < nBaseCaseSize)
 294     {
 295         std::stable_sort(aBegin, aEnd, aComp);
 296         return;
 297     }
 298
 299     using ValueType = typename std::iterator_traits<RandItr>::value_type;
 300     auto pOut = std::make_unique<ValueType[]>(nLen);
 301
 302     const size_t nBins = lcl_tree_array_size(nThreadCount);
 303     assert(nBins >= 1);
 304     const size_t nOverSamplingFactor = std::max(1.0, std::sqrt(static_cast<double>(nLen) / 64));
 305     const size_t nSamples = nOverSamplingFactor * nBins;
 306     auto aSamples = std::make_unique<ValueType[]>(nSamples);
 307     ProfileZone aZoneSampleAnsSort("SampleAndSort");
 308     // Select samples and sort them
 309     Sampler<RandItr>::sample(aBegin, aEnd, aSamples.get(), nSamples, nBins);
 310     std::sort(aSamples.get(), aSamples.get() + nSamples, aComp);
 311     aZoneSampleAnsSort.stop();
 312
 313     if (!aComp(aSamples[0], aSamples[nSamples - 1]))
 314     {
 315         // All samples are equal, fallback to standard sort.
 316         std::sort(aBegin, aEnd, aComp);
 317         return;
 318     }
 319
 320     ProfileZone aZoneBinner("Binner");
 321     // Create and populate bins using pOut from input iterators.
 322     Binner<RandItr, Compare> aBinner(aSamples.get(), nSamples, nBins, bThreaded);
 323     aBinner.label(aBegin, aEnd, aComp);
 324     aBinner.bin(aBegin, aEnd, pOut.get());
 325     aZoneBinner.stop();
 326
 327     ProfileZone aZoneSortBins("SortBins");
 328     ValueType* pOutRaw = pOut.get();
 329     if (bThreaded)
 330     {
 331         ParallelRunner aPRunner;
 332         // Sort the bins separately.
 333         for (size_t nBinIdx = 0, nBinStart = 0; nBinIdx < nBins; ++nBinIdx)
 334         {
 335             size_t nBinEnd = aBinner.maBinEnds[nBinIdx];
 336             aPRunner.enqueue([pOutRaw, nBinStart, nBinEnd, &aComp] {
 337                 std::sort(pOutRaw + nBinStart, pOutRaw + nBinEnd, aComp);
 338             });
 339
 340             nBinStart = nBinEnd;
 341         }
 342
 343         aPRunner.wait();
 344     }
 345     else
 346     {
 347         for (size_t nBinIdx = 0, nBinStart = 0; nBinIdx < nBins; ++nBinIdx)
 348         {
 349             auto nBinEnd = aBinner.maBinEnds[nBinIdx];
 350             std::sort(pOutRaw + nBinStart, pOutRaw + nBinEnd, aComp);
 351             nBinStart = nBinEnd;
 352         }
 353     }
 354
 355     aZoneSortBins.stop();
 356
 357     // Move the sorted array to the array specified by input iterators.
 358     std::move(pOutRaw, pOutRaw + nLen, aBegin);
 359 }
 360
 361 } // anonymous namespace
 362
 363 template <class RandItr, class Compare = std::less<>>
 364 void parallelSort(const RandItr aBegin, const RandItr aEnd, Compare aComp = Compare())
 365 {
 366     assert(aBegin <= aEnd);
 367     s3sort(aBegin, aEnd, aComp);
 368 }
 369
 370 } // namespace comphelper
 371
 372 #endif // INCLUDED_COMPHELPER_PARALLELSORT_HXX
 373
 374 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */