代码

#pragma GCC target("avx2")
#include <immintrin.h>
#include <pthread.h>

#define N (1U << 27)
#define BUCKETS 256U
#define BUFFER_SIZE 16U
#ifndef THREADS
#define THREADS 4U
#endif

typedef struct __attribute__((aligned(4096))) {
    unsigned high_count[BUCKETS];
    unsigned position[BUCKETS];
    unsigned local_count[3][BUCKETS];
    unsigned fill[BUCKETS];
    unsigned buffer[BUCKETS][BUFFER_SIZE] __attribute__((aligned(64)));
    unsigned id;
} Worker;

static unsigned work[N] __attribute__((aligned(4096)));
static unsigned starts[BUCKETS + 1];
static Worker workers[THREADS];
static pthread_barrier_t barrier;
static unsigned next_bucket;
static unsigned *input_array;

static inline void prefix_from(unsigned *count, unsigned base) {
    unsigned sum = base;
    for (unsigned i = 0; i < BUCKETS; ++i) {
        unsigned value = count[i];
        count[i] = sum;
        sum += value;
    }
}

static __attribute__((always_inline)) inline void flush_nt(
        unsigned *destination, const unsigned *source) {
    _mm256_stream_si256((__m256i *)(destination + 0),
                        _mm256_load_si256((const __m256i *)(source + 0)));
    _mm256_stream_si256((__m256i *)(destination + 8),
                        _mm256_load_si256((const __m256i *)(source + 8)));
}

static __attribute__((always_inline)) inline void radix_pass_range(
        Worker *worker, const unsigned *source, unsigned *destination,
        unsigned *position, unsigned shift, unsigned begin, unsigned end) {
    unsigned *fill = worker->fill;
    for (unsigned i = begin; i < end; ++i) {
        unsigned value = source[i];
        unsigned digit = (value >> shift) & 255U;
        unsigned output = position[digit];
        if (output & 7U) {
            destination[output] = value;
            position[digit] = output + 1;
        } else {
            unsigned slot = fill[digit];
            worker->buffer[digit][slot] = value;
            if (slot + 1 == BUFFER_SIZE) {
                flush_nt(destination + output, worker->buffer[digit]);
                position[digit] = output + BUFFER_SIZE;
                fill[digit] = 0;
            } else {
                fill[digit] = slot + 1;
            }
        }
    }
    for (unsigned digit = 0; digit < BUCKETS; ++digit) {
        unsigned amount = fill[digit];
        unsigned output = position[digit];
        for (unsigned j = 0; j < amount; ++j)
            destination[output + j] = worker->buffer[digit][j];
        position[digit] = output + amount;
        fill[digit] = 0;
    }
    _mm_sfence();
}

static void *worker_main(void *argument) {
    Worker *worker = argument;
    unsigned id = worker->id;
    unsigned begin = (unsigned)(((unsigned long long)N * id) / THREADS);
    unsigned end = (unsigned)(((unsigned long long)N * (id + 1)) / THREADS);

    for (unsigned i = begin; i < end; ++i)
        ++worker->high_count[input_array[i] >> 24];
    pthread_barrier_wait(&barrier);

    if (id == 0) {
        unsigned sum = 0;
        for (unsigned digit = 0; digit < BUCKETS; ++digit) {
            starts[digit] = sum;
            for (unsigned thread = 0; thread < THREADS; ++thread) {
                workers[thread].position[digit] = sum;
                sum += workers[thread].high_count[digit];
            }
        }
        starts[BUCKETS] = N;
    }
    pthread_barrier_wait(&barrier);

    radix_pass_range(worker, input_array, work, worker->position, 24,
                     begin, end);
    pthread_barrier_wait(&barrier);

    if (id == 0) next_bucket = 0;
    pthread_barrier_wait(&barrier);
    for (;;) {
        unsigned bucket = __atomic_fetch_add(&next_bucket, 1, __ATOMIC_RELAXED);
        if (bucket >= BUCKETS) break;
        begin = starts[bucket];
        end = starts[bucket + 1];
        __builtin_memset(worker->local_count, 0,
                         sizeof(worker->local_count));
        for (unsigned i = begin; i < end; ++i) {
            unsigned value = work[i];
            ++worker->local_count[0][(unsigned char)value];
            ++worker->local_count[1][(unsigned char)(value >> 8)];
            ++worker->local_count[2][(unsigned char)(value >> 16)];
        }
        prefix_from(worker->local_count[0], begin);
        prefix_from(worker->local_count[1], begin);
        prefix_from(worker->local_count[2], begin);
        radix_pass_range(worker, work, input_array, worker->local_count[0],
                         0, begin, end);
        radix_pass_range(worker, input_array, work, worker->local_count[1],
                         8, begin, end);
        radix_pass_range(worker, work, input_array, worker->local_count[2],
                         16, begin, end);
    }
    _mm256_zeroupper();
    return 0;
}

void sort(unsigned *a, int n) {
    (void)n;
    pthread_t thread[THREADS - 1];
    input_array = a;
    pthread_barrier_init(&barrier, 0, THREADS);
    for (unsigned i = 0; i < THREADS; ++i) workers[i].id = i;
    for (unsigned i = 1; i < THREADS; ++i)
        pthread_create(&thread[i - 1], 0, worker_main, &workers[i]);
    worker_main(&workers[0]);
    for (unsigned i = 1; i < THREADS; ++i)
        pthread_join(thread[i - 1], 0);
    pthread_barrier_destroy(&barrier);
}

评测结果