Конкурс: самый быстрый способ сортировки большого массива данных, распределенных по Гауссу

Учитывая интерес к этому вопросу , я подумал, что было бы интересно сделать ответы более объективными и количественными, предложив конкурс.

Идея проста: я сгенерировал двоичный файл, содержащий 50 миллионов распределенных по Гауссу двойных чисел (в среднем: 0, stdev 1). Цель состоит в том, чтобы сделать программу, которая будет сортировать их в памяти как можно быстрее. Очень простая эталонная реализация в python требует 1m4s. Как низко мы можем пойти?

Правила таковы: ответьте с помощью программы, которая открывает файл "gaussian.dat" и сортирует числа в памяти (не нужно их выводить), а также инструкции по сборке и запуску программы. Программа должна быть в состоянии работать на моем компьютере с Arch Linux (это означает, что вы можете использовать любой язык программирования или библиотеку, которую легко установить в этой системе).

Программа должна быть в достаточной степени читабельной, чтобы я мог убедиться, что она безопасна для запуска (нет решения только для ассемблера, пожалуйста!).

Я буду запускать ответы на своей машине (четырехъядерный процессор, 4 гигабайта оперативной памяти). Самое быстрое решение получит принятый ответ и награду в 100 баллов :)

Программа для генерации чисел:

#!/usr/bin/env python
import random
from array import array
from sys import argv
count=int(argv[1])
a=array('d',(random.gauss(0,1) for x in xrange(count)))
f=open("gaussian.dat","wb")
a.tofile(f)

Простая справочная реализация:

#!/usr/bin/env python
from array import array
from sys import argv
count=int(argv[1])
a=array('d')
a.fromfile(open("gaussian.dat"),count)
print "sorting..."
b=sorted(a)

РЕДАКТИРОВАТЬ: только 4 ГБ оперативной памяти, извините

РЕДАКТИРОВАТЬ # 2: Обратите внимание, что смысл конкурса заключается в том, чтобы увидеть, можем ли мы использовать предварительную информацию о данных . это не должно быть отличным совпадением между различными реализациями языка программирования!

Вот решение в C ++, которое сначала разбивает числа на сегменты с одинаковым ожидаемым количеством элементов, а затем сортирует каждый сегмент отдельно. Он предварительно вычисляет таблицу накопительной функции распределения на основе некоторых формул из Википедии, а затем интерполирует значения из этой таблицы, чтобы получить быстрое приближение.

Несколько шагов выполняются в нескольких потоках, чтобы использовать четыре ядра.

#include <cstdlib>
#include <math.h>
#include <stdio.h>
#include <algorithm>

#include <tbb/parallel_for.h>

using namespace std;

typedef unsigned long long ull;

double signum(double x) {
    return (x<0) ? -1 : (x>0) ? 1 : 0;
}

const double fourOverPI = 4 / M_PI;

double erf(double x) {
    double a = 0.147;
    double x2 = x*x;
    double ax2 = a*x2;
    double f1 = -x2 * (fourOverPI + ax2) / (1 + ax2);
    double s1 = sqrt(1 - exp(f1));
    return signum(x) * s1;
}

const double sqrt2 = sqrt(2);

double cdf(double x) {
    return 0.5 + erf(x / sqrt2) / 2;
}

const int cdfTableSize = 200;
const double cdfTableLimit = 5;
double* computeCdfTable(int size) {
    double* res = new double[size];
    for (int i = 0; i < size; ++i) {
        res[i] = cdf(cdfTableLimit * i / (size - 1));
    }
    return res;
}
const double* const cdfTable = computeCdfTable(cdfTableSize);

double cdfApprox(double x) {
    bool negative = (x < 0);
    if (negative) x = -x;
    if (x > cdfTableLimit) return negative ? cdf(-x) : cdf(x);
    double p = (cdfTableSize - 1) * x / cdfTableLimit;
    int below = (int) p;
    if (p == below) return negative ? -cdfTable[below] : cdfTable[below];
    int above = below + 1;
    double ret = cdfTable[below] +
            (cdfTable[above] - cdfTable[below])*(p - below);
    return negative ? 1 - ret : ret;
}

void print(const double* arr, int len) {
    for (int i = 0; i < len; ++i) {
        printf("%e; ", arr[i]);
    }
    puts("");
}

void print(const int* arr, int len) {
    for (int i = 0; i < len; ++i) {
        printf("%d; ", arr[i]);
    }
    puts("");
}

void fillBuckets(int N, int bucketCount,
        double* data, int* partitions,
        double* buckets, int* offsets) {
    for (int i = 0; i < N; ++i) {
        ++offsets[partitions[i]];
    }

    int offset = 0;
    for (int i = 0; i < bucketCount; ++i) {
        int t = offsets[i];
        offsets[i] = offset;
        offset += t;
    }
    offsets[bucketCount] = N;

    int next[bucketCount];
    memset(next, 0, sizeof(next));
    for (int i = 0; i < N; ++i) {
        int p = partitions[i];
        int j = offsets[p] + next[p];
        ++next[p];
        buckets[j] = data[i];
    }
}

class Sorter {
public:
    Sorter(double* data, int* offsets) {
        this->data = data;
        this->offsets = offsets;
    }

    static void radixSort(double* arr, int len) {
        ull* encoded = (ull*)arr;
        for (int i = 0; i < len; ++i) {
            ull n = encoded[i];
            if (n & signBit) {
                n ^= allBits;
            } else {
                n ^= signBit;
            }
            encoded[i] = n;
        }

        const int step = 11;
        const ull mask = (1ull << step) - 1;
        int offsets[8][1ull << step];
        memset(offsets, 0, sizeof(offsets));

        for (int i = 0; i < len; ++i) {
            for (int b = 0, j = 0; b < 64; b += step, ++j) {
                int p = (encoded[i] >> b) & mask;
                ++offsets[j][p];
            }
        }

        int sum[8] = {0};
        for (int i = 0; i <= mask; i++) {
            for (int b = 0, j = 0; b < 64; b += step, ++j) {
                int t = sum[j] + offsets[j][i];
                offsets[j][i] = sum[j];
                sum[j] = t;
            }
        }

        ull* copy = new ull[len];
        ull* current = encoded;
        for (int b = 0, j = 0; b < 64; b += step, ++j) {
            for (int i = 0; i < len; ++i) {
                int p = (current[i] >> b) & mask;
                copy[offsets[j][p]] = current[i];
                ++offsets[j][p];
            }

            ull* t = copy;
            copy = current;
            current = t;
        }

        if (current != encoded) {
            for (int i = 0; i < len; ++i) {
                encoded[i] = current[i];
            }
        }

        for (int i = 0; i < len; ++i) {
            ull n = encoded[i];
            if (n & signBit) {
                n ^= signBit;
            } else {
                n ^= allBits;
            }
            encoded[i] = n;
        }
    }

    void operator() (tbb::blocked_range<int>& range) const {
        for (int i = range.begin(); i < range.end(); ++i) {
            double* begin = &data[offsets[i]];
            double* end = &data[offsets[i+1]];
            //std::sort(begin, end);
            radixSort(begin, end-begin);
        }
    }

private:
    double* data;
    int* offsets;
    static const ull signBit = 1ull << 63;
    static const ull allBits = ~0ull;
};

void sortBuckets(int bucketCount, double* data, int* offsets) {
    Sorter sorter(data, offsets);
    tbb::blocked_range<int> range(0, bucketCount);
    tbb::parallel_for(range, sorter);
    //sorter(range);
}

class Partitioner {
public:
    Partitioner(int bucketCount, double* data, int* partitions) {
        this->data = data;
        this->partitions = partitions;
        this->bucketCount = bucketCount;
    }

    void operator() (tbb::blocked_range<int>& range) const {
        for (int i = range.begin(); i < range.end(); ++i) {
            double d = data[i];
            int p = (int) (cdfApprox(d) * bucketCount);
            partitions[i] = p;
        }
    }

private:
    double* data;
    int* partitions;
    int bucketCount;
};

const int bucketCount = 512;
int offsets[bucketCount + 1];

int main(int argc, char** argv) {
    if (argc != 2) {
        printf("Usage: %s N\n N = the size of the input\n", argv[0]);
        return 1;
    }

    puts("initializing...");
    int N = atoi(argv[1]);
    double* data = new double[N];
    double* buckets = new double[N];
    memset(offsets, 0, sizeof(offsets));
    int* partitions = new int[N];

    puts("loading data...");
    FILE* fp = fopen("gaussian.dat", "rb");
    if (fp == 0 || fread(data, sizeof(*data), N, fp) != N) {
        puts("Error reading data");
        return 1;
    }
    //print(data, N);

    puts("assigning partitions...");
    tbb::parallel_for(tbb::blocked_range<int>(0, N),
            Partitioner(bucketCount, data, partitions));

    puts("filling buckets...");
    fillBuckets(N, bucketCount, data, partitions, buckets, offsets);
    data = buckets;

    puts("sorting buckets...");
    sortBuckets(bucketCount, data, offsets);

    puts("done.");

    /*
    for (int i = 0; i < N-1; ++i) {
        if (data[i] > data[i+1]) {
            printf("error at %d: %e > %e\n", i, data[i], data[i+1]);
        }
    }
    */

    //print(data, N);

    return 0;
}

Чтобы скомпилировать и запустить его, используйте эту команду:

g++ -O3 -ltbb -o gsort gsort.cpp && time ./gsort 50000000

РЕДАКТИРОВАТЬ: Все сегменты теперь помещены в один и тот же массив, чтобы избавить от необходимости копировать сегменты обратно в массив. Также был уменьшен размер таблицы с предварительно вычисленными значениями, поскольку значения достаточно точны. Тем не менее, если я изменю количество сегментов выше 256, программа будет работать дольше, чем с таким количеством сегментов.

РЕДАКТИРОВАТЬ: один и тот же алгоритм, другой язык программирования. Я использовал C ++ вместо Java, и время работы на моей машине сократилось с ~ 3.2 с до 2.35 с. Оптимальное количество сегментов по-прежнему около 256 (опять же на моем компьютере).

Кстати, tbb действительно круто.

РЕДАКТИРОВАТЬ: Я был вдохновлен отличным решением Александру и заменил std :: sort на последнем этапе модифицированной версией его radix sort. Я использовал другой метод для работы с положительными / отрицательными числами, даже если ему нужно больше проходов через массив. Я также решил точно отсортировать массив и удалить сортировку вставки. Позже я потрачу некоторое время на тестирование того, как эти изменения влияют на производительность и, возможно, возвращают их. Однако с помощью радикальной сортировки время сократилось с ~ 2,35 до ~ 1,63 с.

Без умного, просто чтобы обеспечить намного более быстрый наивный сортировщик, вот один в C, который должен быть в значительной степени эквивалентен вашему Python:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
int cmp(const void* av, const void* bv) {
    double a = *(const double*)av;
    double b = *(const double*)bv;
    return a < b ? -1 : a > b ? 1 : 0;
}
int main(int argc, char** argv) {
    if (argc <= 1)
        return puts("No argument!");
    unsigned count = atoi(argv[1]);

    double *a = malloc(count * sizeof *a);

    FILE *f = fopen("gaussian.dat", "rb");
    if (fread(a, sizeof *a, count, f) != count)
        return puts("fread failed!");
    fclose(f);

    puts("sorting...");
    double *b = malloc(count * sizeof *b);
    memcpy(b, a, count * sizeof *b);
    qsort(b, count, sizeof *b, cmp);
    return 0;
}

Скомпилировано с gcc -O3, на моей машине это занимает более минуты меньше, чем Python: около 11 с по сравнению с 87 с.

Я разделил на сегменты на основе стандартного отклонения, которое лучше всего разбить на 4-е. Изменить: переписать на раздел на основе значения x в http://en.wikipedia.org/wiki/Error_function#Table_of_values

http://www.wolframalpha.com/input/?i=percentages+by++normal+distribution

Я пытался использовать меньшие сегменты, но, как мне показалось, это оказало небольшое влияние на 2 * сверх количества доступных ядер. Без каких-либо параллельных коллекций это займет 37 секунд на моем боксе и 24 с параллельными коллекциями. При разделении с помощью распределения вы не можете просто использовать массив, так что есть некоторые дополнительные издержки. Я не ясно, когда значение будет упаковано / распаковано в Scala.

Я использую Scala 2.9, для параллельной коллекции. Вы можете просто скачать дистрибутив tar.gz.

Для компиляции: scalac SortFile.scala (я просто скопировал его прямо в папку scala / bin.

Для запуска: JAVA_OPTS = "- Xmx4096M" ./scala SortFile (я запустил его с 2 гигабайтами оперативной памяти и получил примерно в то же время)

Изменить: Удален allocateDirect, медленнее, чем просто распределить. Убрано заполнение начального размера для буферов массива. На самом деле сделал это прочитать все 50000000 значений. Переписано, чтобы избежать проблем с автобоксом (все еще медленнее, чем наивный c)

import java.io.FileInputStream;
import java.nio.ByteBuffer
import java.nio.ByteOrder
import scala.collection.mutable.ArrayBuilder


object SortFile {

//used partition numbers from Damascus' solution
val partList = List(0, 0.15731, 0.31864, 0.48878, 0.67449, 0.88715, 1.1503, 1.5341)

val listSize = partList.size * 2;
val posZero = partList.size;
val neg = partList.map( _ * -1).reverse.zipWithIndex
val pos = partList.map( _ * 1).zipWithIndex.reverse

def partition(dbl:Double): Int = { 

//for each partition, i am running through the vals in order
//could make this a binary search to be more performant... but our list size is 4 (per side)

  if(dbl < 0) { return neg.find( dbl < _._1).get._2  }
  if(dbl > 0) { return posZero  + pos.find( dbl > _._1).get._2  }
      return posZero; 

}

  def main(args: Array[String])
    { 

    var l = 0
    val dbls = new Array[Double](50000000)
    val partList = new Array[Int](50000000)
    val pa = Array.fill(listSize){Array.newBuilder[Double]}
    val channel = new FileInputStream("../../gaussian.dat").getChannel()
    val bb = ByteBuffer.allocate(50000000 * 8)
    bb.order(ByteOrder.LITTLE_ENDIAN)
    channel.read(bb)
    bb.rewind
    println("Loaded" + System.currentTimeMillis())
    var dbl = 0.0
    while(bb.hasRemaining)
    { 
      dbl = bb.getDouble
      dbls.update(l,dbl) 

      l+=1
    }
    println("Beyond first load" + System.currentTimeMillis());

    for( i <- (0 to 49999999).par) { partList.update(i, partition(dbls(i)))}

    println("Partition computed" + System.currentTimeMillis() )
    for(i <- (0 to 49999999)) { pa(partList(i)) += dbls(i) }
    println("Partition completed " + System.currentTimeMillis())
    val toSort = for( i <- pa) yield i.result()
    println("Arrays Built" + System.currentTimeMillis());
    toSort.par.foreach{i:Array[Double] =>scala.util.Sorting.quickSort(i)};

    println("Read\t" + System.currentTimeMillis());

  }
}

Просто поместите это в файл cs и скомпилируйте его с помощью csc в теории: (Требуется моно)

using System;
using System.IO;
using System.Threading;

namespace Sort
{
    class Program
    {
        const int count = 50000000;
        static double[][] doubles;
        static WaitHandle[] waiting = new WaitHandle[4];
        static AutoResetEvent[] events = new AutoResetEvent[4];

        static double[] Merge(double[] left, double[] right)
        {
            double[] result = new double[left.Length + right.Length];
            int l = 0, r = 0, spot = 0;
            while (l < left.Length && r < right.Length)
            {
                if (right[r] < left[l])
                    result[spot++] = right[r++];
                else
                    result[spot++] = left[l++];
            }
            while (l < left.Length) result[spot++] = left[l++];
            while (r < right.Length) result[spot++] = right[r++];
            return result;
        }

        static void ThreadStart(object data)
        {
            int index = (int)data;
            Array.Sort(doubles[index]);
            events[index].Set();
        }

        static void Main(string[] args)
        {
            System.Diagnostics.Stopwatch watch = new System.Diagnostics.Stopwatch();
            watch.Start();
            byte[] bytes = File.ReadAllBytes(@"..\..\..\SortGuassian\Data.dat");
            doubles = new double[][] { new double[count / 4], new double[count / 4], new double[count / 4], new double[count / 4] };
            for (int i = 0; i < 4; i++)
            {
                for (int j = 0; j < count / 4; j++)
                {
                    doubles[i][j] = BitConverter.ToDouble(bytes, i * count/4 + j * 8);
                }
            }
            Thread[] threads = new Thread[4];
            for (int i = 0; i < 4; i++)
            {
                threads[i] = new Thread(ThreadStart);
                waiting[i] = events[i] = new AutoResetEvent(false);
                threads[i].Start(i);
            }
            WaitHandle.WaitAll(waiting);
            double[] left = Merge(doubles[0], doubles[1]);
            double[] right = Merge(doubles[2], doubles[3]);
            double[] result = Merge(left, right);
            watch.Stop();
            Console.WriteLine(watch.Elapsed.ToString());
            Console.ReadKey();
        }
    }
}

Поскольку вы знаете, что такое распределение, вы можете использовать прямую индексацию O (N). (Если вам интересно, что это такое, предположим, что у вас есть колода из 52 карт, и вы хотите ее отсортировать. Просто наберите 52 ячейки и бросьте каждую карточку в отдельную ячейку.)

У вас есть 5e7 дублей. Выделите массив результатов R из 5e7 double. Возьми каждый номер xи получи i = phi(x) * 5e7. В основном делать R[i] = x. Иметь способ обработки коллизий, таких как перемещение числа, с которым он может столкнуться (как в простом хэш-кодировании). В качестве альтернативы, вы можете сделать R в несколько раз больше, заполнить уникальным пустым значением. В конце вы просто подметаете элементы R.

phiэто просто гауссовская функция распределения. Он преобразует гауссово распределенное число между +/- бесконечностью в равномерное распределенное число между 0 и 1. Простой способ его вычисления - поиск в таблице и интерполяция.

Вот еще одно последовательное решение:

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <ctime>

typedef unsigned long long ull;

int size;
double *dbuf, *copy;
int cnt[8][1 << 16];

void sort()
{
  const int step = 10;
  const int start = 24;
  ull mask = (1ULL << step) - 1;

  ull *ibuf = (ull *) dbuf;
  for (int i = 0; i < size; i++) {
    for (int w = start, v = 0; w < 64; w += step, v++) {
      int p = (~ibuf[i] >> w) & mask;
      cnt[v][p]++;
    }
  }

  int sum[8] = { 0 };
  for (int i = 0; i <= mask; i++) {
    for (int w = start, v = 0; w < 64; w += step, v++) {
      int tmp = sum[v] + cnt[v][i];
      cnt[v][i] = sum[v];
      sum[v] = tmp;
    }
  }

  for (int w = start, v = 0; w < 64; w += step, v++) {
    ull *ibuf = (ull *) dbuf;
    for (int i = 0; i < size; i++) {
      int p = (~ibuf[i] >> w) & mask;
      copy[cnt[v][p]++] = dbuf[i];
    }

    double *tmp = copy;
    copy = dbuf;
    dbuf = tmp;
  }

  for (int p = 0; p < size; p++)
    if (dbuf[p] >= 0.) {
      std::reverse(dbuf + p, dbuf + size);
      break;
    }

  // Insertion sort
  for (int i = 1; i < size; i++) {
    double value = dbuf[i];
    if (value < dbuf[i - 1]) {
      dbuf[i] = dbuf[i - 1];
      int p = i - 1;
      for (; p > 0 && value < dbuf[p - 1]; p--)
        dbuf[p] = dbuf[p - 1];
      dbuf[p] = value;
    }
  }
}

int main(int argc, char **argv) {
  size = atoi(argv[1]);
  dbuf = new double[size];
  copy = new double[size];

  FILE *f = fopen("gaussian.dat", "r");
  fread(dbuf, size, sizeof(double), f);
  fclose(f);

  clock_t c0 = clock();
  sort();
  printf("Finished after %.3f\n", (double) ((clock() - c0)) / CLOCKS_PER_SEC);
  return 0;
}

Я сомневаюсь, что это превосходит многопоточное решение, но время на моем ноутбуке i7 таково (stdsort - это решение C ++, представленное в другом ответе):

$ g++ -O3 mysort.cpp -o mysort && ./mysort 50000000
Finished after 2.10
$ g++ -O3 stdsort.cpp -o stdsort && ./stdsort
Finished after 7.12

Обратите внимание, что это решение имеет линейную сложность по времени (потому что оно использует специальное представление двойников).

РЕДАКТИРОВАТЬ : Исправлен порядок увеличения элементов.

РЕДАКТИРОВАТЬ : Улучшена скорость почти на полсекунды.

РЕДАКТИРОВАТЬ : Улучшена скорость еще на 0,7 секунды. Сделал алгоритм более дружественным к кешу.

РЕДАКТИРОВАТЬ : Улучшена скорость еще на 1 секунду. Так как там только 50.000.000 элементов, я могу частично отсортировать мантиссу и использовать вставку сортировки (которая дружественна кешу) для исправления неуместных элементов. Эта идея удаляет около двух итераций из последнего цикла сортировки по основанию.

РЕДАКТИРОВАТЬ : 0,16 меньше секунд. Первый std :: reverse можно удалить, если порядок сортировки обратный.

Взять решение Кристиана Аммера и распараллелить его с многопоточными блоками Intel

#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <ctime>
#include <tbb/parallel_sort.h>

int main(void)
{
    std::ifstream ifs("gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<double> values;
    values.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
    values.push_back(d);
    clock_t c0 = clock();
    tbb::parallel_sort(values.begin(), values.end());
    std::cout << "Finished after "
              << static_cast<double>((clock() - c0)) / CLOCKS_PER_SEC
              << std::endl;
}

Если у вас есть доступ к библиотеке Intel Performance Primitives (IPP), вы можете использовать ее основную сортировку. Просто замени

#include <tbb/parallel_sort.h>

с участием

#include "ipps.h"

а также

tbb::parallel_sort(values.begin(), values.end());

с участием

std::vector<double> copy(values.size());
ippsSortRadixAscend_64f_I(&values[0], &copy[0], values.size());

На моем двухъядерном ноутбуке время

C               16.4 s
C#              20 s
C++ std::sort   7.2 s
C++ tbb         5 s
C++ ipp         4.5 s
python          too long

Как насчет реализации параллельной быстрой сортировки, которая выбирает свои основные значения на основе статистики распределения, тем самым обеспечивая разделы равного размера? Первый круг будет в среднем (в данном случае ноль), следующая пара будет в 25-м и 75-м процентилях (+/- -0,67449 стандартных отклонений) и т. Д., При этом каждый раздел делит оставшийся набор данных пополам вдвое больше или менее идеально.

Очень некрасиво (зачем использовать массивы, когда я могу использовать переменные, заканчивающиеся цифрами), но быстрый код (моя первая попытка std :: threads), все время (реальное время) в моей системе 1,8 с (по сравнению с std :: sort () 4,8 с), скомпилировать с помощью g ++ -std = c ++ 0x -O3 -march = native -pthread Просто передать данные через стандартный ввод (работает только для 50M).

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <thread>
using namespace std;
const size_t size=50000000;

void pivot(double* start,double * end, double middle,size_t& koniec){
    double * beg=start;
    end--;
    while (start!=end){
        if (*start>middle) swap (*start,*end--);
        else start++;
    }
    if (*end<middle) start+=1;
    koniec= start-beg;
}
void s(double * a, double* b){
    sort(a,b);
}
int main(){
    double *data=new double[size];
    FILE *f = fopen("gaussian.dat", "rb");
    fread(data,8,size,f);
    size_t end1,end2,end3,temp;
    pivot(data, data+size,0,end2);
    pivot(data, data+end2,-0.6745,end1);
    pivot(data+end2,data+size,0.6745,end3);
    end3+=end2;
    thread ts1(s,data,data+end1);
    thread ts2(s,data+end1,data+end2);
    thread ts3(s,data+end2,data+end3);
    thread ts4(s,data+end3,data+size);
    ts1.join(),ts2.join(),ts3.join(),ts4.join();
    //for (int i=0; i<size-1; i++){
    //  if (data[i]>data[i+1]) cerr<<"BLAD\n";
    //}
    fclose(f);
    //fwrite(data,8,size,stdout);
}

// Изменить, чтобы прочитать файл gaussian.dat.

Решение на C ++, использующее std::sort(в конечном итоге быстрее, чем qsort, относительно производительности qsort vs std :: sort )

#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <ctime>

int main(void)
{
    std::ifstream ifs("C:\\Temp\\gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<double> values;
    values.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
        values.push_back(d);
    clock_t c0 = clock();
    std::sort(values.begin(), values.end());
    std::cout << "Finished after "
              << static_cast<double>((clock() - c0)) / CLOCKS_PER_SEC
              << std::endl;
}

Я не могу с уверенностью сказать, сколько времени это займет, потому что у меня есть только 1 ГБ на моей машине, и с данным кодом Python я мог бы сделать gaussian.datфайл только с двойными 25 млн. (Без ошибки памяти). Но мне очень интересно, как долго работает алгоритм std :: sort.

Вот смесь радикальной сортировки Александру и умного поворота Жарека. Скомпилируйте это с

g++ -std=c++0x -pthread -O3 -march=native sorter_gaussian_radix.cxx -o sorter_gaussian_radix

Вы можете изменить размер радиуса, определив STEP (например, добавить -DSTEP = 11). Я нашел лучший для моего ноутбука 8 (по умолчанию).

По умолчанию он разбивает проблему на 4 части и запускает ее в нескольких потоках. Вы можете изменить это, передав параметр глубины в командную строку. Так что, если у вас есть два ядра, запустите его как

sorter_gaussian_radix 50000000 1

и если у вас есть 16 ядер

sorter_gaussian_radix 50000000 4

Максимальная глубина сейчас составляет 6 (64 потока). Если вы ставите слишком много уровней, вы просто замедляете код.

Одна вещь, которую я также попробовал, была сортировка radix из библиотеки Intel Performance Primitives (IPP). Реализация Александру явно затрагивает IPP, причем IPP работает примерно на 30% медленнее. Это изменение также включено сюда (закомментировано).

#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <ctime>
#include <iostream>
#include <thread>
#include <vector>
#include <boost/cstdint.hpp>
// #include "ipps.h"

#ifndef STEP
#define STEP 8
#endif

const int step = STEP;
const int start_step=24;
const int num_steps=(64-start_step+step-1)/step;
int size;
double *dbuf, *copy;

clock_t c1, c2, c3, c4, c5;

const double distrib[]={-2.15387,
                        -1.86273,
                        -1.67594,
                        -1.53412,
                        -1.4178,
                        -1.31801,
                        -1.22986,
                        -1.15035,
                        -1.07752,
                        -1.00999,
                        -0.946782,
                        -0.887147,
                        -0.830511,
                        -0.776422,
                        -0.724514,
                        -0.67449,
                        -0.626099,
                        -0.579132,
                        -0.53341,
                        -0.488776,
                        -0.445096,
                        -0.40225,
                        -0.36013,
                        -0.318639,
                        -0.27769,
                        -0.237202,
                        -0.197099,
                        -0.157311,
                        -0.11777,
                        -0.0784124,
                        -0.0391761,
                        0,
                        0.0391761,
                        0.0784124,
                        0.11777,
                        0.157311,
                        0.197099,
                        0.237202,
                        0.27769,
                        0.318639,
                        0.36013,
                        0.40225,
                        0.445097,
                        0.488776,
                        0.53341,
                        0.579132,
                        0.626099,
                        0.67449,
                        0.724514,
                        0.776422,
                        0.830511,
                        0.887147,
                        0.946782,
                        1.00999,
                        1.07752,
                        1.15035,
                        1.22986,
                        1.31801,
                        1.4178,
                        1.53412,
                        1.67594,
                        1.86273,
                        2.15387};


class Distrib
{
  const int value;
public:
  Distrib(const double &v): value(v) {}

  bool operator()(double a)
  {
    return a<value;
  }
};


void recursive_sort(const int start, const int end,
                    const int index, const int offset,
                    const int depth, const int max_depth)
{
  if(depth<max_depth)
    {
      Distrib dist(distrib[index]);
      const int middle=std::partition(dbuf+start,dbuf+end,dist) - dbuf;

      // const int middle=
      //   std::partition(dbuf+start,dbuf+end,[&](double a)
      //                  {return a<distrib[index];})
      //   - dbuf;

      std::thread lower(recursive_sort,start,middle,index-offset,offset/2,
                        depth+1,max_depth);
      std::thread upper(recursive_sort,middle,end,index+offset,offset/2,
                        depth+1,max_depth);
      lower.join(), upper.join();
    }
  else
    {
  // ippsSortRadixAscend_64f_I(dbuf+start,copy+start,end-start);

      c1=clock();

      double *dbuf_local(dbuf), *copy_local(copy);
      boost::uint64_t mask = (1 << step) - 1;
      int cnt[num_steps][mask+1];

      boost::uint64_t *ibuf = reinterpret_cast<boost::uint64_t *> (dbuf_local);

      for(int i=0;i<num_steps;++i)
        for(uint j=0;j<mask+1;++j)
          cnt[i][j]=0;

      for (int i = start; i < end; i++)
        {
          for (int w = start_step, v = 0; w < 64; w += step, v++)
            {
              int p = (~ibuf[i] >> w) & mask;
              (cnt[v][p])++;
            }
        }

      c2=clock();

      std::vector<int> sum(num_steps,0);
      for (uint i = 0; i <= mask; i++)
        {
          for (int w = start_step, v = 0; w < 64; w += step, v++)
            {
              int tmp = sum[v] + cnt[v][i];
              cnt[v][i] = sum[v];
              sum[v] = tmp;
            }
        }

      c3=clock();

      for (int w = start_step, v = 0; w < 64; w += step, v++)
        {
          ibuf = reinterpret_cast<boost::uint64_t *>(dbuf_local);

          for (int i = start; i < end; i++)
            {
              int p = (~ibuf[i] >> w) & mask;
              copy_local[start+((cnt[v][p])++)] = dbuf_local[i];
            }
          std::swap(copy_local,dbuf_local);
        }

      // Do the last set of reversals
      for (int p = start; p < end; p++)
        if (dbuf_local[p] >= 0.)
          {
            std::reverse(dbuf_local+p, dbuf_local + end);
            break;
          }

      c4=clock();

      // Insertion sort
      for (int i = start+1; i < end; i++) {
        double value = dbuf_local[i];
        if (value < dbuf_local[i - 1]) {
          dbuf_local[i] = dbuf_local[i - 1];
          int p = i - 1;
          for (; p > 0 && value < dbuf_local[p - 1]; p--)
            dbuf_local[p] = dbuf_local[p - 1];
          dbuf_local[p] = value;
        }
      }
      c5=clock();

    }
}


int main(int argc, char **argv) {
  size = atoi(argv[1]);
  copy = new double[size];

  dbuf = new double[size];
  FILE *f = fopen("gaussian.dat", "r");
  fread(dbuf, size, sizeof(double), f);
  fclose(f);

  clock_t c0 = clock();

  const int max_depth= (argc > 2) ? atoi(argv[2]) : 2;

  // ippsSortRadixAscend_64f_I(dbuf,copy,size);

  recursive_sort(0,size,31,16,0,max_depth);

  if(num_steps%2==1)
    std::swap(dbuf,copy);

  // for (int i=0; i<size-1; i++){
  //   if (dbuf[i]>dbuf[i+1])
  //     std::cout << "BAD "
  //               << i << " "
  //               << dbuf[i] << " "
  //               << dbuf[i+1] << " "
  //               << "\n";
  // }

  std::cout << "Finished after "
            << (double) (c1 - c0) / CLOCKS_PER_SEC << " "
            << (double) (c2 - c1) / CLOCKS_PER_SEC << " "
            << (double) (c3 - c2) / CLOCKS_PER_SEC << " "
            << (double) (c4 - c3) / CLOCKS_PER_SEC << " "
            << (double) (c5 - c4) / CLOCKS_PER_SEC << " "
            << "\n";

  // delete [] dbuf;
  // delete [] copy;
  return 0;
}

РЕДАКТИРОВАТЬ : я реализовал улучшения кэша Александру, и это сэкономило около 30% времени на моей машине.

РЕДАКТИРОВАТЬ : Это реализует рекурсивную сортировку, поэтому он должен хорошо работать на 16-ядерном компьютере Александру. Он также использует последнее улучшение Александру и удаляет одно из обратного. Для меня это дало улучшение на 20%.

РЕДАКТИРОВАТЬ : Исправлена ​​ошибка со знаком, которая приводила к неэффективности при наличии более 2 ядер.

РЕДАКТИРОВАТЬ : Удалена лямбда, поэтому он будет компилироваться с более старыми версиями gcc. Он включает закомментированное изменение кода IPP. Я также исправил документацию для работы на 16 ядрах. Насколько я могу судить, это самая быстрая реализация.

РЕДАКТИРОВАТЬ : Исправлена ​​ошибка, когда STEP не 8. Увеличено максимальное количество потоков до 64. Добавлена ​​некоторая информация о времени.

Я думаю, это действительно зависит от того, что вы хотите сделать. Если вы хотите отсортировать группу гауссов, то это вам не поможет. Но если вы хотите кучу отсортированных гауссиан, это будет. Даже если это немного не решит проблему, я думаю, что будет интересно сравнить ее с реальными процедурами сортировки.

Если хочешь, чтобы что-то было быстрым, делай меньше.

Вместо генерации группы случайных выборок из нормального распределения и последующей сортировки можно создать группу выборок из нормального распределения в отсортированном порядке.

Вы можете использовать решение здесь для создания п однородных случайных чисел в отсортированном порядке. Затем вы можете использовать обратный cdf (scipy.stats.norm.ppf) нормального распределения, чтобы превратить единообразные случайные числа в числа из нормального распределения посредством выборки обратного преобразования .

import scipy.stats
import random

# slightly modified from linked stackoverflow post
def n_random_numbers_increasing(n):
  """Like sorted(random() for i in range(n))),                                
  but faster because we avoid sorting."""
  v = 1.0
  while n:
    v *= random.random() ** (1.0 / n)
    yield 1 - v
    n -= 1

def n_normal_samples_increasing(n):
  return map(scipy.stats.norm.ppf, n_random_numbers_increasing(n))

Если вы хотите сделать свои руки более грязными, я думаю, вы могли бы ускорить многие обратные вычисления cdf, используя какой-то итерационный метод и используя предыдущий результат в качестве первоначального предположения. Поскольку предположения будут очень близки, вероятно, одна итерация даст вам большую точность.

Попробуйте изменить решение Guvante с помощью функции Main (), она начнет сортировку, как только будет выполнено чтение 1/4 ввода-вывода, и в моем тесте она будет быстрее:

    static void Main(string[] args)
    {
        FileStream filestream = new FileStream(@"..\..\..\gaussian.dat", FileMode.Open, FileAccess.Read);
        doubles = new double[][] { new double[count / 4], new double[count / 4], new double[count / 4], new double[count / 4] };
        Thread[] threads = new Thread[4];

        for (int i = 0; i < 4; i++)
        {
            byte[] bytes = new byte[count * 4];
            filestream.Read(bytes, 0, count * 4);

            for (int j = 0; j < count / 4; j++)
            {
                doubles[i][j] = BitConverter.ToDouble(bytes, i * count/4 + j * 8);
            }

            threads[i] = new Thread(ThreadStart);
            waiting[i] = events[i] = new AutoResetEvent(false);
            threads[i].Start(i);    
        }

        WaitHandle.WaitAll(waiting);
        double[] left = Merge(doubles[0], doubles[1]);
        double[] right = Merge(doubles[2], doubles[3]);
        double[] result = Merge(left, right);
        Console.ReadKey();
    }
}

Поскольку вы знаете распределение, моя идея состоит в том, чтобы создать k блоков, каждый из которых имеет одинаковое ожидаемое количество элементов (поскольку вы знаете распределение, вы можете вычислить это). Затем за O (n) время сметает массив и помещает элементы в их сегменты.

Затем одновременно сортируйте ведра. Предположим, у вас есть k блоков и n элементов. Сортировка займет (n / k) lg (n / k) время. Теперь предположим, что у вас есть p процессоров, которые вы можете использовать. Так как ведра могут быть отсортированы независимо, у вас есть множитель ceil (k / p), чтобы иметь дело с. Это дает окончательное время выполнения n + ceil (k / p) * (n / k) lg (n / k), что должно быть намного быстрее, чем n lg n, если вы выбираете k хорошо.

Одна из низкоуровневых идей оптимизации состоит в том, чтобы поместить два двойных в регистр SSE, чтобы каждый поток работал с двумя элементами одновременно. Это может быть сложно сделать для некоторых алгоритмов.

Еще одна вещь, которую нужно сделать, это отсортировать массив в кеширующие части, а затем объединить результаты. Следует использовать два уровня: например, сначала 4 КБ для L1, затем 64 КБ для L2.

Это должно быть очень дружественным к кешу, так как сортировка сегментов не будет выходить за пределы кэша, а окончательное слияние будет последовательно обходить память.

В наши дни вычисления намного дешевле, чем доступ к памяти. Однако у нас есть большое количество элементов, поэтому трудно сказать, какой размер массива, когда глупая сортировка с поддержкой кэша выполняется медленнее, чем версия с низкой сложностью, не поддерживающая кэширование.

Но я не буду предоставлять реализацию вышеупомянутого, поскольку я сделал бы это в Windows (VC ++).

Вот реализация сортировки с линейным сканированием. Я думаю, что это быстрее, чем все текущие однопоточные реализации, за исключением сортировки по основанию. Он должен иметь ожидаемое линейное время выполнения, если я достаточно точно оцениваю cdf (я использую линейную интерполяцию значений, которые я нашел в Интернете) и не допустил ошибок, которые могли бы вызвать чрезмерное сканирование:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <ctime>

using std::fill;

const double q[] = {
  0.0,
  9.865E-10,
  2.8665150000000003E-7,
  3.167E-5,
  0.001349898,
  0.022750132,
  0.158655254,
  0.5,
  0.8413447460000001,
  0.9772498679999999,
  0.998650102,
  0.99996833,
  0.9999997133485,
  0.9999999990134999,
  1.0,
};
int main(int argc, char** argv) {
  if (argc <= 1)
    return puts("No argument!");
  unsigned count = atoi(argv[1]);
  unsigned count2 = 3 * count;

  bool *ba = new bool[count2 + 1000];
  fill(ba, ba + count2 + 1000, false);
  double *a = new double[count];
  double *c = new double[count2 + 1000];

  FILE *f = fopen("gaussian.dat", "rb");
  if (fread(a, 8, count, f) != count)
    return puts("fread failed!");
  fclose(f);

  int i;
  int j;
  bool s;
  int t;
  double z;
  double p;
  double d1;
  double d2;
  for (i = 0; i < count; i++) {
    s = a[i] < 0;
    t = a[i];
    if (s) t--;
    z = a[i] - t;
    t += 7;
    if (t < 0) {
      t = 0;
      z = 0;
    } else if (t >= 14) {
      t = 13;
      z = 1;
    }
    p = q[t] * (1 - z) + q[t + 1] * z;
    j = count2 * p;
    while (ba[j] && c[j] < a[i]) {
      j++;
    }
    if (!ba[j]) {
      ba[j] = true;
      c[j] = a[i];
    } else {
      d1 = c[j];
      c[j] = a[i];
      j++;
      while (ba[j]) {
        d2 = c[j];
        c[j] = d1;
        d1 = d2;
        j++;
      }
      c[j] = d1;
      ba[j] = true;
    }
  }
  i = 0;
  int max = count2 + 1000;
  for (j = 0; j < max; j++) {
    if (ba[j]) {
      a[i++] = c[j];
    }
  }
  // for (i = 0; i < count; i += 1) {
  //   printf("here %f\n", a[i]);
  // }
  return 0;
}

Я не знаю, почему я не могу отредактировать свой предыдущий пост, поэтому вот новая версия, на 0,2 секунды быстрее (но примерно на 1,5 с быстрее во время процессора (пользователь)). Это решение имеет 2 программы, которые предварительно рассчитывают квантили для нормального распределения для сортировки сегментов и сохраняют их в таблице, t [double * scale] = индекс сегмента, где scale - произвольное число, которое делает возможным удвоение приведения. Тогда основная программа может использовать эти данные, чтобы поместить двойники в правильное ведро. У этого есть один недостаток: если данные не гауссовы, они не будут работать правильно (а также почти нулевой шанс работать некорректно для нормального распределения), но модификация для особого случая проста и быстра (только количество проверок сегментов и падение до стандартного значения) ::Сортировать()).

Компиляция: g ++ => http://pastebin.com/WG7pZEzH вспомогательная программа

g ++ -std = c ++ 0x -O3 -march = native -pthread => http://pastebin.com/T3yzViZP основная программа сортировки

Вот еще одно последовательное решение. Этот использует тот факт, что элементы распределены нормально, и я думаю, что идея в целом применима для сортировки, близкой к линейному времени.

Алгоритм такой:

• Приблизительный CDF (см. phi()Функцию в реализации)
• Для всех элементов вычисляем примерную позицию в отсортированном массиве: size * phi(x)
• Поместите элементы в новый массив близко к их конечной позиции
  • В моей реализации целевой массив имеет некоторые пробелы, поэтому мне не нужно сдвигать слишком много элементов при вставке.
• Используйте inserttsort для сортировки конечных элементов (inserttsort является линейным, если расстояние до конечной позиции меньше константы).

К сожалению, скрытая константа довольно велика, и это решение в два раза медленнее алгоритма сортировки по основанию.

Мой личный фаворит с использованием многопоточных строительных блоков Intel уже был опубликован, но вот грубое параллельное решение с использованием JDK 7 и его нового API-интерфейса fork / join:

import java.io.FileInputStream;
import java.nio.channels.FileChannel;
import java.util.Arrays;
import java.util.concurrent.*;
import static java.nio.channels.FileChannel.MapMode.READ_ONLY;
import static java.nio.ByteOrder.LITTLE_ENDIAN;


/**
 * 
 * Original Quicksort: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel
 *
 */
public class ForkJoinQuicksortTask extends RecursiveAction {

    public static void main(String[] args) throws Exception {

        double[] array = new double[Integer.valueOf(args[0])];

        FileChannel fileChannel = new FileInputStream("gaussian.dat").getChannel();
        fileChannel.map(READ_ONLY, 0, fileChannel.size()).order(LITTLE_ENDIAN).asDoubleBuffer().get(array);

        ForkJoinPool mainPool = new ForkJoinPool();

        System.out.println("Starting parallel computation");

        mainPool.invoke(new ForkJoinQuicksortTask(array));        
    }

    private static final long serialVersionUID = -642903763239072866L;
    private static final int SERIAL_THRESHOLD = 0x1000;

    private final double a[];
    private final int left, right;

    public ForkJoinQuicksortTask(double[] a) {this(a, 0, a.length - 1);}

    private ForkJoinQuicksortTask(double[] a, int left, int right) {
        this.a = a;
        this.left = left;
        this.right = right;
    }

    @Override
    protected void compute() {
        if (right - left < SERIAL_THRESHOLD) {
            Arrays.sort(a, left, right + 1);
        } else {
            int pivotIndex = partition(a, left, right);
            ForkJoinTask<Void> t1 = null;

            if (left < pivotIndex)
                t1 = new ForkJoinQuicksortTask(a, left, pivotIndex).fork();
            if (pivotIndex + 1 < right)
                new ForkJoinQuicksortTask(a, pivotIndex + 1, right).invoke();

            if (t1 != null)
                t1.join();
        }
    }

    public static int partition(double[] a, int left, int right) {
        // chose middle value of range for our pivot
        double pivotValue = a[left + (right - left) / 2];

        --left;
        ++right;

        while (true) {
            do
                ++left;
            while (a[left] < pivotValue);

            do
                --right;
            while (a[right] > pivotValue);

            if (left < right) {
                double tmp = a[left];
                a[left] = a[right];
                a[right] = tmp;
            } else {
                return right;
            }
        }
    }    
}

Важный отказ от ответственности : я использовал адаптацию быстрой сортировки для fork / join: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel

Для этого вам понадобится бета-версия JDK 7 (http://jdk7.java.net/download.html).

На моем 2.93 ГГц Quad Core i7 (OS X):

Python ссылка

time python sort.py 50000000
sorting...

real    1m13.885s
user    1m11.942s
sys     0m1.935s

Java JDK 7 форк / соединение

time java ForkJoinQuicksortTask 50000000
Starting parallel computation

real    0m2.404s
user    0m10.195s
sys     0m0.347s

Я также попытался поэкспериментировать с параллельным чтением и преобразованием байтов в удвоения, но я не увидел там никакой разницы.

Обновить:

Если кто-то хочет поэкспериментировать с параллельной загрузкой данных, версия с параллельной загрузкой приведена ниже. Теоретически это может сделать его еще немного быстрее, если ваше устройство ввода-вывода имеет достаточную параллельную емкость (обычно это делают SSD). Существуют также некоторые издержки при создании дублей из байтов, так что они потенциально могут работать быстрее параллельно. На моих системах (Ubuntu 10.10 / Nehalem Quad / Intel X25M SSD и OS X 10.6 / i7 Quad / Samsung SSD) я не увидел никакой реальной разницы.

import static java.nio.ByteOrder.LITTLE_ENDIAN;
import static java.nio.channels.FileChannel.MapMode.READ_ONLY;

import java.io.FileInputStream;
import java.nio.DoubleBuffer;
import java.nio.channels.FileChannel;
import java.util.Arrays;
import java.util.concurrent.ForkJoinPool;
import java.util.concurrent.ForkJoinTask;
import java.util.concurrent.RecursiveAction;


/**
 *
 * Original Quicksort: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel
 *
 */
public class ForkJoinQuicksortTask extends RecursiveAction {

   public static void main(String[] args) throws Exception {

       ForkJoinPool mainPool = new ForkJoinPool();

       double[] array = new double[Integer.valueOf(args[0])];
       FileChannel fileChannel = new FileInputStream("gaussian.dat").getChannel();
       DoubleBuffer buffer = fileChannel.map(READ_ONLY, 0, fileChannel.size()).order(LITTLE_ENDIAN).asDoubleBuffer();

       mainPool.invoke(new ReadAction(buffer, array, 0, array.length));
       mainPool.invoke(new ForkJoinQuicksortTask(array));
   }

   private static final long serialVersionUID = -642903763239072866L;
   private static final int SERIAL_THRESHOLD = 0x1000;

   private final double a[];
   private final int left, right;

   public ForkJoinQuicksortTask(double[] a) {this(a, 0, a.length - 1);}

   private ForkJoinQuicksortTask(double[] a, int left, int right) {
       this.a = a;
       this.left = left;
       this.right = right;
   }

   @Override
   protected void compute() {
       if (right - left < SERIAL_THRESHOLD) {
           Arrays.sort(a, left, right + 1);
       } else {
           int pivotIndex = partition(a, left, right);
           ForkJoinTask<Void> t1 = null;

           if (left < pivotIndex)
               t1 = new ForkJoinQuicksortTask(a, left, pivotIndex).fork();
           if (pivotIndex + 1 < right)
               new ForkJoinQuicksortTask(a, pivotIndex + 1, right).invoke();

           if (t1 != null)
               t1.join();
       }
   }

   public static int partition(double[] a, int left, int right) {
       // chose middle value of range for our pivot
       double pivotValue = a[left + (right - left) / 2];

       --left;
       ++right;

       while (true) {
           do
               ++left;
           while (a[left] < pivotValue);

           do
               --right;
           while (a[right] > pivotValue);

           if (left < right) {
               double tmp = a[left];
               a[left] = a[right];
               a[right] = tmp;
           } else {
               return right;
           }
       }
   }

}

class ReadAction extends RecursiveAction {

   private static final long serialVersionUID = -3498527500076085483L;

   private final DoubleBuffer buffer;
   private final double[] array;
   private final int low, high;

   public ReadAction(DoubleBuffer buffer, double[] array, int low, int high) {
       this.buffer = buffer;
       this.array = array;
       this.low = low;
       this.high = high;
   }

   @Override
   protected void compute() {
       if (high - low < 100000) {
           buffer.position(low);
           buffer.get(array, low, high-low);
       } else {
           int middle = (low + high) >>> 1;

           invokeAll(new ReadAction(buffer.slice(), array, low, middle),  new ReadAction(buffer.slice(), array, middle, high));
       }
   }
}

Update2:

Я выполнил код на одной из наших 12-ядерных машин с небольшой модификацией, чтобы установить фиксированное количество ядер. Это дало следующие результаты:

Cores  Time
1      7.568s
2      3.903s
3      3.325s
4      2.388s
5      2.227s
6      1.956s
7      1.856s
8      1.827s
9      1.682s
10     1.698s
11     1.620s
12     1.503s

В этой системе я также попробовал версию Python, которая занимала 1m2.994s, и версию Cjjarek C ++, которая занимала 1.925s (почему-то версия Zjarek C ++, кажется, работает относительно быстрее на компьютере static_rtti).

Я также попробовал, что случилось, если я удвоил размер файла до 100 000 000 удваивается:

Cores  Time
1      15.056s
2      8.116s
3      5.925s
4      4.802s
5      4.430s
6      3.733s
7      3.540s
8      3.228s
9      3.103s
10     2.827s
11     2.784s
12     2.689s

В этом случае версия Zjarek для C ++ заняла 3.968 с. Питон просто занял слишком много времени здесь.

150 000 000 пар:

Cores  Time
1      23.295s
2      12.391s
3      8.944s
4      6.990s
5      6.216s
6      6.211s
7      5.446s
8      5.155s
9      4.840s
10     4.435s
11     4.248s
12     4.174s

В этом случае версия Zjarek's C ++ была 6.044s. Я даже не пытался Python.

Версия C ++ очень согласуется с ее результатами, где Java немного качается. Сначала она становится немного более эффективной, когда проблема становится больше, но затем снова становится менее эффективной.

Версия с использованием традиционных pthreads. Код для слияния скопирован из ответа Гуванте. Компилировать с g++ -O3 -pthread.

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>
#include <algorithm>

static unsigned int nthreads = 4;
static unsigned int size = 50000000;

typedef struct {
  double *array;
  int size;
} array_t;


void 
merge(double *left, int leftsize,
      double *right, int rightsize,
      double *result)
{
  int l = 0, r = 0, insertat = 0;
  while (l < leftsize && r < rightsize) {
    if (left[l] < right[r])
      result[insertat++] = left[l++];
    else
      result[insertat++] = right[r++];
  }

  while (l < leftsize) result[insertat++] = left[l++];
  while (r < rightsize) result[insertat++] = right[r++];
}


void *
run_thread(void *input)
{
  array_t numbers = *(array_t *)input;
  std::sort(numbers.array, numbers.array+numbers.size); 
  pthread_exit(NULL);
}

int 
main(int argc, char **argv) 
{
  double *numbers = (double *) malloc(size * sizeof(double));

  FILE *f = fopen("gaussian.dat", "rb");
  if (fread(numbers, sizeof(double), size, f) != size)
    return printf("Reading gaussian.dat failed");
  fclose(f);

  array_t worksets[nthreads];
  int worksetsize = size / nthreads;
  for (int i = 0; i < nthreads; i++) {
    worksets[i].array=numbers+(i*worksetsize);
    worksets[i].size=worksetsize;
  }

  pthread_attr_t attributes;
  pthread_attr_init(&attributes);
  pthread_attr_setdetachstate(&attributes, PTHREAD_CREATE_JOINABLE);

  pthread_t threads[nthreads];
  for (int i = 0; i < nthreads; i++) {
    pthread_create(&threads[i], &attributes, &run_thread, &worksets[i]);
  }

  for (int i = 0; i < nthreads; i++) {
    pthread_join(threads[i], NULL);
  }

  double *tmp = (double *) malloc(size * sizeof(double));
  merge(numbers, worksetsize, numbers+worksetsize, worksetsize, tmp);
  merge(numbers+(worksetsize*2), worksetsize, numbers+(worksetsize*3), worksetsize, tmp+(size/2));
  merge(tmp, worksetsize*2, tmp+(size/2), worksetsize*2, numbers);

  /*
  printf("Verifying result..\n");
  for (int i = 0; i < size - 1; i++) {
    if (numbers[i] > numbers[i+1])
      printf("Result is not correct\n");
  }
  */

  pthread_attr_destroy(&attributes);
  return 0;
}  

На моем ноутбуке я получаю следующие результаты:

real    0m6.660s
user    0m9.449s
sys     0m1.160s

Вот последовательная реализация C99, которая пытается реально использовать известный дистрибутив. Он в основном выполняет один раунд сортировки сегментов с использованием информации о распределении, затем несколько циклов быстрой сортировки в каждом сегменте, предполагая равномерное распределение в пределах сегмента, и, наконец, измененную сортировку выбора для копирования данных обратно в исходный буфер. Быстрая сортировка запоминает точки разделения, поэтому сортировка выбора должна работать только с небольшими порциями. И несмотря на (потому что?) Всю эту сложность, это даже не очень быстро.

Чтобы сделать оценку Φ быстрой, значения выбираются в нескольких точках, а позже используется только линейная интерполяция. На самом деле не имеет значения, точно ли оценивается Φ, если приближение строго монотонно.

Размеры бункера выбираются таким образом, чтобы вероятность переполнения бункера была незначительной. Точнее, с текущими параметрами вероятность того, что набор данных из 50000000 элементов вызовет переполнение ячейки, равна 3.65e-09. (Это может быть вычислено с использованием функции выживания в распределении Пуассона ) .

Для компиляции, пожалуйста, используйте

gcc -std=c99 -msse3 -O3 -ffinite-math-only

Поскольку вычислений значительно больше, чем в других решениях, эти флаги компилятора необходимы для того, чтобы сделать его хотя бы достаточно быстрым. Без -msse3преобразований, doubleчтобы intстать очень медленно. Если ваша архитектура не поддерживает SSE3, эти преобразования также можно выполнить с помощью lrint()функции.

Код довольно уродливый - не уверен, что он соответствует требованию «разумно читаемого» ...

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <math.h>

#define N 50000000
#define BINSIZE 720
#define MAXBINSIZE 880
#define BINCOUNT (N / BINSIZE)
#define SPLITS 64
#define PHI_VALS 513

double phi_vals[PHI_VALS];

int bin_index(double x)
{
    double y = (x + 8.0) * ((PHI_VALS - 1) / 16.0);
    int interval = y;
    y -= interval;
    return (1.0 - y) * phi_vals[interval] + y * phi_vals[interval + 1];
}

double bin_value(int bin)
{
    int left = 0;
    int right = PHI_VALS - 1;
    do
    {
        int centre = (left + right) / 2;
        if (bin < phi_vals[centre])
            right = centre;
        else
            left = centre;
    } while (right - left > 1);
    double frac = (bin - phi_vals[left]) / (phi_vals[right] - phi_vals[left]);
    return (left + frac) * (16.0 / (PHI_VALS - 1)) - 8.0;
}

void gaussian_sort(double *restrict a)
{
    double *b = malloc(BINCOUNT * MAXBINSIZE * sizeof(double));
    double **pos = malloc(BINCOUNT * sizeof(double*));
    for (size_t i = 0; i < BINCOUNT; ++i)
        pos[i] = b + MAXBINSIZE * i;
    for (size_t i = 0; i < N; ++i)
        *pos[bin_index(a[i])]++ = a[i];
    double left_val, right_val = bin_value(0);
    for (size_t bin = 0, i = 0; bin < BINCOUNT; ++bin)
    {
        left_val = right_val;
        right_val = bin_value(bin + 1);
        double *splits[SPLITS + 1];
        splits[0] = b + bin * MAXBINSIZE;
        splits[SPLITS] = pos[bin];
        for (int step = SPLITS; step > 1; step >>= 1)
            for (int left_split = 0; left_split < SPLITS; left_split += step)
            {
                double *left = splits[left_split];
                double *right = splits[left_split + step] - 1;
                double frac = (double)(left_split + (step >> 1)) / SPLITS;
                double pivot = (1.0 - frac) * left_val + frac * right_val;
                while (1)
                {
                    while (*left < pivot && left <= right)
                        ++left;
                    while (*right >= pivot && left < right)
                        --right;
                    if (left >= right)
                        break;
                    double tmp = *left;
                    *left = *right;
                    *right = tmp;
                    ++left;
                    --right;
                }
                splits[left_split + (step >> 1)] = left;
            }
        for (int left_split = 0; left_split < SPLITS; ++left_split)
        {
            double *left = splits[left_split];
            double *right = splits[left_split + 1] - 1;
            while (left <= right)
            {
                double *min = left;
                for (double *tmp = left + 1; tmp <= right; ++tmp)
                    if (*tmp < *min)
                        min = tmp;
                a[i++] = *min;
                *min = *right--;
            }
        }
    }
    free(b);
    free(pos);
}

int main()
{
    double *a = malloc(N * sizeof(double));
    FILE *f = fopen("gaussian.dat", "rb");
    assert(fread(a, sizeof(double), N, f) == N);
    fclose(f);
    for (int i = 0; i < PHI_VALS; ++i)
    {
        double x = (i * (16.0 / PHI_VALS) - 8.0) / sqrt(2.0);
        phi_vals[i] =  (erf(x) + 1.0) * 0.5 * BINCOUNT;
    }
    gaussian_sort(a);
    free(a);
}

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <algorithm>

// maps [-inf,+inf] to (0,1)
double normcdf(double x) {
        return 0.5 * (1 + erf(x * M_SQRT1_2));
}

int calcbin(double x, int bins) {
        return (int)floor(normcdf(x) * bins);
}

int *docensus(int bins, int n, double *arr) {
        int *hist = calloc(bins, sizeof(int));
        int i;
        for(i = 0; i < n; i++) {
                hist[calcbin(arr[i], bins)]++;
        }
        return hist;
}

void partition(int bins, int *orig_counts, double *arr) {
        int *counts = malloc(bins * sizeof(int));
        memcpy(counts, orig_counts, bins*sizeof(int));
        int *starts = malloc(bins * sizeof(int));
        int b, i;
        starts[0] = 0;
        for(i = 1; i < bins; i++) {
                starts[i] = starts[i-1] + counts[i-1];
        }
        for(b = 0; b < bins; b++) {
                while (counts[b] > 0) {
                        double v = arr[starts[b]];
                        int correctbin;
                        do {
                                correctbin = calcbin(v, bins);
                                int swappos = starts[correctbin];
                                double tmp = arr[swappos];
                                arr[swappos] = v;
                                v = tmp;
                                starts[correctbin]++;
                                counts[correctbin]--;
                        } while (correctbin != b);
                }
        }
        free(counts);
        free(starts);
}


void sortbins(int bins, int *counts, double *arr) {
        int start = 0;
        int b;
        for(b = 0; b < bins; b++) {
                std::sort(arr + start, arr + start + counts[b]);
                start += counts[b];
        }
}


void checksorted(double *arr, int n) {
        int i;
        for(i = 1; i < n; i++) {
                if (arr[i-1] > arr[i]) {
                        printf("out of order at %d: %lf %lf\n", i, arr[i-1], arr[i]);
                        exit(1);
                }
        }
}


int main(int argc, char *argv[]) {
        if (argc == 1 || argv[1] == NULL) {
                printf("Expected data size as argument\n");
                exit(1);
        }
        int n = atoi(argv[1]);
        const int cachesize = 128 * 1024; // a guess
        int bins = (int) (1.1 * n * sizeof(double) / cachesize);
        if (argc > 2) {
                bins = atoi(argv[2]);
        }
        printf("Using %d bins\n", bins);
        FILE *f = fopen("gaussian.dat", "rb");
        if (f == NULL) {
                printf("Couldn't open gaussian.dat\n");
                exit(1);
        }
        double *arr = malloc(n * sizeof(double));
        fread(arr, sizeof(double), n, f);
        fclose(f);

        int *counts = docensus(bins, n, arr);
        partition(bins, counts, arr);
        sortbins(bins, counts, arr);
        checksorted(arr, n);

        return 0;
}

Это использует erf () для правильного размещения каждого элемента в ячейке, а затем сортирует каждую ячейку. Он сохраняет массив полностью на месте.

Первый проход: docensus () подсчитывает количество элементов в каждом бине.

Второй проход: partition () переставляет массив, помещая каждый элемент в соответствующую корзину

Третий этап: sortbins () выполняет qsort для каждого бина.

Это наивно и вызывает дорогую функцию erf () дважды для каждого значения. Первый и третий проходы потенциально распараллеливаются. Второй является очень последовательным и, вероятно, замедлен из-за его очень случайных образцов доступа к памяти. Также может быть целесообразно кэшировать номер каждой двойной ячейки, в зависимости от соотношения мощности процессора и скорости памяти.

Эта программа позволяет вам выбрать количество корзин для использования. Просто добавьте второй номер в командной строке. Я скомпилировал его с помощью gcc -O3, но моя машина настолько слабая, что не могу сказать вам хороших показателей производительности.

Редактировать: Пуф! Моя программа на C волшебным образом превратилась в программу на C ++ с использованием std :: sort!

Взгляните на реализацию сортировки по основанию от Michael Herf ( Radix Tricks ). На моей машине сортировка была в 5 раз быстрее по сравнению с std::sortалгоритмом в моем первом ответе. Имя функции сортировки является RadixSort11.

int main(void)
{
    std::ifstream ifs("C:\\Temp\\gaussian.dat", std::ios::binary | std::ios::in);
    std::vector<float> v;
    v.reserve(50000000);
    double d;
    while (ifs.read(reinterpret_cast<char*>(&d), sizeof(double)))
        v.push_back(static_cast<float>(d));
    std::vector<float> vres(v.size(), 0.0);
    clock_t c0 = clock();
    RadixSort11(&v[0], &vres[0], v.size());
    std::cout << "Finished after: "
              << static_cast<double>(clock() - c0) / CLOCKS_PER_SEC << std::endl;
    return 0;
}