Как можно сортировать сортировку?

Ответ 1

Вот решение, использующее (тип экспериментального) модуля Base.Threads.

Решение с использованием pmap (и т.д.) для распределенного parallelism будет аналогичным. Хотя я думаю, что накладные расходы между процессами могут навредить вам.

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

Затем происходит объединение этих предварительно отсортированных блоков.

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

И не забудьте настроить мультипоточность, установив переменную окружения JULIA_NUM_THREADS, прежде чем начать.

Вот мой код:

using Base.Threads

function blockranges(nblocks, total_len)
    rem = total_len % nblocks
    main_len = div(total_len, nblocks)

    starts=Int[1]
    ends=Int[]
    for ii in 1:nblocks
        len = main_len
        if rem>0
            len+=1
            rem-=1
        end
        push!(ends, starts[end]+len-1)
        push!(starts, ends[end] + 1)
    end
    @assert ends[end] == total_len
    starts[1:end-1], ends
end

function threadedsort!(data::Vector)
    starts, ends = blockranges(nthreads(), length(data))

    # Sort each block
    @threads for (ss, ee) in collect(zip(starts, ends))
        @inbounds sort!(@view data[ss:ee])
    end


    # Go through each sorted block taking out the smallest item and putting it in the new array
    # This code could maybe be optimised. see https://stackoverflow.com/a/22057372/179081
    ret = similar(data) # main bit of allocation right here. avoiding it seems expensive.
    # Need to not overwrite data we haven't read yet
    @inbounds for ii in eachindex(ret)
        minblock_id = 1
        ret[ii]=data[starts[1]]
        @inbounds for blockid in 2:endof(starts) # findmin allocates a lot for some reason, so do the find by hand. (maybe use findmin! ?)
            ele = data[starts[blockid]]
            if ret[ii] > ele
                ret[ii] = ele
                minblock_id = blockid
            end
        end
        starts[minblock_id]+=1 # move the start point forward
        if starts[minblock_id] > ends[minblock_id]
            deleteat!(starts, minblock_id)
            deleteat!(ends, minblock_id)
        end
    end
    data.=ret  # copy back into orignal as we said we would do it inplace
    return data
end

Я провел несколько тестов:

using Plots
function evaluate_timing(range)
    sizes = Int[]
    threadsort_times = Float64[]
    sort_times = Float64[]
        for sz in 2.^collect(range)
            data_orig = rand(Int, sz)
            push!(sizes, sz)

            data = copy(data_orig)
            push!(sort_times,       @elapsed sort!(data))

            data = copy(data_orig)
            push!(threadsort_times, @elapsed threadedsort!(data))

            @show (sz, sort_times[end], threadsort_times[end])
    end
    return sizes, threadsort_times, sort_times
end

sizes, threadsort_times, sort_times = evaluate_timing(0:28)
plot(sizes, [threadsort_times sort_times]; title="Sorting Time", ylabel="time(s)", xlabel="number of elements", label=["threadsort!" "sort!"])
plot(sizes, [threadsort_times sort_times]; title="Sorting Time", ylabel="time(s)", xlabel="number of elements", label=["threadsort!" "sort!"], xscale=:log10, yscale=:log10)

Мои результаты: , используя 8 потоков.

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

Как ни странно, эти результаты не воспроизводятся при использовании BenchmarkTools. Инструменты проверки производительности остановили бы подсчет первоначального времени. Но они очень последовательно воспроизводятся при использовании нормального кода времени, как и в приведенном выше контрольном коде. Я предполагаю, что он делает что-то, что убивает многопоточность, как

Большое спасибо @xiaodai, который указал на ошибку в моем аналитическом коде

Ответ 2

Я еще тестировал, если только 1% элементов уникальны и также выборка из 1:1_000_000. Результаты ниже

  function valu_timing_w_repeats (диапазон)       размеры = Int []       threadsort_times = Float64 []       sort_times = Float64 []           для sz в 2. ^ собирать (диапазон)               data_orig = rand (rand (Int, sz ÷ 100), sz)               push (размеры, sz)

            data = copy(data_orig)
            push!(sort_times,       @elapsed sort!(data))

            data = copy(data_orig)
            push!(threadsort_times, @elapsed threadedsort!(data))

            @show (sz, sort_times[end], threadsort_times[end])
    end
    return sizes, threadsort_times, sort_times
end

sizes, threadsort_times, sort_times = evaluate_timing_w_repeats(7:28)
plot(sizes, [threadsort_times sort_times]; title="Sorting Time", ylabel="time(s)", xlabel="number of elements", label=["threadsort!" "sort!"])
plot(sizes, [threadsort_times sort_times]; title="Sorting Time", ylabel="time(s)", xlabel="number of elements", label=["threadsort!" "sort!"], xscale=:log10, yscale=:log10)
savefig("sort_with_repeats.png")

function evaluate_timing1m(range)
    sizes = Int[]
    threadsort_times = Float64[]
    sort_times = Float64[]
        for sz in 2.^collect(range)
            data_orig = rand(1:1_000_000, sz)
            push!(sizes, sz)

            data = copy(data_orig)
            push!(sort_times,       @elapsed sort!(data))

            data = copy(data_orig)
            push!(threadsort_times, @elapsed threadedsort!(data))

            @show (sz, sort_times[end], threadsort_times[end])
    end
    return sizes, threadsort_times, sort_times
end

sizes, threadsort_times, sort_times = evaluate_timing1m(7:28)
plot(sizes, [threadsort_times sort_times]; title="Sorting Time", ylabel="time(s)", xlabel="number of elements", label=["threadsort!" "sort!"])
plot(sizes, [threadsort_times sort_times]; title="Sorting Time sampel from 1:1_000_000", ylabel="time(s)", xlabel="number of elements", label=["threadsort!" "sort!"], xscale=:log10, yscale=:log10)
savefig("sort1m.png")