代码之家 › 专栏 › 技术社区 › BlueLemon

为什么这个版本的矩阵副本这么慢?

julia copy compilation performance

5

BlueLemon · 技术社区 · 8 年前

我注意到朱莉娅在矩阵复制时有一种奇怪的行为。

考虑以下三个功能:

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 .

此外,对于第一个函数,垃圾收集器消耗的时间变化很大。如果不执行垃圾收集,则较大的差异会变小,但仍然存在。在第3版和第2版之间仍然是2.5倍。

那么为什么“memcopy”版本不如“assignment”版本高效呢?

2 回复 | 直到 8 年前

1

5

Colin T Bowers 8 年前

首先,您的代码包含一个bug。运行以下命令:

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

即使修复了bug,您仍然会注意到 虽然 循环解决方案更快。这是因为julia中的数组切片创建临时数组。所以在这一行:

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

右侧操作创建一个临时nxn数组,然后将临时数组分配给左侧。

如果希望避免临时数组分配,可以改为写:

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

您会注意到这两种方法的计时非常接近,尽管 虽然 循环速度仍然略快。一般来说,Julia中的循环速度很快(就像C fast一样),显式写出循环通常会得到最优化的编译代码。我仍然希望显式循环比 view 方法

至于垃圾收集,这只是你计时方法的结果。使用起来更好 @btime 从包裹里 BenchmarkTools ,它使用各种技巧来避免陷阱,如定时垃圾收集等。

2

3

BlueLemon 8 年前

为什么 A[1:n,1:n] = view(B, 1:n, 1:n) 或者它的变体,比一组while循环慢?让我们看看什么 A[1:n,1:n]=视图(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 指数呢 我 从…起 A .在每个步骤中 reindex 从视图计算 五、 指数呢 我 获取元素的索引 B .它叫 r 我们还了它。最后 R 已写入 A. .

每次复印后 inc 增加索引 我 到中的下一个元素 A. 并测试是否完成。请注意,代码来自v0。63但在 master 这或多或少是一样的。

原则上,代码可以简化为一组while循环,但它更通用。它适用于任意的视图 B 以及表格的任意切片 a:b:c 对于任意数量的矩阵维数。大人物 N 这就是我们的情况 2 .

由于函数更复杂,编译器也不会对它们进行优化。也就是说,有一个建议是编译器应该内联它们,但它没有这样做。这表明所显示的函数是非平凡的。

对于一组循环,编译器将最里面的循环减少为三个加法(每个加法用于指向 A. 和 B 一个用于循环索引)和一条复制指令。

tl;博士 企业内部的呼叫链 A[1:n,1:n]=视图(B,1:n,1:n) 再加上多个分派是非常重要的,可以处理一般情况。这会导致开销。一组while循环已经针对一种特殊情况进行了优化。

请注意,性能取决于编译器。如果我们看一维的情况 A[1:n] = view(B, 1:n) ,它比while循环更快,因为它将代码矢量化。但对于更高的维度 N >2 差别越来越大。