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

Question

在矩阵复制过程中,我注意到了julia的奇怪行为.

考虑以下三个功能:

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,对于大内存块来说很难被击败.

第二个版本具有指针arithmetic(getindex),分支条件(<=)和每个赋值中的边界检查的开销.然而,每项任务都需要~3 ns.

此外,垃圾收集器消耗的时间对于第一个功能而言变化很大.如果没有执行垃圾收集,则大的差异会变小,但仍然存在.在版本3和版本2之间,它仍然是~2.5倍.

那么为什么"memcopy"版本不如"赋值"版本效率高？

Answer 1

首先,您的代码包含一个错误.运行这个:

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中使用更好,它使用各种技巧来避免陷阱,如计时垃圾收集等.