Kontes: cara tercepat untuk mengurutkan array besar data yang didistribusikan Gaussian

Mengikuti minat pada pertanyaan ini , saya pikir akan menarik untuk membuat jawaban yang sedikit lebih objektif dan kuantitatif dengan mengusulkan sebuah kontes.

Idenya sederhana: Saya telah menghasilkan file biner yang berisi 50 juta gaussian ganda yang didistribusikan (rata-rata: 0, stdev 1). Tujuannya adalah membuat program yang akan mengurutkan ini dalam memori secepat mungkin. Implementasi referensi yang sangat sederhana dalam python membutuhkan 1m4s untuk selesai. Seberapa rendah kita bisa pergi?

Aturannya adalah sebagai berikut: jawab dengan program yang membuka file "gaussian.dat" dan urutkan angka dalam memori (tidak perlu untuk menampilkannya), dan instruksi untuk membangun dan menjalankan program. Program ini harus dapat bekerja pada mesin Arch Linux saya (artinya Anda dapat menggunakan bahasa pemrograman atau pustaka yang mudah diinstal pada sistem ini).

Program ini harus mudah dibaca, sehingga saya bisa memastikannya aman untuk diluncurkan (tolong jangan ada solusi assembler saja!).

Saya akan menjalankan jawaban di mesin saya (quad core, 4 Gigabytes RAM). Solusi tercepat akan mendapatkan jawaban yang diterima dan hadiah 100 poin :)

Program yang digunakan untuk menghasilkan angka:

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

Implementasi referensi sederhana:

#!/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: hanya 4 GB RAM, maaf

EDIT # 2: Perhatikan bahwa tujuan dari kontes ini adalah untuk melihat apakah kita dapat menggunakan informasi sebelumnya tentang data . itu tidak seharusnya menjadi pertandingan kencing antara implementasi bahasa pemrograman yang berbeda!

Berikut ini adalah solusi dalam C ++ yang pertama-tama mempartisi angka-angka ke dalam ember dengan jumlah elemen yang diharapkan yang sama dan kemudian mengurutkan masing-masing ember secara terpisah. Itu precomputes tabel dari fungsi distribusi kumulatif berdasarkan pada beberapa formula dari Wikipedia dan kemudian interpolasi nilai dari tabel ini untuk mendapatkan perkiraan cepat.

Beberapa langkah dijalankan dalam beberapa utas untuk memanfaatkan keempat inti.

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

Untuk mengkompilasi dan menjalankannya, gunakan perintah ini:

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

EDIT: Semua ember sekarang ditempatkan ke dalam array yang sama untuk menghilangkan kebutuhan untuk menyalin kembali ember ke dalam array. Juga ukuran tabel dengan nilai yang dikomputasi berkurang, karena nilainya cukup akurat. Namun, jika saya mengubah jumlah ember di atas 256, program ini membutuhkan waktu lebih lama untuk dijalankan daripada dengan jumlah ember itu.

EDIT: Algoritma yang sama, bahasa pemrograman yang berbeda. Saya menggunakan C ++ sebagai ganti Java dan waktu operasi berkurang dari ~ 3.2s menjadi ~ 2.35s pada mesin saya. Jumlah bucket optimal masih sekitar 256 (sekali lagi, di komputer saya).

Omong-omong, tbb benar-benar hebat.

EDIT: Saya terinspirasi oleh solusi hebat Alexandru dan menggantikan std :: sort pada fase terakhir dengan versi modifikasi dari jenis radix-nya. Saya memang menggunakan metode yang berbeda untuk menangani angka positif / negatif, meskipun perlu melewati array lebih banyak. Saya juga memutuskan untuk mengurutkan array dengan tepat dan menghapus jenis penyisipan. Saya nantinya akan meluangkan waktu menguji bagaimana perubahan ini mempengaruhi kinerja dan mungkin mengembalikannya. Namun, dengan menggunakan jenis radix, waktu berkurang dari ~ 2,35 ke ~ 1,63.

Tanpa menjadi pintar, hanya untuk menyediakan penyortir naif yang jauh lebih cepat, inilah satu di C yang seharusnya cukup banyak setara dengan Python Anda:

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

Dikompilasi dengan gcc -O3, pada mesin saya ini membutuhkan waktu lebih dari satu menit kurang dari Python: sekitar 11 detik dibandingkan dengan 87 detik.

Saya mempartisi ke dalam segmen berdasarkan standar deviasi yang harus dipecah menjadi 4s. Sunting: Ditulis ulang ke partisi berdasarkan nilai x di http://en.wikipedia.org/wiki/Error_function#Table_of_values

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

Saya mencoba menggunakan ember yang lebih kecil, tetapi tampaknya memiliki sedikit efek 2 * melebihi jumlah core yang tersedia. Tanpa koleksi paralel, akan butuh 37 detik di kotak saya dan 24 dengan koleksi paralel. Jika mempartisi melalui distribusi, Anda tidak bisa hanya menggunakan array, jadi ada beberapa overhead lagi. Saya tidak jelas kapan nilai akan kotak / unboxed di scala.

Saya menggunakan scala 2.9, untuk koleksi paralel. Anda bisa mengunduh distribusi tar.gz itu.

Untuk mengkompilasi: scalac SortFile.scala (Saya baru saja menyalinnya langsung ke folder scala / bin.

Untuk menjalankan: JAVA_OPTS = "- Xmx4096M" ./scala SortFile (Saya menjalankannya dengan 2 pertunjukan ram dan mendapatkan waktu yang hampir bersamaan)

Sunting: Dihapus mengalokasikan Direct, lebih lambat dari hanya mengalokasikan. Priming yang dihapus dari ukuran awal untuk buffer array. Sebenarnya membuatnya membaca seluruh nilai 50000000. Menulis ulang untuk mudah-mudahan menghindari masalah autoboxing (masih lebih lambat dari naif 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());

  }
}

Masukkan saja ini ke dalam file cs dan kompilasi dengan teori csc: (Membutuhkan 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();
        }
    }
}

Karena Anda tahu apa distribusinya, Anda bisa menggunakan pengindeksan O (N) langsung. (Jika Anda bertanya-tanya apa itu, anggaplah Anda memiliki setumpuk 52 kartu dan Anda ingin menyortirnya. Hanya memiliki 52 nampan dan melemparkan setiap kartu ke nampan itu sendiri.)

Anda memiliki 5e7 ganda. Alokasikan array hasil R dari 5e7 ganda. Ambil setiap nomor xdan dapatkan i = phi(x) * 5e7. Pada dasarnya lakukan R[i] = x. Memiliki cara untuk menangani tabrakan, seperti memindahkan nomor yang mungkin bertabrakan dengannya (seperti dalam pengkodean hash sederhana). Atau, Anda bisa membuat R beberapa kali lebih besar, diisi dengan nilai kosong yang unik . Pada akhirnya, Anda hanya menyapu unsur-unsur R.

phihanyalah fungsi distribusi kumulatif gaussian. Ini mengubah angka terdistribusi gaussian antara +/- tak terhingga menjadi angka terdistribusi seragam antara 0 dan 1. Cara sederhana untuk menghitungnya adalah dengan pencarian tabel dan interpolasi.

Berikut ini adalah solusi berurutan lainnya:

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

Saya ragu itu mengalahkan solusi multi-threaded, tetapi timing pada laptop i7 saya (stdsort adalah solusi C ++ yang disediakan dalam jawaban lain):

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

Perhatikan bahwa solusi ini memiliki kompleksitas waktu linier (karena menggunakan representasi khusus ganda).

EDIT : Memperbaiki urutan elemen yang akan meningkat.

EDIT : Peningkatan kecepatan hampir setengah detik.

EDIT : Peningkatan kecepatan oleh 0,7 detik lainnya. Jadikan algoritma lebih ramah cache.

EDIT : Peningkatan kecepatan 1 detik lagi. Karena hanya ada 50.000.000 elemen yang saya dapat mengurutkan sebagian mantissa dan menggunakan menyortir jenis (yang ramah cache) untuk memperbaiki elemen out-of-place. Ide ini menghilangkan sekitar dua iterasi dari loop penyortiran radix terakhir.

EDIT : 0,16 lebih sedikit detik. Pertama std :: membalikkan bisa dihilangkan jika urutan penyortiran terbalik.

Mengambil solusi Christian Ammer dan memparalelkannya dengan Blok Bangunan Berulir 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;
}

Jika Anda memiliki akses ke perpustakaan Intel Performance Primitives (IPP), Anda dapat menggunakan jenis radix-nya. Ganti saja

#include <tbb/parallel_sort.h>

dengan

#include "ipps.h"

dan

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

dengan

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

Pada laptop dual core saya, waktunya adalah

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

Bagaimana dengan penerapan quicksort paralel yang memilih nilai pivot berdasarkan statistik distribusi, sehingga memastikan partisi berukuran sama? Poros pertama akan berada pada rata-rata (nol dalam kasus ini), pasangan berikutnya akan berada pada persentil ke-25 dan ke-75 (+/- -0,67449 standar deviasi), dan seterusnya, dengan masing-masing partisi mengurangi separuh dari set data yang tersisa lebih atau kurang sempurna.

Sangat jelek (mengapa menggunakan array ketika saya bisa menggunakan variabel yang diakhiri dengan angka), tetapi kode cepat (percobaan pertama saya ke std :: threads), sepanjang waktu (real time) pada sistem saya 1,8 s (dibandingkan dengan std :: sort () 4,8 dtk), kompilasi dengan g ++ -std = c ++ 0x -O3 -march = asli -pthread Cukup lewat data melalui stdin (hanya berfungsi untuk 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);
}

// Edit diubah untuk membaca file gaussian.dat.

Solusi C ++ menggunakan std::sort(akhirnya lebih cepat dari qsort, mengenai Kinerja 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;
}

Saya tidak dapat diandalkan untuk mengatakan berapa lama karena saya hanya memiliki 1GB di mesin saya dan dengan kode Python yang diberikan saya hanya bisa membuat gaussian.datfile dengan hanya 25mio ganda (tanpa mendapatkan kesalahan memori). Tapi saya sangat tertarik berapa lama algoritma std :: sort berjalan.

Berikut adalah campuran dari jenis radix Alexandru dengan pivot cerdas berulir Zjarek. Kompilasi dengan

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

Anda dapat mengubah ukuran radix dengan mendefinisikan LANGKAH (misalnya, tambahkan -DSTEP = 11). Saya menemukan yang terbaik untuk laptop saya adalah 8 (default).

Secara default, ini membagi masalah menjadi 4 bagian dan menjalankannya pada beberapa utas. Anda dapat mengubahnya dengan mengirimkan parameter kedalaman ke baris perintah. Jadi, jika Anda memiliki dua inti, jalankan sebagai

sorter_gaussian_radix 50000000 1

dan jika Anda memiliki 16 core

sorter_gaussian_radix 50000000 4

Kedalaman maks sekarang adalah 6 (64 utas). Jika Anda memasukkan terlalu banyak level, Anda hanya akan memperlambat kodenya.

Satu hal yang juga saya coba adalah radix sort dari perpustakaan Intel Performance Primitives (IPP). Implementasi Alexandru benar-benar merusak IPP, dengan IPP sekitar 30% lebih lambat. Variasi itu juga termasuk di sini (dikomentari).

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

EDIT : Saya menerapkan perbaikan cache Alexandru, dan itu mencukur sekitar 30% dari waktu di mesin saya.

EDIT : Ini mengimplementasikan semacam rekursif, sehingga harus bekerja dengan baik pada mesin 16 core Alexandru. Itu juga menggunakan perbaikan terakhir Alexandru dan menghapus salah satu yang sebaliknya. Bagi saya, ini memberikan peningkatan 20%.

EDIT : Memperbaiki bug tanda yang menyebabkan inefisiensi ketika ada lebih dari 2 core.

EDIT : Menghapus lambda, sehingga akan dikompilasi dengan versi gcc yang lebih lama. Ini termasuk variasi kode IPP yang dikomentari. Saya juga memperbaiki dokumentasi untuk berjalan pada 16 core. Sejauh yang saya tahu, ini adalah implementasi tercepat.

EDIT : Memperbaiki bug ketika LANGKAH tidak 8. Meningkatkan jumlah maksimum utas menjadi 64. Menambahkan beberapa info waktu.

Saya kira ini sangat tergantung pada apa yang ingin Anda lakukan. Jika Anda ingin mengurutkan sekelompok orang Gaussi, maka ini tidak akan membantu Anda. Tetapi jika Anda ingin sekelompok orang Gaussias yang diurutkan, ini akan. Bahkan jika ini sedikit merindukan masalahnya, saya pikir akan menarik untuk membandingkan vs rutinitas penyortiran yang sebenarnya.

Jika Anda ingin sesuatu menjadi cepat, lakukan lebih sedikit.

Alih-alih menghasilkan banyak sampel acak dari distribusi normal dan kemudian menyortir, Anda dapat menghasilkan banyak sampel dari distribusi normal dalam urutan diurutkan.

Anda dapat menggunakan solusi di sini untuk menghasilkan n nomor acak yang seragam dalam urutan diurutkan. Kemudian Anda dapat menggunakan invers cdf (scipy.stats.norm.ppf) dari distribusi normal untuk mengubah angka acak seragam menjadi angka dari distribusi normal melalui inverse transform sampling .

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

Jika Anda ingin mendapatkan tangan Anda lebih kotor, saya kira Anda mungkin bisa mempercepat banyak perhitungan cdf terbalik dengan menggunakan beberapa jenis metode berulang, dan menggunakan hasil sebelumnya sebagai tebakan awal Anda. Karena tebakannya akan sangat dekat, mungkin satu iterasi tunggal akan memberi Anda akurasi tinggi.

Coba ubah solusi Guvante ini dengan Main ini (), ia mulai mengurutkan begitu 1/4 IO selesai, lebih cepat dalam pengujian saya:

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

Karena Anda tahu distribusinya, ide saya adalah membuat k ember, masing-masing dengan jumlah elemen yang diharapkan (karena Anda tahu distribusinya, Anda dapat menghitung ini). Kemudian dalam waktu O (n), sapu array dan masukkan elemen ke dalam ember mereka.

Kemudian secara bersamaan menyortir ember. Misalkan Anda memiliki k k, dan n elemen. Sebuah ember akan mengambil (n / k) lg (n / k) waktu untuk menyortir. Sekarang anggaplah Anda memiliki prosesor p yang dapat Anda gunakan. Karena bucket dapat disortir secara independen, Anda memiliki pengali langit-langit (k / p) untuk ditangani. Ini memberikan runtime akhir n + ceil (k / p) * (n / k) lg (n / k), yang seharusnya jauh lebih cepat daripada n lg n jika Anda memilih k dengan baik.

Satu ide optimasi tingkat rendah adalah mencocokkan dua ganda dalam register SSE, sehingga setiap utas akan bekerja dengan dua item sekaligus. Ini mungkin rumit untuk dilakukan untuk beberapa algoritma.

Hal lain yang harus dilakukan adalah mengurutkan array dalam potongan cache-friendly, lalu menggabungkan hasilnya. Dua level harus digunakan: misalnya pertama 4 KB untuk L1 kemudian 64 KB untuk L2.

Ini harus sangat ramah-cache, karena jenis bucket tidak akan keluar dari cache, dan penggabungan akhir akan memunculkan memori secara berurutan.

Perhitungan hari ini jauh lebih murah daripada akses memori. Namun kami memiliki sejumlah besar item, jadi sulit untuk mengetahui mana ukuran array ketika jenis cache-aware lebih lambat daripada versi non-cache-aware dengan kompleksitas rendah.

Tapi saya tidak akan memberikan implementasi di atas karena saya akan melakukannya di Windows (VC ++).

Berikut adalah implementasi sortir ember pemindaian linier. Saya pikir ini lebih cepat dari semua implementasi single-threaded saat ini kecuali untuk jenis radix. Seharusnya linear waktu yang diharapkan berjalan jika saya memperkirakan cdf cukup akurat (Saya menggunakan interpolasi linear dari nilai yang saya temukan di web) dan tidak membuat kesalahan yang akan menyebabkan pemindaian berlebihan:

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

Saya tidak tahu, mengapa saya tidak bisa mengedit posting saya sebelumnya, jadi ini versi baru, 0,2 detik lebih cepat (tetapi sekitar 1,5 detik lebih cepat dalam waktu CPU (pengguna)). Solusi ini memiliki 2 program, pertama menghitung kuantil untuk distribusi normal untuk jenis bucket, dan menyimpannya dalam tabel, t [skala ganda *] = indeks bucket, di mana skala adalah beberapa angka acak yang memungkinkan casting untuk menggandakan kemungkinan. Kemudian program utama dapat menggunakan data ini untuk menempatkan ganda di ember yang benar. Ini memiliki satu kelemahan, jika data tidak gaussian itu tidak akan berfungsi dengan benar (dan juga hampir tidak ada peluang untuk bekerja secara salah untuk distribusi normal), tetapi modifikasi untuk case khusus mudah dan cepat (hanya jumlah bucket yang diperiksa dan jatuh ke std ::menyortir()).

Kompilasi: g ++ => http://pastebin.com/WG7pZEzH program bantuan

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

Berikut ini adalah solusi berurutan lainnya. Yang ini menggunakan fakta bahwa elemen-elemennya terdistribusi normal, dan saya pikir idenya secara umum berlaku untuk menyortir mendekati waktu linier.

Algoritmanya seperti ini:

• Perkiraan CDF (lihat phi()fungsi dalam implementasi)
• Untuk semua elemen hitung perkiraan posisi dalam array yang diurutkan: size * phi(x)
• Letakkan elemen dalam array baru dekat dengan posisi akhir mereka
  • Dalam array tujuan implementasi saya ada beberapa celah di dalamnya sehingga saya tidak perlu menggeser terlalu banyak elemen saat menyisipkan.
• Gunakan insertsort untuk mengurutkan elemen akhir (insertsort linier jika jarak ke posisi akhir lebih kecil dari konstanta).

Sayangnya, konstanta tersembunyi cukup besar dan solusi ini dua kali lebih lambat dari algoritma sortir radix.

Favorit pribadi saya menggunakan Blok Bangunan Berulir Intel telah diposting, tetapi berikut ini adalah solusi paralel paralel menggunakan JDK 7 dan fork / join API baru:

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

Penafian penting : Saya melakukan adaptasi cepat untuk fork / gabung dari: https://github.com/pmbauer/parallel/tree/master/src/main/java/pmbauer/parallel

Untuk menjalankan ini, Anda memerlukan versi beta JDK 7 (http://jdk7.java.net/download.html).

Pada core i7 2.93Ghz Quad core saya (OS X):

Referensi Python

time python sort.py 50000000
sorting...

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

Fork JDK 7 Java / gabung

time java ForkJoinQuicksortTask 50000000
Starting parallel computation

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

Saya juga mencoba melakukan beberapa percobaan dengan membaca paralel dan mengubah byte menjadi dua kali lipat, tetapi saya tidak melihat perbedaan di sana.

Memperbarui:

Jika ada yang ingin bereksperimen dengan pemuatan paralel data, versi pemuatan paralel di bawah. Secara teori ini bisa membuatnya sedikit lebih cepat, jika perangkat IO Anda memiliki kapasitas paralel yang cukup (SSD biasanya melakukannya). Ada juga beberapa overhead dalam membuat Doubles from bytes, sehingga berpotensi menjadi lebih cepat secara paralel juga. Pada sistem saya (SSD Ubuntu 10.10 / Nehalem Quad / Intel X25M, dan OS X 10.6 / i7 Quad / Samsung SSD) saya tidak melihat perbedaan nyata.

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

Pembaruan2:

Saya mengeksekusi kode pada salah satu dari 12 mesin dev inti kami dengan sedikit modifikasi untuk menetapkan jumlah inti yang tetap. Ini memberikan hasil sebagai berikut:

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

Pada sistem ini saya juga mencoba versi Python yang mengambil 1m2.994s dan versi C ++ Zjarek yang mengambil 1.925s (untuk beberapa alasan versi C ++ Zjarek tampaknya berjalan relatif lebih cepat di komputer static_rtti).

Saya juga mencoba apa yang terjadi jika saya menggandakan ukuran file menjadi 100.000.000 ganda:

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

Dalam hal ini, versi C ++ Zjarek mengambil 3,968s. Python terlalu lama di sini.

150.000.000 ganda:

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

Dalam hal ini, versi C ++ Zjarek adalah 6.044. Saya bahkan tidak mencoba Python.

Versi C ++ sangat konsisten dengan hasilnya, di mana Java sedikit berayun. Pertama itu menjadi sedikit lebih efisien ketika masalahnya menjadi lebih besar, tetapi kemudian kurang efisien lagi.

Versi menggunakan pthreads tradisional. Kode untuk menggabungkan disalin dari jawaban Guvante. Kompilasi dengan 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;
}  

Di laptop saya, saya mendapatkan hasil berikut:

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

Berikut ini adalah implementasi C99 berurutan yang mencoba untuk benar-benar memanfaatkan distribusi yang dikenal. Ini pada dasarnya melakukan satu putaran bucket sortir menggunakan informasi distribusi, kemudian beberapa putaran quicksort pada setiap bucket dengan asumsi distribusi seragam dalam batas bucket dan akhirnya sortir seleksi yang dimodifikasi untuk menyalin data kembali ke buffer asli. Quicksort menghafal titik perpecahan, jadi pemilihan semacam hanya perlu beroperasi pada potongan kecil. Dan terlepas dari (karena?) Semua kerumitan itu, itu bahkan tidak terlalu cepat.

Untuk membuat evaluasi Φ cepat, nilai sampel dalam beberapa poin dan kemudian hanya interpolasi linier yang digunakan. Sebenarnya tidak masalah jika Φ dievaluasi dengan tepat, asalkan perkiraannya benar-benar monoton.

Ukuran nampan dipilih sedemikian rupa sehingga kemungkinan nampan bin diabaikan. Lebih tepatnya, dengan parameter saat ini, kemungkinan dataset 50000000 elemen akan menyebabkan overflow bin adalah 3,65e-09. (Ini dapat dihitung menggunakan fungsi survival dari distribusi Poisson .)

Untuk mengkompilasi, silakan gunakan

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

Karena ada perhitungan yang jauh lebih banyak daripada solusi lain, flag-flag compiler ini diperlukan untuk membuatnya setidaknya cukup cepat. Tanpa -msse3konversi dari doubleke intmenjadi sangat lambat. Jika arsitektur Anda tidak mendukung SSE3, konversi ini juga dapat dilakukan menggunakan lrint()fungsi.

Kode ini agak jelek - tidak yakin apakah ini memenuhi persyaratan "cukup mudah dibaca" ...

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

Ini menggunakan erf () untuk menempatkan setiap elemen secara tepat ke dalam nampan, lalu mengurutkan masing-masing nampan. Itu membuat array sepenuhnya di tempat.

Pass pertama: docensus () menghitung jumlah elemen dalam setiap bin.

Lulus kedua: partisi () memungkinkan array, menempatkan setiap elemen ke dalam bin yang tepat

Lulus ketiga: sortbins () melakukan qsort pada setiap nampan.

Agak naif, dan memanggil fungsi erf () mahal dua kali untuk setiap nilai. Pass pertama dan ketiga berpotensi diparalelkan. Yang kedua sangat serial dan mungkin diperlambat oleh pola akses memori yang sangat acak. Mungkin juga bermanfaat untuk men-cache nomor bin setiap dobel, tergantung pada rasio kecepatan CPU terhadap memori.

Program ini memungkinkan Anda memilih jumlah sampah untuk digunakan. Tambahkan saja angka kedua ke baris perintah. Saya mengkompilasinya dengan gcc -O3, tetapi mesin saya sangat lemah sehingga saya tidak bisa memberi tahu Anda angka kinerja yang baik.

Edit: Poof! Program C saya secara ajaib berubah menjadi program C ++ menggunakan std :: sort!

Lihatlah implementasi radix sort oleh Michael Herf ( Radix Tricks ). Di mesin saya penyortiran 5 kali lebih cepat dibandingkan dengan std::sortalgoritma dalam jawaban pertama saya. Nama fungsi penyortiran adalah 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;
}