Wettbewerb: schnellste Methode zum Sortieren einer großen Anzahl von Gauß-verteilten Daten

Angesichts des Interesses an dieser Frage hielt ich es für interessant, die Antworten durch einen Wettbewerbsvorschlag etwas objektiver und quantitativer zu gestalten.

Die Idee ist einfach: Ich habe eine Binärdatei generiert, die 50 Millionen Gauß-verteilte Double enthält (Durchschnitt: 0, stdev 1). Das Ziel ist es, ein Programm zu erstellen, das diese so schnell wie möglich im Speicher sortiert. Eine sehr einfache Referenzimplementierung in Python benötigt 1m4s. Wie tief können wir gehen?

Die Regeln lauten wie folgt: Beantworten Sie die Frage mit einem Programm, das die Datei "gaussian.dat" öffnet und die Nummern im Speicher sortiert (keine Ausgabe erforderlich) sowie Anweisungen zum Erstellen und Ausführen des Programms. Das Programm muss auf meinem Arch Linux-Computer ausgeführt werden können (dh Sie können eine beliebige Programmiersprache oder Bibliothek verwenden, die auf diesem System leicht zu installieren ist).

Das Programm muss einigermaßen lesbar sein, damit ich sicherstellen kann, dass es sicher gestartet werden kann (bitte keine Nur-Assembler-Lösung!).

Ich werde die Antworten auf meinem Computer ausführen (Quadcore, 4 Gigabyte RAM). Die schnellste Lösung erhält die akzeptierte Antwort und ein Kopfgeld von 100 Punkten :)

Das Programm zur Generierung der Zahlen:

#!/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)

Die einfache Referenzimplementierung:

#!/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)

EDIT: nur 4 GB RAM, sorry

EDIT # 2: Beachten Sie, dass der Punkt des Wettbewerbs ist, um zu sehen, ob wir vorherige Informationen über die Daten verwenden können . Es soll kein Pissing Match zwischen verschiedenen Programmiersprachen-Implementierungen sein!

Hier ist eine Lösung in C ++, bei der zuerst die Zahlen in Gruppen mit der gleichen erwarteten Anzahl von Elementen aufgeteilt werden und dann jede Gruppe separat sortiert wird. Es berechnet eine Tabelle der kumulativen Verteilungsfunktion anhand einiger Formeln aus Wikipedia vor und interpoliert dann Werte aus dieser Tabelle, um eine schnelle Annäherung zu erhalten.

In mehreren Threads werden mehrere Schritte ausgeführt, um die vier Kerne zu nutzen.

#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;
}

Verwenden Sie zum Kompilieren und Ausführen diesen Befehl:

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

BEARBEITEN: Alle Buckets werden jetzt in dasselbe Array gestellt, damit die Buckets nicht wieder in das Array kopiert werden müssen. Auch die Größe der Tabelle mit vorberechneten Werten wurde reduziert, da die Werte genau genug sind. Wenn ich jedoch die Anzahl der Buckets über 256 ändere, dauert die Ausführung des Programms länger als bei dieser Anzahl von Buckets.

EDIT: Gleicher Algorithmus, andere Programmiersprache. Ich habe C ++ anstelle von Java verwendet und die Laufzeit auf meinem Computer von ~ 3,2 s auf ~ 2,35 s reduziert. Die optimale Anzahl von Eimern liegt immer noch bei 256 (wieder auf meinem Computer).

Übrigens, tbb ist wirklich großartig.

EDIT: Ich war von Alexandru's großartiger Lösung inspiriert und habe die std :: sort in der letzten Phase durch eine modifizierte Version seiner radix Sort ersetzt. Ich habe eine andere Methode verwendet, um mit den positiven / negativen Zahlen umzugehen, obwohl mehr Durchgänge durch das Array erforderlich sind. Ich habe auch beschlossen, das Array genau zu sortieren und die Einfügesortierung zu entfernen. Ich werde später einige Zeit damit verbringen, zu testen, wie sich diese Änderungen auf die Leistung auswirken und sie möglicherweise rückgängig machen. Bei Verwendung der Radix-Sortierung verringerte sich die Zeit jedoch von ~ 2,35 s auf ~ 1,63 s.

Ohne schlau zu werden, nur um einen viel schnelleren naiven Sortierer bereitzustellen, ist hier einer in C, der Ihrem Python-Sortierer ziemlich ähnlich sein sollte:

#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;
}

Kompiliert mit gcc -O3, auf meinem Rechner dauert dies mehr als eine Minute weniger als auf dem Python: ca. 11 s im Vergleich zu 87 s.

Ich habe anhand der Standardabweichung in Segmente unterteilt, die es am besten in Vierteln aufteilen sollten. Bearbeiten: Auf Partition basierend auf dem x-Wert in http://en.wikipedia.org/wiki/Error_function#Table_of_values umgeschrieben

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

Ich habe versucht, kleinere Eimer zu verwenden, aber es schien nur eine geringe Wirkung zu haben, sobald die Anzahl der verfügbaren Kerne überschritten wurde. Ohne parallele Sammlungen würde es 37 Sekunden auf meiner Box und 24 Sekunden bei den parallelen Sammlungen dauern. Wenn Sie über eine Distribution partitionieren, können Sie nicht nur ein Array verwenden, es entsteht also ein zusätzlicher Aufwand. Mir ist nicht klar, wann ein Wert in Scala ein- oder ausgepackt wird.

Ich verwende Scala 2.9 für die parallele Sammlung. Sie können einfach die tar.gz-Distribution herunterladen.

Zum Kompilieren: scalac SortFile.scala (Ich habe es gerade direkt in den Ordner scala / bin kopiert.

Zum Ausführen: JAVA_OPTS = "- Xmx4096M" ./scala SortFile (Ich habe es mit 2 Gigs RAM ausgeführt und ungefähr zur gleichen Zeit erhalten)

Bearbeiten: allocateDirect wurde entfernt, langsamer als nur das Zuweisen. Priming der Anfangsgröße für Array-Puffer entfernt. Eigentlich hat es die ganzen 50000000 Werte gelesen. Neu geschrieben, um hoffentlich Autoboxing-Probleme zu vermeiden (immer noch langsamer als naiv 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());

  }
}

Schreiben Sie dies einfach in eine cs-Datei und kompilieren Sie es theoretisch mit csc: (Benötigt Mono)

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();
        }
    }
}

Da Sie die Verteilung kennen, können Sie eine direkte Indizierung nach O (N) verwenden. (Wenn Sie sich fragen, was das ist, nehmen Sie an, Sie haben ein Kartenspiel mit 52 Karten und möchten es sortieren. Haben Sie nur 52 Fächer und werfen Sie jede Karte in ihr eigenes Fach.)

Sie haben 5e7 Doppel. Ordnen Sie ein Ergebnisarray R mit 5e7-Doppelwerten zu. Nimm jede Zahl xund hol sie dir i = phi(x) * 5e7. Grundsätzlich tun R[i] = x. Sie können mit Kollisionen umgehen, z. B. indem Sie die Nummer verschieben, mit der sie möglicherweise kollidiert (wie bei der einfachen Hash-Codierung). Alternativ können Sie R ein paar Mal größer machen und mit einem eindeutigen leeren Wert füllen . Am Ende fegen Sie einfach die Elemente von R.

phiist nur die Gaußsche kumulative Verteilungsfunktion. Es konvertiert eine gaußsche verteilte Zahl zwischen +/- unendlich in eine gleichmäßige verteilte Zahl zwischen 0 und 1. Eine einfache Methode zur Berechnung ist das Nachschlagen und Interpolieren von Tabellen.

Hier ist eine andere sequentielle Lösung:

#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;
}

Ich bezweifle, dass es die Multithread-Lösung schlägt, aber die Timings auf meinem i7-Laptop sind (stdsort ist die C ++ - Lösung, die in einer anderen Antwort bereitgestellt wird):

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

Beachten Sie, dass diese Lösung eine lineare Zeitkomplexität aufweist (da die spezielle Darstellung von Doppelwerten verwendet wird).

BEARBEITEN : Die Reihenfolge der Elemente wurde korrigiert.

EDIT : Geschwindigkeit um fast eine halbe Sekunde verbessert.

BEARBEITEN : Geschwindigkeit um weitere 0,7 Sekunden verbessert. Der Algorithmus wurde cachefreundlicher.

EDIT : Geschwindigkeit um 1 Sekunde erhöht. Da es nur 50.000.000 Elemente gibt, kann ich die Mantisse teilweise sortieren und Insert-Sort (das cachefreundlich ist) verwenden, um fehl am Platze liegende Elemente zu reparieren. Diese Idee entfernt ungefähr zwei Iterationen aus der letzten Radix-Sortierschleife.

EDIT : 0,16 Sekunden weniger. First std :: reverse kann bei umgekehrter Sortierreihenfolge entfallen.

Nehmen Sie die Lösung von Christian Ammer und parallelisieren Sie sie mit den Threaded Building Blocks von 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;
}

Wenn Sie Zugriff auf die IPP-Bibliothek (Performance Primitives) von Intel haben, können Sie deren Radix-Sortierung verwenden. Einfach austauschen

#include <tbb/parallel_sort.h>

mit

#include "ipps.h"

und

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

mit

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

Auf meinem Dual-Core-Laptop sind die Timings

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

Wie wäre es mit einer Implementierung von parallelem QuickSort , bei der die Pivot-Werte basierend auf den Statistiken der Verteilung ausgewählt werden und dabei gleich große Partitionen sichergestellt werden? Der erste Pivot wäre der Mittelwert (in diesem Fall Null), das nächste Paar wäre das 25. und das 75. Perzentil (+/- 0,67449 Standardabweichungen) usw., wobei jede Partition den verbleibenden Datensatz mehr oder halbiert weniger perfekt.

Sehr hässlich (warum Arrays verwenden, wenn ich Variablen verwenden kann, die mit Zahlen enden), aber schneller Code (mein erster Versuch zu std :: threads), ganze Zeit (Zeit real) auf meinem System 1,8 s (im Vergleich zu std :: sort) () 4,8 s), kompiliere mit g ++ -std = c ++ 0x -O3 -march = native -pthread Übergebe einfach Daten über stdin (funktioniert nur für 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);
}

// Bearbeiten geändert, um die Datei gaussian.dat zu lesen.

Eine C ++ - Lösung mit std::sort(möglicherweise schneller als qsort, in Bezug auf die Leistung von qsort gegenüber 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;
}

Ich kann nicht zuverlässig sagen, wie lange es dauert, da ich nur 1 GB auf meinem Computer habe und mit dem angegebenen Python-Code nur eine gaussian.datDatei mit nur 25 Millionen Doppeln erstellen konnte (ohne einen Speicherfehler zu erhalten). Aber ich bin sehr interessiert, wie lange der std :: sort-Algorithmus läuft.

Hier ist eine Mischung aus Alexandru's Radix-Sorte mit Zjareks gewundenem Smart-Pivoting. Kompiliere es mit

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

Sie können die Radixgröße ändern, indem Sie STEP definieren (z. B. -DSTEP = 11 hinzufügen). Ich fand das Beste für meinen Laptop 8 (die Standardeinstellung).

Standardmäßig wird das Problem in vier Teile aufgeteilt und auf mehreren Threads ausgeführt. Sie können dies ändern, indem Sie einen Tiefenparameter an die Befehlszeile übergeben. Also, wenn Sie zwei Kerne haben, führen Sie es als

sorter_gaussian_radix 50000000 1

und wenn du 16 Kerne hast

sorter_gaussian_radix 50000000 4

Die maximale Tiefe beträgt derzeit 6 (64 Fäden). Wenn Sie zu viele Ebenen setzen, verlangsamen Sie den Code.

Eine Sache, die ich auch ausprobiert habe, war die Radix-Sortierung aus der Intel Performance Primitives (IPP) -Bibliothek. Die Implementierung von Alexandru stützt IPP auf solide Weise, wobei IPP etwa 30% langsamer ist. Diese Variante ist auch hier enthalten (auskommentiert).

#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;
}

BEARBEITEN : Ich habe die Cache-Verbesserungen von Alexandru implementiert und das hat ungefähr 30% der Zeit auf meinem Computer gekürzt.

BEARBEITEN : Dies implementiert eine rekursive Sortierung, so dass es auf Alexandru's 16-Kern-Maschine gut funktionieren sollte. Es benutzt auch Alexandru's letzte Verbesserung und entfernt eine der Umkehrungen. Für mich ergab sich eine Verbesserung von 20%.

BEARBEITEN : Es wurde ein Vorzeichenfehler behoben, der zu Ineffizienz führte, wenn mehr als 2 Kerne vorhanden waren.

BEARBEITEN : Lambda wurde entfernt, so dass es mit älteren Versionen von gcc kompiliert werden kann. Es enthält die auskommentierte IPP-Codevariante. Ich habe auch die Dokumentation für das Laufen auf 16 Kernen korrigiert. Soweit ich das beurteilen kann, ist dies die schnellste Implementierung.

BEARBEITEN : Ein Fehler wurde behoben, wenn STEP nicht 8 ist. Die maximale Anzahl der Threads wurde auf 64 erhöht. Einige Timing-Informationen wurden hinzugefügt.

Ich denke, das hängt wirklich davon ab, was Sie tun möchten. Wenn Sie eine Reihe von Gaußschen sortieren möchten, hilft Ihnen das nicht weiter. Aber wenn Sie einen Haufen sortierter Gaußscher wollen, ist dies der Fall. Auch wenn dies das Problem ein wenig verfehlt, halte ich es für interessant, die tatsächlichen Sortierroutinen mit denen zu vergleichen.

Wenn Sie schnell sein möchten, tun Sie weniger.

Anstatt eine Reihe von Zufallsstichproben aus der Normalverteilung zu generieren und anschließend zu sortieren, können Sie eine Reihe von Stichproben aus der Normalverteilung in sortierter Reihenfolge generieren.

Sie können die Lösung hier verwenden , um n einheitliche Zufallszahlen in sortierter Reihenfolge zu generieren. Dann können Sie die inverse cdf (scipy.stats.norm.ppf) der Normalverteilung verwenden, um die einheitlichen Zufallszahlen durch inverse Transformationsabtastung in Zahlen aus der Normalverteilung umzuwandeln .

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))

Wenn Sie Ihre Hände schmutziger machen möchten, können Sie die vielen inversen cdf-Berechnungen möglicherweise beschleunigen, indem Sie eine iterative Methode verwenden und das vorherige Ergebnis als erste Schätzung verwenden. Da die Vermutungen sehr nahe beieinander liegen werden, erhalten Sie mit einer einzigen Iteration wahrscheinlich eine große Genauigkeit.

Probieren Sie diese wechselnde Guvante-Lösung mit dieser Main () aus. Sie beginnt zu sortieren, sobald das 1/4 IO-Lesen abgeschlossen ist. In meinem Test ist sie schneller:

    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();
    }
}

Da Sie die Verteilung kennen, wäre meine Idee, k Buckets mit jeweils der gleichen erwarteten Anzahl von Elementen zu erstellen (da Sie die Verteilung kennen, können Sie diese berechnen). Fegen Sie dann in O (n) -Zeit das Array und legen Sie die Elemente in ihre Eimer.

Sortieren Sie dann gleichzeitig die Eimer. Angenommen, Sie haben k Eimer und n Elemente. Das Sortieren eines Eimers dauert (n / k) lg (n / k). Angenommen, Sie haben p Prozessoren, die Sie verwenden können. Da Eimer unabhängig voneinander sortiert werden können, müssen Sie einen Multiplikator für die Obergrenze (k / p) festlegen. Dies ergibt eine endgültige Laufzeit von n + ceil (k / p) * (n / k) lg (n / k), die viel schneller sein sollte als n lg n, wenn Sie k gut wählen.

Eine einfache Optimierungsidee besteht darin, zwei Double in ein SSE-Register einzufügen, sodass jeder Thread mit zwei Elementen gleichzeitig arbeiten würde. Dies kann für einige Algorithmen kompliziert sein.

Sie können das Array auch in cachefreundliche Blöcke sortieren und die Ergebnisse dann zusammenführen. Es sollten zwei Ebenen verwendet werden: zum Beispiel zuerst 4 KB für L1 und dann 64 KB für L2.

Dies sollte sehr cachefreundlich sein, da die Bucket-Sortierung nicht über den Cache hinausgeht und die endgültige Zusammenführung den Speicher sequentiell durchläuft.

Heutzutage ist die Berechnung viel billiger als Speicherzugriffe. Wir haben jedoch eine große Anzahl von Elementen, sodass es schwierig ist, die Arraygröße zu bestimmen, wenn die dumme cachebewusste Sortierung langsamer ist als eine nicht cachebewusste Version mit geringer Komplexität.

Aber ich werde keine Implementierung des oben genannten bereitstellen, da ich es in Windows (VC ++) tun würde.

Hier ist eine lineare Scan-Bucket-Sortierung. Ich denke, es ist schneller als alle aktuellen Single-Thread-Implementierungen mit Ausnahme der Radix-Sortierung. Die erwartete Laufzeit sollte linear sein, wenn ich die PDF-Datei genau genug einschätze (ich verwende die lineare Interpolation von Werten, die ich im Web gefunden habe) und keine Fehler gemacht habe, die zu übermäßigem Scannen führen würden:

#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;
}

Ich weiß nicht, warum ich meinen vorherigen Beitrag nicht bearbeiten kann. Deshalb gibt es hier eine neue Version, die 0,2 Sekunden schneller ist (aber ungefähr 1,5 Sekunden schneller in der CPU-Zeit (Benutzer)). Diese Lösung hat 2 Programme, berechnet zuerst Quantile für die Normalverteilung für die Bucket-Sortierung vor und speichert sie in einer Tabelle, t [double * scale] = Bucket-Index, wobei scale eine willkürliche Zahl ist, die das Casting auf double ermöglicht. Das Hauptprogramm kann diese Daten dann verwenden, um Doppelsätze in den richtigen Eimer zu legen. Es hat einen Nachteil: Wenn die Daten nicht gaußsch sind, funktionieren sie nicht richtig (und es gibt auch fast keine Chance, dass sie bei normaler Verteilung falsch funktionieren), aber die Änderung für Sonderfälle ist einfach und schnell (nur die Anzahl der Buckets wird überprüft und fällt auf std.) ::Sortieren()).

Kompilieren: g ++ => Hilfsprogramm http://pastebin.com/WG7pZEzH

g ++ -std = c ++ 0x -O3 -march = native -pthread => http://pastebin.com/T3yzViZP Hauptsortierprogramm

Hier ist eine andere sequentielle Lösung. Dieser nutzt die Tatsache, dass die Elemente normalverteilt sind, und ich denke, die Idee ist allgemein anwendbar, um eine Sortierung nahe der linearen Zeit zu erreichen.

Der Algorithmus sieht folgendermaßen aus:

• Ungefähre CDF (siehe phi()Funktion in der Implementierung)
• Berechnen Sie für alle Elemente die ungefähre Position im sortierten Array: size * phi(x)
• Platzieren Sie die Elemente in einem neuen Array nahe ihrer endgültigen Position
  • In meiner Implementierung weist das Ziel-Array einige Lücken auf, damit ich beim Einfügen nicht zu viele Elemente verschieben muss.
• Verwenden Sie Insertsort, um die endgültigen Elemente zu sortieren (Insertsort ist linear, wenn der Abstand zur endgültigen Position kleiner als eine Konstante ist).

Leider ist die versteckte Konstante ziemlich groß und diese Lösung ist doppelt so langsam wie der Radix-Sortieralgorithmus.

Mein persönlicher Favorit unter Verwendung der Threaded Building Blocks von Intel wurde bereits veröffentlicht, aber hier ist eine einfache parallele Lösung unter Verwendung von JDK 7 und seiner neuen Fork / Join-API:

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;
            }
        }
    }    
}

Wichtiger Haftungsausschluss : Ich habe die Quick-Sort-Anpassung für fork / join von https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel übernommen

Um dies auszuführen, benötigen Sie eine Betaversion von JDK 7 (http://jdk7.java.net/download.html).

Auf meinem 2,93 GHz Quad Core i7 (OS X):

Python-Referenz

time python sort.py 50000000
sorting...

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

Java JDK 7 Fork / Join

time java ForkJoinQuicksortTask 50000000
Starting parallel computation

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

Ich habe auch versucht, mit parallelem Lesen und Konvertieren der Bytes in Doppelbytes zu experimentieren, aber ich habe dort keinen Unterschied festgestellt.

Aktualisieren:

Wenn jemand mit dem parallelen Laden der Daten experimentieren möchte, finden Sie unten die Version zum parallelen Laden. Theoretisch könnte dies dazu führen, dass es noch ein bisschen schneller geht, wenn Ihr IO-Device über genügend parallele Kapazität verfügt (SSDs normalerweise). Das Erstellen von Doubles aus Bytes ist außerdem mit einem gewissen Mehraufwand verbunden, sodass dies möglicherweise auch parallel schneller gehen kann. Auf meinen Systemen (Ubuntu 10.10 / Nehalem Quad / Intel X25M SSD und OS X 10.6 / i7 Quad / Samsung SSD) habe ich keinen wirklichen Unterschied festgestellt.

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:

Ich habe den Code auf einem unserer 12-Kern-Entwicklungscomputer mit einer geringfügigen Änderung ausgeführt, um eine feste Anzahl von Kernen festzulegen. Dies ergab die folgenden Ergebnisse:

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

Auf diesem System habe ich auch die Python-Version mit 1m2.994s und die C ++ - Version von Zjarek mit 1.925s ausprobiert (aus irgendeinem Grund scheint die C ++ - Version von Zjarek auf dem Computer von static_rtti relativ schneller zu laufen).

Ich habe auch versucht, was passiert ist, wenn ich die Dateigröße auf 100.000.000 verdoppelt habe:

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

In diesem Fall hat Zjareks C ++ - Version 3.968s gedauert. Python hat hier einfach zu lange gedauert.

150.000.000 Doppelbetten:

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

In diesem Fall war Zjareks C ++ - Version 6.044s. Ich habe es nicht einmal mit Python versucht.

Die C ++ - Version ist sehr konsistent mit ihren Ergebnissen, bei denen Java ein wenig schwankt. Zuerst wird es ein bisschen effizienter, wenn das Problem größer wird, aber dann wieder weniger effizient.

Eine Version mit traditionellen pthreads. Code zum Zusammenführen aus Guvantes Antwort kopiert. Kompilieren mit 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;
}  

Auf meinem Laptop erhalte ich folgende Ergebnisse:

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

Hier ist eine sequentielle C99-Implementierung, die versucht, die bekannte Distribution wirklich zu nutzen. Grundsätzlich wird eine einzelne Runde der Bucket-Sortierung anhand der Verteilungsinformationen durchgeführt. Anschließend werden einige Runden der Quick-Sortierung für jeden Bucket durchgeführt, wobei eine gleichmäßige Verteilung innerhalb der Grenzen des Bucket und schließlich eine geänderte Auswahlsortierung vorausgesetzt werden, um die Daten zurück in den ursprünglichen Puffer zu kopieren. Die Schnellsortierung speichert die Teilungspunkte, sodass die Auswahlsortierung nur für kleine Gruppen ausgeführt werden muss. Und trotz all dieser Komplexität (wegen?) Ist es nicht einmal wirklich schnell.

Um Φ schnell auswerten zu können, werden die Werte in wenigen Punkten abgetastet und später wird nur eine lineare Interpolation verwendet. Es ist eigentlich egal, ob Φ genau ausgewertet wird, solange die Approximation streng monoton ist.

Die Behältergrößen werden so gewählt, dass die Wahrscheinlichkeit eines Behälterüberlaufs vernachlässigbar ist. Genauer gesagt beträgt die Wahrscheinlichkeit, dass ein Datensatz mit 50000000 Elementen einen Bin-Überlauf verursacht, bei den aktuellen Parametern 3,65e-09. (Dies kann mit der Überlebensfunktion der Poisson-Verteilung berechnet werden .)

Zum Kompilieren bitte verwenden

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

Da es erheblich mehr Berechnungen gibt als in den anderen Lösungen, werden diese Compiler-Flags benötigt, um es zumindest einigermaßen schnell zu machen. Ohne -msse3die Konvertierungen von werden wirklich langsam. Wenn Ihre Architektur SSE3 nicht unterstützt, können diese Konvertierungen auch mit der Funktion durchgeführt werden.doubleintlrint()

Der Code ist ziemlich hässlich - ich bin mir nicht sicher, ob dies die Anforderung erfüllt, "einigermaßen lesbar" zu sein ...

#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;
}

Dies verwendet erf (), um jedes Element entsprechend in eine Bin zu platzieren, und sortiert dann jede Bin. Es hält das Array vollständig an Ort und Stelle.

Erster Durchgang: docensus () zählt die Anzahl der Elemente in jeder Bin.

Zweiter Durchgang: partition () durchläuft das Array und platziert jedes Element in seinem richtigen Bin

Dritter Durchgang: sortbins () führt eine Q-Sortierung für jeden Bin durch.

Es ist ein bisschen naiv und ruft die teure erf () - Funktion zweimal für jeden Wert auf. Der erste und dritte Durchgang sind möglicherweise parallelisierbar. Die zweite ist sehr seriell und wird wahrscheinlich durch ihre sehr zufälligen Speicherzugriffsmuster verlangsamt. Abhängig vom Verhältnis von CPU-Leistung zu Speichergeschwindigkeit kann es auch sinnvoll sein, die Bin-Nummern der einzelnen Double zwischenzuspeichern.

Mit diesem Programm können Sie die Anzahl der zu verwendenden Fächer auswählen. Fügen Sie einfach eine zweite Zahl in die Befehlszeile ein. Ich habe es mit gcc -O3 kompiliert, aber meine Maschine ist so schwach, dass ich Ihnen keine guten Leistungszahlen nennen kann.

Bearbeiten: Poof! Mein C-Programm hat sich mit std :: sort auf magische Weise in ein C ++ - Programm verwandelt!

Schauen Sie sich die radix sort Implementierung von Michael Herf ( Radix Tricks ) an. Auf meiner Maschine war die Sortierung im Vergleich zum std::sortAlgorithmus in meiner ersten Antwort fünfmal schneller . Der Name der Sortierfunktion lautet 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;
}