Почему эта версия матричной копии такая медленная? - PullRequest
Новая
       48

Почему эта версия матричной копии такая медленная?

0 голосов
/ 16 мая 2018

Я заметил странное поведение Джулии во время матричной копии.

Рассмотрим следующие три функции: 

function priv_memcopyBtoA!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  A[1:n,1:n] = B[1:n,1:n]
  return nothing
end

function priv_memcopyBtoA2!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  ii = 1; jj = 1;
  while ii <= n
    jj = 1 #(*)
    while jj <= n
      A[jj,ii] = B[jj,ii]
      jj += 1
    end
    ii += 1
  end
  return nothing
end

function priv_memcopyBtoA3!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  A[1:n,1:n] = view(B, 1:n, 1:n)
  return nothing
end

Редактировать : 1) Я проверял, будет ли код выдавать BoundsError, поэтому строка, отмеченная jj = 1 #(*), отсутствовала в исходном коде. Результаты тестирования были уже из фиксированной версии, поэтому они остаются без изменений. 2) Я добавил вариант просмотра, спасибо @Colin T Bowers за решение обеих проблем.

Кажется, что обе функции должны приводить к более или менее одинаковому коду. Все же я получаю за 

A = fill!(Matrix{Int32}(2^12,2^12),2); B = Int32.(eye(2^12));

результаты 

@timev priv_memcopyBtoA!(A,B, 2000)
  0.178327 seconds (10 allocations: 15.259 MiB, 85.52% gc time)
elapsed time (ns): 178326537
gc time (ns):      152511699
bytes allocated:   16000304
pool allocs:       9
malloc() calls:    1
GC pauses:         1

и 

@timev priv_memcopyBtoA2!(A,B, 2000)
 0.015760 seconds (4 allocations: 160 bytes)
elapsed time (ns): 15759742
bytes allocated:   160
pool allocs:       4

и 

@timev priv_memcopyBtoA3!(A,B, 2000)
  0.043771 seconds (7 allocations: 224 bytes)
elapsed time (ns): 43770978
bytes allocated:   224
pool allocs:       7

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

Вторая версия имеет издержки из арифметики указателя (getindex), условия ветвления (<=) и проверки границ в каждом назначении. Все же каждое назначение занимает всего ~3 ns.

Кроме того, время, которое потребляет сборщик мусора, сильно зависит от первой функции. Если сборка мусора не выполняется, большая разница становится небольшой, но она остается. Это все еще в ~ 2,5 раза между версией 3 и 2.

Так почему же версия "memcopy" не так эффективна, как версия "присваивания"?

Ответы [ 2 ]

0 голосов
/ 20 мая 2018

Почему A[1:n,1:n] = view(B, 1:n, 1:n) или его варианты медленнее, чем набор циклов while?Давайте посмотрим, что делает A[1:n,1:n] = view(B, 1:n, 1:n).

view возвращает итератор, который содержит указатель на родительский элемент B и информацию о том, как вычислять индексы, которые следует скопировать.A[1:n,1:n] = ... анализируется на вызов _setindex!(...).После этого, и несколько вызовов по цепочке вызовов, основная работа выполняется следующим образом: 

.\abstractarray.jl:883;
# In general, we simply re-index the parent indices by the provided ones
function getindex(V::SlowSubArray{T,N}, I::Vararg{Int,N}) where {T,N}
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds r = V.parent[reindex(V, V.indexes, I)...]
    r
end
#.\multidimensional.jl:212;
@inline function next(iter::CartesianRange{I}, state) where I<:CartesianIndex
    state, I(inc(state.I, iter.start.I, iter.stop.I))
end
@inline inc(::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
@inline inc(state::Tuple{Int}, start::Tuple{Int}, stop::Tuple{Int}) = (state[1]+1,)
@inline function inc(state, start, stop)
    if state[1] < stop[1]
        return (state[1]+1,tail(state)...)
    end
    newtail = inc(tail(state), tail(start), tail(stop))
    (start[1], newtail...)
end

getindex принимает представление V и индекс I.Мы получаем представление из B и индекс I из A.На каждом этапе reindex вычисляет из представления V и индексов индекса I, чтобы получить элемент в B.Он называется r и мы его возвращаем.Наконец, r записывается в A.

. После каждой копии inc увеличивает индекс I до следующего элемента в A и проверяет, выполнено ли оно.Обратите внимание, что код взят из v0.63, но в master он более или менее одинаков.

В принципе код можно сократить до набора циклов while, но он более общий.Он работает для произвольных представлений B и произвольных срезов вида a:b:c и для произвольного числа размеров матрицы.Большое N в нашем случае 2.

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

Для набора циклов компилятор сокращает самый внутренний цикл до трех сложений (каждый для указателя на A и B и один для индекса цикла)и одна инструкция копирования.

tl; dr Внутренняя цепочка вызовов A[1:n,1:n] = view(B, 1:n, 1:n) в сочетании с многократной диспетчеризацией нетривиальна и обрабатывает общий случай.Это вызывает накладные расходы.Набор циклов while уже оптимизирован для специального случая.

Обратите внимание, что производительность зависит от компилятора.Если взглянуть на одномерный случай A[1:n] = view(B, 1:n), это быстрее, чем цикл while, потому что он векторизует код.Но для более высоких измерений N >2 разница возрастает.

0 голосов
/ 17 мая 2018

Во-первых, ваш код содержит ошибку. Запустите это: 

A = [1 2 ; 3 4]
B = [5 6 ; 7 8]
priv_memcopyBtoA2!(A, B, 2)

тогда: 

julia> A
2×2 Array{Int64,2}:
 5  2
 7  4

Вам необходимо переназначить jj обратно на 1 в конце каждого внутреннего цикла while, то есть: 

function priv_memcopyBtoA2!(A::Matrix{Int}, B::Matrix{Int}, n::Int)
  ii = 1 
  while ii <= n
    jj = 1
    while jj <= n
      A[jj,ii] = B[jj,ii]
      jj += 1
    end
    ii += 1
  end
  return nothing
end

Даже с исправлением ошибки вы все равно заметите, что решение для цикла while работает быстрее. Это связано с тем, что фрагменты массива в julia создают временные массивы. Итак, в этой строке: 

A[1:n,1:n] = B[1:n,1:n]

правая операция создает временный массив nxn, а затем назначает временный массив левой стороне.

Если вы хотите избежать выделения временного массива, вы должны написать: 

A[1:n,1:n] = view(B, 1:n, 1:n)

и вы заметите, что синхронизация двух методов теперь довольно близка, хотя цикл while все еще немного быстрее. Как правило, циклы в Julia быстрые (как в C fast), и явное написание цикла обычно дает вам наиболее оптимизированный скомпилированный код. Я все еще ожидал бы, что явный цикл будет быстрее, чем метод view.

Что касается сбора мусора, то это всего лишь результат вашего метода синхронизации. Гораздо лучше использовать @btime из пакета BenchmarkTools, который использует различные приемы, чтобы избежать ловушек, таких как сборка мусора и т. Д.

...