cpuRadixSorter.h

Go to the documentation of this file.

/*--------------------------------*- VM2D -*-----------------*---------------*\
| ##  ## ##   ##  ####  #####   |                            | Version 1.14   |
| ##  ## ### ### ##  ## ##  ##  |  VM2D: Vortex Method       | 2026/03/06     |
| ##  ## ## # ##    ##  ##  ##  |  for 2D Flow Simulation    *----------------*
|  ####  ##   ##   ##   ##  ##  |  Open Source Code                           |
|   ##   ##   ## ###### #####   |  https://www.github.com/vortexmethods/VM2D  |
|                                                                             |
| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova   |
*-----------------------------------------------------------------------------*
| File name: cpuRadixSorter.h                                                 |
| Info: Source code of VM2D                                                   |
|                                                                             |
| This file is part of VM2D.                                                  |
| VM2D is free software: you can redistribute it and/or modify it             |
| under the terms of the GNU General Public License as published by           |
| the Free Software Foundation, either version 3 of the License, or           |
| (at your option) any later version.                                         |
|                                                                             |
| VM2D is distributed in the hope that it will be useful, but WITHOUT         |
| ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       |
| FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License       |
| for more details.                                                           |
|                                                                             |
| You should have received a copy of the GNU General Public License           |
| along with VM2D.  If not, see <http://www.gnu.org/licenses/>.               |
\*---------------------------------------------------------------------------*/
 
#ifndef CPURADIXSORT_H
#define CPURADIXSORT_H
 
#include <iostream>
#include <vector>
#include <algorithm>
#include <type_traits>
#include <omp.h>
#include <string>
 
namespace VM2D
{
 
    // SFINAE äëÿ ïðîâåðêè òðèâèàëüíîé êîïèðóåìîñòè (îïòèìèçàöèÿ)
    template<typename T>
    using enable_if_trivially_copyable = typename std::enable_if<std::is_trivially_copyable<T>::value>::type;
 
    template<typename T>
    using enable_if_not_trivially_copyable = typename std::enable_if<!std::is_trivially_copyable<T>::value>::type;
 
    // Îñíîâíîé øàáëîííûé êëàññ
    template<typename ValueType>
    class OptimizedRadixSorter {
    private:
        // Ïåðåèñïîëüçóåìûå áóôåðû
        std::vector<int> keyBuffer;
        std::vector<ValueType> valueBuffer;
 
        std::vector<int> globalHist;
        std::vector<std::vector<int>> histThreads;
        std::vector<std::vector<int>> localStarts;
        std::vector<std::vector<int>> localCounters;
 
        int numThreads;
        size_t currentCapacity;
        bool initialized;
 
        constexpr static int RADIX = 256;
        constexpr static int MASK = 0xFF;
        constexpr static int PASSES = 4; // äëÿ 32-bit int
 
        void initializeThreadData(int threads) {
            numThreads = threads;
 
            // Âûäåëÿåì ïàìÿòü äëÿ ãèñòîãðàìì ïîòîêîâ
            histThreads.resize(numThreads);
            for (int t = 0; t < numThreads; ++t) {
                histThreads[t].assign(RADIX, 0);
            }
 
            // Âûäåëÿåì ïàìÿòü äëÿ ëîêàëüíûõ ñòàðòîâûõ ïîçèöèé
            localStarts.resize(numThreads);
            localCounters.resize(numThreads);
            for (int t = 0; t < numThreads; ++t) {
                localStarts[t].assign(RADIX, 0);
                localCounters[t].assign(RADIX, 0);
            }
 
            // Ãëîáàëüíàÿ ãèñòîãðàììà
            globalHist.assign(RADIX, 0);
 
            initialized = true;
        }
 
        void ensureCapacity(size_t n) {
            if (n > currentCapacity) {
                // Óâåëè÷èâàåì áóôåðû ñ çàïàñîì 10% äëÿ áóäóùèõ ñîðòèðîâîê
                keyBuffer.resize(n + n / 10);
                valueBuffer.resize(n + n / 10);
                currentCapacity = keyBuffer.size();
            }
        }
 
        void resetHistograms() {
            // Áûñòðûé ñáðîñ ãèñòîãðàìì (ïåðåèñïîëüçóåì ïàìÿòü)
            for (int t = 0; t < numThreads; ++t) {
                std::fill(histThreads[t].begin(), histThreads[t].end(), 0);
            }
            std::fill(globalHist.begin(), globalHist.end(), 0);
        }
 
        // Îïòèìèçèðîâàííîå êîïèðîâàíèå äëÿ òðèâèàëüíî êîïèðóåìûõ òèïîâ
        template<typename T = ValueType>
        void copyValues(T* dest, const T* src, size_t n, enable_if_trivially_copyable<T>* = nullptr) {
            std::memcpy(dest, src, n * sizeof(T));
        }
 
        // Êîïèðîâàíèå äëÿ íåòðèâèàëüíî êîïèðóåìûõ òèïîâ
        template<typename T = ValueType>
        void copyValues(T* dest, const T* src, size_t n, enable_if_not_trivially_copyable<T>* = nullptr) {
            for (size_t i = 0; i < n; ++i) {
                dest[i] = src[i];
            }
        }
 
    public:
        OptimizedRadixSorter() : numThreads(0), currentCapacity(0), initialized(false) {
#pragma omp parallel
            {
#pragma omp single
                numThreads = omp_get_num_threads();
            }
            if (numThreads == 0) numThreads = 1;
 
            initializeThreadData(numThreads);
        }
 
        // Êîíñòðóêòîð ñ ïðåäâàðèòåëüíûì âûäåëåíèåì ïàìÿòè
        explicit OptimizedRadixSorter(size_t maxSize) : numThreads(0), currentCapacity(0), initialized(false) {
#pragma omp parallel
            {
#pragma omp single
                numThreads = omp_get_num_threads();
            }
            if (numThreads == 0) numThreads = 1;
 
            initializeThreadData(numThreads);
            ensureCapacity(maxSize);
        }
 
        // Îñíîâíîé ìåòîä ñîðòèðîâêè
        void sort(int* keys, ValueType* values, size_t n) {
            if (n < 2) return;
 
            ensureCapacity(n);
            resetHistograms();
 
            for (int pass = 0; pass < PASSES; ++pass) {
                int shift = pass * 8;
 
                // 1. Ïàðàëëåëüíûé ïîäñ÷¸ò ãèñòîãðàìì
#pragma omp parallel for schedule(static)
                for (int i = 0; i < n; ++i) {
                    int tid = omp_get_thread_num();
                    int bucket = (keys[i] >> shift) & MASK;
                    histThreads[tid][bucket]++;
                }
 
                // 2. Îáúåäèíåíèå ãèñòîãðàìì
                for (int t = 0; t < numThreads; ++t) {
                    const auto& hist = histThreads[t];
                    auto& global = globalHist;
                    for (int b = 0; b < RADIX; ++b) {
                        global[b] += hist[b];
                    }
                }
 
                // 3. Âû÷èñëåíèå ïðåôèêñíûõ ñóìì
                std::vector<int> prefix(RADIX);
                int sum = 0;
                for (int b = 0; b < RADIX; ++b) {
                    prefix[b] = sum;
                    sum += globalHist[b];
                }
 
                // 4. Ðàñ÷¸ò ñòàðòîâûõ ïîçèöèé äëÿ êàæäîãî ïîòîêà
                std::copy(prefix.begin(), prefix.end(), localStarts[0].begin());
 
                for (int t = 1; t < numThreads; ++t) {
                    auto& start = localStarts[t];
                    const auto& prevStart = localStarts[t - 1];
                    const auto& prevHist = histThreads[t - 1];
 
                    for (int b = 0; b < RADIX; ++b) {
                        start[b] = prevStart[b] + prevHist[b];
                    }
                }
 
                // 5. Êîïèðóåì ñòàðòîâûå ïîçèöèè â ñ÷¸ò÷èêè
                for (int t = 0; t < numThreads; ++t) {
                    std::copy(localStarts[t].begin(), localStarts[t].end(),
                        localCounters[t].begin());
                }
 
                // 6. Ïàðàëëåëüíàÿ ðàçíîñêà ýëåìåíòîâ
#pragma omp parallel for schedule(static)
                for (int i = 0; i < n; ++i) {
                    int tid = omp_get_thread_num();
                    int bucket = (keys[i] >> shift) & MASK;
                    int pos = localCounters[tid][bucket]++;
 
                    keyBuffer[pos] = keys[i];
                    valueBuffer[pos] = values[i];
                }
 
                // 7. Êîïèðîâàíèå îáðàòíî â èñõîäíûå ìàññèâû
#pragma omp parallel for schedule(static)
                for (int i = 0; i < n; ++i) {
                    keys[i] = keyBuffer[i];
                    values[i] = valueBuffer[i];
                }
 
                // 8. Ñáðîñ ãèñòîãðàìì
                for (int t = 0; t < numThreads; ++t) {
                    std::fill(histThreads[t].begin(), histThreads[t].end(), 0);
                }
                std::fill(globalHist.begin(), globalHist.end(), 0);
            }
        }
 
        // Ïåðåãðóæåííûé ìåòîä äëÿ ðàáîòû ñ std::vector
        void sort(std::vector<int>& keys, std::vector<ValueType>& values) {
            if (keys.size() != values.size()) {
                throw std::runtime_error("Key and value arrays must have same size");
            }
            sort(keys.data(), values.data(), keys.size());
        }
 
        // Ìåòîä äëÿ ñîðòèðîâêè ñ óêàçàíèåì äèàïàçîíà
        void sort(int* keys, ValueType* values, size_t start, size_t end) {
            if (start >= end) return;
            sort(keys + start, values + start, end - start);
        }
 
        void releaseMemory() {
            std::vector<int>().swap(keyBuffer);
            std::vector<ValueType>().swap(valueBuffer);
            std::vector<int>().swap(globalHist);
            std::vector<std::vector<int>>().swap(histThreads);
            std::vector<std::vector<int>>().swap(localStarts);
            std::vector<std::vector<int>>().swap(localCounters);
            currentCapacity = 0;
            initialized = false;
        }
 
        size_t getMemoryUsage() const {
            size_t total = keyBuffer.capacity() * sizeof(int);
            total += valueBuffer.capacity() * sizeof(ValueType);
            total += globalHist.capacity() * sizeof(int);
            for (const auto& hist : histThreads) {
                total += hist.capacity() * sizeof(int);
            }
            for (const auto& start : localStarts) {
                total += start.capacity() * sizeof(int);
            }
            for (const auto& counter : localCounters) {
                total += counter.capacity() * sizeof(int);
            }
            return total;
        }
 
        // Ïîëó÷èòü ÷èñëî ïîòîêîâ
        int getNumThreads() const { return numThreads; }
 
        // Óñòàíîâèòü ÷èñëî ïîòîêîâ (äîëæíî áûòü âûçâàíî ÄÎ ïåðâîé ñîðòèðîâêè)
        void setNumThreads(int threads) {
            if (!initialized) {
                numThreads = threads;
                initializeThreadData(numThreads);
            }
        }
    };
 
    // ×àñòè÷íàÿ ñïåöèàëèçàöèÿ äëÿ óêàçàòåëåé (åñëè íóæíî õðàíèòü óêàçàòåëè)
    template<typename ValueType>
    class OptimizedRadixSorter<ValueType*> {
    private:
        std::vector<int> keyBuffer;
        std::vector<ValueType*> valueBuffer;
        // ... îñòàëüíàÿ ðåàëèçàöèÿ àíàëîãè÷íà
    };
 
    // Ïðèìåðû èñïîëüçîâàíèÿ
    struct ComplexData {
        double x;
        double y;
        int id;
 
        ComplexData() : x(0), y(0), id(0) {}
        ComplexData(double _x, double _y, int _id) : x(_x), y(_y), id(_id) {}
 
        bool operator==(const ComplexData& other) const {
            return x == other.x && y == other.y && id == other.id;
        }
    };
 
}//VM2D
 
#endif