Julia 提供了一个基于消息传递的多处理器环境,能够同时在多处理器上使用独立的内存空间运行程序。
Julia 的消息传递与 MPI [1] 等环境不同。Julia 中的通信是“单边”的,即程序员只需要管理双处理器运算中的一个处理器即可。
Julia 中的并行编程基于两个原语:remote references 和 remote calls 。remote reference 对象,用于从任意的处理器,查阅指定处理器上存储的对象。 remote call 请求,用于一个处理器对另一个(也有可能是同一个)处理器调用某个函数处理某些参数。 remote call 返回 remote reference 对象。 remote call 是立即返回的;调用它的处理器继续执行下一步操作,而 remote call 继续在某处执行。可以对 remote reference 调用 wait
,以等待 remote call 执行完毕,然后通过 fetch
获取结果的完整值。使用 put
可将值存储到 remote reference 。
通过 julia -p n
启动,可以在本地机器上提供 n
个处理器。一般 n
等于机器上 CPU 内核个数:
$ ./julia -p 2 |
julia> r = remotecall(2, rand, 2, 2) |
RemoteRef(2,1,5) |
julia> fetch(r) |
2x2 Float64 Array: |
0.60401 0.501111 |
0.174572 0.157411 |
julia> s = @spawnat 2 1 .+ fetch(r) |
RemoteRef(2,1,7) |
julia> fetch(s) |
2x2 Float64 Array: |
1.60401 1.50111 |
1.17457 1.15741 |
remote_call
的第一个参数是要进行这个运算的处理器索引值。Julia 中大部分并行编程不查询特定的处理器或可用处理器的个数,但可认为 remote_call
是个为精细控制所提供的低级接口。第二个参数是要调用的函数,剩下的参数是该函数的参数。此例中,我们先让处理器 2 构造一个 2x2 的随机矩阵,然后我们在结果上加 1 。两个计算的结果保存在两个 remote reference 中,即 r
和 s
。 @spawnat
宏在由第一个参数指明的处理器上,计算第二个参数中的表达式。
remote_call_fetch
函数可以立即获取要在远端计算的值。它等价于 fetch(remote_call(...))
,但比之更高效:
julia> remotecall_fetch(2, getindex, r, 1, 1) |
0.10824216411304866 |
getindex(r,1,1)
:ref:等价于 <man-array-indexing>
r[1,1]
,因此,这个调用获取 remote reference 对象 r
的第一个元素。
remote_call
语法不太方便。 @spawn
宏简化了这件事儿,它对表达式而非函数进行操作,并自动选取在哪儿进行计算:
julia> r = @spawn rand(2,2) |
RemoteRef(1,1,0) |
julia> s = @spawn 1 .+ fetch(r) |
RemoteRef(1,1,1) |
julia> fetch(s) |
1.10824216411304866 1.13798233877923116 |
1.12376292706355074 1.18750497916607167 |
注意,此处用 1 .+ fetch(r)
而不是 1 .+ r
。这是因为我们不知道代码在何处运行,而 fetch
会将需要的 r
移到做加法的处理器上。此例中, @spawn
很聪明,它知道在有 r
对象的处理器上进行计算,因而 fetch
将不做任何操作。
( @spawn
不是内置函数,而是 Julia 定义的 :ref:宏 <man-macros>
)
所有执行程序代码的处理器上,都必须能获得程序代码。例如,输入:
julia> function rand2(dims...) |
return 2*rand(dims...) |
end |
julia> rand2(2,2) |
2x2 Float64 Array: |
0.153756 0.368514 |
1.15119 0.918912 |
julia> @spawn rand2(2,2) |
RemoteRef(1,1,1) |
julia> @spawn rand2(2,2) |
RemoteRef(2,1,2) |
julia> exception on 2: in anonymous: rand2 not defined |
rand2
函数,但进程 2 不知道。 require
函数自动在当前所有可用的处理器上载入源文件,使所有的处理器都能运行代码:
julia> require("myfile") |
@everywhere
宏在所有处理器上执行命令:
julia> @everywhere id = myid() |
julia> remotecall_fetch(2, ()->id) |
2 |
@everywhere include("defs.jl") |
julia -p <n> -L file1.jl -L file2.jl driver.jl |
workers
。当只有一个进程的时候,进程 1 就被当做一个 worker。否则,worker 就是指除了进程 1 之外的所有进程。
Julia 内置有对于两种集群的支持:
—p
选项。
--machinefile
选项。它使用一个无密码的 ssh
登来在指定的机器上启动 julia 工作进程(以相同的路径作为当前主机)。
函数 addprocs
,rmprocs
,workers
,当然还有其他的在一个集群中可用的以可编程的方式进行添加,删除和查询的函数。
其他类型的集群可以通过编写自己的自定义 ClusterManager。请参阅 ClusterManagers 部分。
并行计算中,消息传递和数据移动是最大的开销。减少这两者的数量,对性能至关重要。
fetch
是显式的数据移动操作,它直接要求将对象移动到当前机器。 @spawn
(及相关宏)也进行数据移动,但不是显式的,因而被称为隐式数据移动操作。对比如下两种构造随机矩阵并计算其平方的方法
# method 1 |
A = rand(1000,1000) |
Bref = @spawn A^2 |
... |
fetch(Bref) |
# method 2 |
Bref = @spawn rand(1000,1000)^2 |
... |
fetch(Bref) |
方法 1 中,本地构造了一个随机矩阵,然后将其传递给做平方计算的处理器。方法 2 中,在同一处理器构造随机矩阵并进行平方计算。因此,方法 2 比方法 1 移动的数据少得多。
大部分并行计算不需要移动数据。最常见的是蒙特卡罗仿真。下例使用 @spawn
在两个处理器上仿真投硬币。先在 count_heads.jl
中写如下函数:
function count_heads(n) |
c::Int = 0 |
for i=1:n |
c += randbool() |
end |
c |
end |
require("count_heads") |
a = @spawn count_heads(100000000) |
b = @spawn count_heads(100000000) |
fetch(a)+fetch(b) |
上例中,我们显式调用了两个 @spawn
语句,它将并行计算限制在两个处理器上。要在任意个数的处理器上运行,应使用 并行 for 循环 ,它在 Julia 中应写为:
nheads = @parallel (+) for i=1:200000000 |
int(randbool()) |
end |
这个构造实现了给多处理器分配迭代的模式,并且使用特定约简来综合结果(此例中为 (+)
)。
注意,尽管并行 for 循环看起来和一组 for 循环差不多,但它们的行为有很大区别。第一,循环不是按顺序进行的。第二,写进变量或数组的值不是全局可见的,因为迭代运行在不同的处理器上。并行循环内使用的所有变量都会被复制、广播到每个处理器。
下列代码并不会按照预想运行:
a = zeros(100000) |
@parallel for i=1:100000 |
a[i] = i |
end |
a
的所有元素,因为每个处理器上都只有独立的一份儿。应避免类似的并行 for 循环。但是我们可以使用分布式数组来规避这种情形,后面我们会讲。
如果“外部”变量是只读的,就可以在并行循环中使用它:
a = randn(1000) |
@parallel (+) for i=1:100000 |
f(a[randi(end)]) |
end |
pmap
函数。下例中并行计算几个大随机矩阵的奇异值:
M = {rand(1000,1000) for i=1:10} |
pmap(svd, M) |
pmap
,反之,则使用 @parallel for
。
Julia 的平行编程平台使用任务(也成为协程) ,其可在多个计算中切换。每当代码执行一个通信操作,例如 fetch
或者 wait
,当前任务便暂停同时调度器会选择另一个任务运行。在事件等待完成后,任务会重新启动。
对于很多问题,没必要直接考虑任务。然而,由于提供了动态调度,可以同时等待多个事件。在动态调度中,一个程序决定计算什么和在哪计算,这是基于其他工作何时完成的。这是被不可预知的或不可平衡的工作荷载所需要的,只有当他们结束当前任务我们才能分配更多的工作进程。
作为一个例子,考虑计算不同大小的矩阵的奇异值:
M = {rand(800,800), rand(600,600), rand(800,800), rand(600,600)} |
pmap(svd, M) |
pmap
的实现过程中可以看到这个:
function pmap(f, lst) |
np = nprocs() # determine the number of processes available |
n = length(lst) |
results = cell(n) |
i = 1 |
# function to produce the next work item from the queue. |
# in this case it's just an index. |
nextidx() = (idx=i; i+=1; idx) |
@sync begin |
for p=1:np |
if p != myid() || np == 1 |
@async begin |
while true |
idx = nextidx() |
if idx > n |
break |
end |
results[idx] = remotecall_fetch(p, f, lst[idx]) |
end |
end |
end |
end |
end |
results |
end |
只有在本地运行任务的过程中,@async
才与 @spawn
类似。我们使用它来为每个流程创建一个“供给”的任务。每个任务选择下一个需要被计算的指数,然后等待它的进程完成,接着一直重复到用完指数。
注意,“供给”任务只有当主要任务到达 @sync
块结束时才开始执行,此时它放弃控制并等待所有的本地任务在从函数返回之前完成。供给任务可以通过 nextidx()
共享状态,因为它们都在相同的进程上运行。
这个过程不需要锁定,因为线程是实时进行调度的而不是一成不变。这意味着内容的切换只发生在定义好的时候:在这种情况下,当 remotecall_fetch
会被调用。
并行计算综合使用多个机器上的内存资源,因而可以使用在一个机器上不能实现的大数组。这时,可使用分布式数组,每个处理器仅对它所拥有的那部分数组进行操作。
分布式数组(或 全局对象 )逻辑上是个单数组,但它分为很多块儿,每个处理器上保存一块儿。但对整个数组的运算与在本地数组的运算是一样的,并行计算是隐藏的。
分布式数组是用 DArray
类型来实现的。 DArray
的元素类型和维度与 Array
一样。 DArray
的数据的分布,是这样实现的:它把索引空间在每个维度都分成一些小块。
一些常用分布式数组可以使用 d
开头的函数来构造:
dzeros(100,100,10) |
dones(100,100,10) |
drand(100,100,10) |
drandn(100,100,10) |
dfill(x, 100,100,10) |
x
来初始化。这些函数自动选取某个分布。如果要指明使用哪个进程,如何分布数据,应这样写:
dzeros((100,100), [1:4], [1,4]) |
DArrays
放在一个进程的子集中,该进程允许多个 DArray
同时计算,并且每个进程拥有更高比例的通信工作。
第三个参数指定了一个分布;数组第 n 个元素指定了应该分成多少个块。在本例中,第一个维度不会分割,而第二个维度将分为四块。因此每个局部块的大小为 (100,25)
。注意,分布式数组必须与进程数量相符。
distribute(a::Array)
可用来将本地数组转换为分布式数组。
localpart(a::DArray)
可用来获取 DArray
本地存储的部分。
localindexes(a::DArray)
返回本地进程所存储的维度索引值范围多元组。
convert(Array, a::DArray)
将所有数据综合到本地进程上。
使用索引值范围来索引 DArray
(方括号)时,会创建 SubArray
对象,但不复制数据。
DArray
的构造函数是 darray
,它的声明如下:
DArray(init, dims[, procs, dist]) |
init
函数的参数,是索引值范围多元组。这个函数在本地声名一块分布式数组,并用指定索引值来进行初始化。 dims
是整个分布式数组的维度。 procs
是可选的,指明一个存有要使用的进程 ID 的向量 。 dist
是一个整数向量,指明分布式数组在每个维度应该被分成几块。
最后俩参数是可选的,忽略的时候使用默认值。
下例演示如果将本地数组 fill
的构造函数更改为分布式数组的构造函数:
dfill(v, args...) = DArray(I->fill(v, map(length,I)), args...) |
此例中 init
函数仅对它构造的本地块的维度调用 fill
。
在这个时候,分布式数组没有太多的功能。主要功能是通过数组索引来允许进行通信,这对许多问题来说都很方便。
作为一个例子,考虑实现“生活”细胞自动机,每个单元网格中的细胞根据其邻近的细胞进行更新。每个进程需要其本地块中直接相邻的细胞才能计算一个迭代的结果。下面的代码可以实现这个功能:
function life_step(d::DArray) |
DArray(size(d),procs(d)) do I |
top = mod(first(I[1])-2,size(d,1))+1 |
bot = mod( last(I[1]) ,size(d,1))+1 |
left = mod(first(I[2])-2,size(d,2))+1 |
right = mod( last(I[2]) ,size(d,2))+1 |
old = Array(Bool, length(I[1])+2, length(I[2])+2) |
old[1 , 1 ] = d[top , left] # left side |
old[2:end-1, 1 ] = d[I[1], left] |
old[end , 1 ] = d[bot , left] |
old[1 , 2:end-1] = d[top , I[2]] |
old[2:end-1, 2:end-1] = d[I[1], I[2]] # middle |
old[end , 2:end-1] = d[bot , I[2]] |
old[1 , end ] = d[top , right] # right side |
old[2:end-1, end ] = d[I[1], right] |
old[end , end ] = d[bot , right] |
life_rule(old) |
end |
end |
old
。注意,do
块语法方便 init
函数传递给 DArray
构造函数。接下来,连续函数 life_rule
被调用以提供数据的更新规则,产生所需的 DArray
块。 life_rule
与 DArray-specific
没有关系,但为了完整性,我们在此仍将它列出:
function life_rule(old) |
m, n = size(old) |
new = similar(old, m-2, n-2) |
for j = 2:n-1 |
for i = 2:m-1 |
nc = +(old[i-1,j-1], old[i-1,j], old[i-1,j+1], |
old[i ,j-1], old[i ,j+1], |
old[i+1,j-1], old[i+1,j], old[i+1,j+1]) |
new[i-1,j-1] = (nc == 3 || nc == 2 && old[i,j]) |
end |
end |
new |
end |
共享阵列使用在许多进程中共享内存来映射相同数组的系统。虽然与 DArray
有一些相似之处,但是 SharedArray
的行为是完全不同的。在一个 DArray
中,每个进程只能本地访问一块数据,并且两个进程共享同一块;相比之下,在 SharedArray
中,每个“参与”的进程能够访问整个数组。
当你想要在同一台机器上大量数据共同访问两个或两个以上的进程时, SharedArray
是一个不错的选择。
SharedArray
索引(分配和访问值)与常规数组一样工作,并且是非常高效的,因为其底层内存可用于本地进程。因此,大多数算法自然地在 SharedArrays
上运行,即使在单进程模式中。当某个算法必须在一个 Array
输入的情况下,可以从 SharedArray
检索底层数组通过调用 sdata(S)
取回。对于其他 AbstractArray
类型, sdata
返回对象本身,所以在任何数组类型下使用 sdata
都是很安全的。
共享数字构造函数数的形式:
SharedArray(T::Type, dims::NTuple; init=false, pids=Int[]) |
创建一个被 pids
进程指定的,bitstype 为 T
并且大小为 dims
的共享数组。与分布式阵列不同,共享数组只能用于这些参与人员指定的以 pid
命名的参数(如果在同一个主机上,创建过程也同样如此)。
如果一个签名为 initfn(S::SharedArray)
的 init
函数被指定,它会被所有参与人员调用。你可以控制它,每个工人可以在数组的不同部分运行 init
函数,因此进行并行的初始化。
这里有一个简单的例子:
julia> addprocs(3) |
3-element Array{Any,1}: |
2 |
3 |
4 |
julia> S = SharedArray(Int, (3,4), init = S -> S[localindexes(S)] = myid()) |
3x4 SharedArray{Int64,2}: |
2 2 3 4 |
2 3 3 4 |
2 3 4 4 |
julia> S[3,2] = 7 |
7 |
julia> S |
3x4 SharedArray{Int64,2}: |
2 2 3 4 |
2 3 3 4 |
2 7 4 4 |
localindexes
提供不相交的一维索引的范围,它有时方便进程之间的任务交流。当然,你可以按你希望的方式来划分工作:
julia> S = SharedArray(Int, (3,4), init = S -> S[myid()-1:nworkers():length(S)] = myid()) |
3x4 SharedArray{Int64,2}: |
2 2 2 2 |
3 3 3 3 |
4 4 4 4 |
@sync begin |
for p in workers() |
@async begin |
remotecall_wait(p, fill!, S, p) |
end |
end |
end |
pid
来充满整个数组,无论最后执行的是哪一个进程(任何特定元素 S
)都将保留他的 pid
。