multithreading - 如何获得多线程dot()函数？

Question

在具有8个逻辑CPU内核的64位Windows 8.1 PC上，在Julia v0.3.7的REPL中执行此命令时:

blas_set_num_threads(CPU_CORES)
const v=ones(Float64,100000)
@time for k=1:1000000;s=dot(v,v);end

我在任务管理器或Process Explorer的CPU表中观察到仅使用了12.5％的CPU(1个逻辑CPU内核)。在Windows 7和Windows 8.1上，Julia v0.3.5也观察到相同的情况。通过在命令行上启动“Julia -p 8”，我也观察到了相同的行为。返回运行不带“-p 8”命令行选项的Julia REPL，我尝试了以下测试:

blas_set_num_threads(CPU_CORES)
@time peakflops(10000)

在这种情况下，CPU仪表显示100％的CPU使用率。

由于dot()和peakflops()都使用BLAS(在我的情况下为OpenBLAS)，因此我希望它们都使用blas_set_num_threads()指定的线程数。但是，仅后者功能实际上起作用。 dot()的行为是否是由于错误(可能是OpenBLAS中的错误)引起的？

我试图通过使用矩阵乘法功能来解决Julia的缺点。但是，我正在对GB大小的2D数组的子向量上的dot()操作进行迭代，其中子向量使用连续的内存。矩阵乘法迫使我转置每个向量，从而创建了一个副本。在内循环中这是昂贵的。因此，对我来说，选择似乎是学习如何使用Julia的并行处理命令/宏，还是回到Python(英特尔的MKL BLAS如ddot()所预期的那样运行)。由于dot()是在我尝试编写的例程中消耗99％CPU的函数，因此我希望OpenBLAS在Julia中为我提供用户友好的最佳解决方案，而不必担心引擎盖下所有并行处理的复杂性。但是也许还不错...

我可以使用一些创建多线程的dot()函数。一些示例代码将是最好的帮助。当所有线程都在同一台计算机上运行时，我需要使用SharedArray吗？如果是这样，从Array转换为SharedArray是否会创建一个副本？由于我的数组很大，因此我不想同时在内存中复制两个副本。因为我的向量长度约为100,000，并且来自一个数组的向量以不可预测的顺序使用，所以对我而言，最好的多线程解决方案是dot()函数，该函数将任务拆分为可用的核心，并对结果求和每个核心。如何在Julia中像BLAS一样有效地做到这一点？

我在这里看到了蒂姆·霍利(Tim Holy)的回答:
BLAS v. parallel updates for Julia SharedArray objects
但是，他的示例(我认为)在一个内核上执行了整个dot()函数，并且无法回答我的其他问题。

编辑1:
我尝试使用“-p 8”命令行选项运行Julia，并从此处将上面示例中的dot()替换为innersimd():
http://docs.julialang.org/en/release-0.3/manual/performance-tips/
仍仅使用1个核心。我修改了innersimd()，将其参数键入为::Array{Float64, 1}，然后键入::SharedArray{Float64, 1}，但这仍然使用1个核心。 :(

编辑2:
我测试了Julia的矩阵乘法(BLAS的gemm!()函数):

blas_set_num_threads(CPU_CORES)
const A=ones(Float64,(4,100000))
const B=ones(Float64,(100000,4))
@time for k=1:100000;s=A*B;end

即使Julia开始时没有使用“-p”命令行选项，Julia也仅使用一个内核。

编辑3:
这是Python的类似测试代码:

import numpy as np
from scipy.linalg.blas import ddot
from timeit import default_timer as timer
v = np.ones(100000)
start = timer()
for k in range(1000000):
    s = ddot(v,v)
exec_time=(timer() - start)
print
print("Execution took", str(round(exec_time, 3)), "seconds")

在64位Anaconda3 v2.1.0和WinPython上，我得到了相同的结果:7.5秒。与我的问题中的第一个示例进行比较，将Julia 0.3.7与OpenBLAS一起使用在28秒内执行。这使Python的速度比Julia快4倍，这可以通过OpenBLAS的ddot()的单线程实现来解释。

编辑4:
我在Python中进行了一些测试，以有效地在内部循环(N = 100000)中执行(4xN)*(Nx2)矩阵乘法，结果发现它的速度要慢两倍以上。我怀疑这是因为缓存需求现在增加了8倍，所以缓存命中率更差。为了保持Julia在缓存中的需求与Python相同，Julia的问题是:如何将我的100000个长度的向量分成4个块，并并行执行OpenBLAS的单线程ddot()(并对结果求和)？该长度可能不完全是4的倍数。由于OpenBLAS是单线程的，因此我可以在命令行中使用“-p 8”命令行参数运行Julia，并在同一 session 中对其他自定义函数进行多线程处理。

编辑5:
使用“-p 8”命令行参数运行Julia v0.3.7，我观察到OpenBLAS gemm!()函数确实以多线程运行(对于某些矩阵尺寸):

blas_set_num_threads(CPU_CORES)
const a = rand(10000, 10000)
@time a * a

Answer 1

原因是OpenBLAS的ddot似乎不是多线程的。我认为它们主要为BLAS-3函数(例如xgemm)实现多线程，因此，如果可能的话，最好的办法是使用矩阵乘法编写问题。

如果您可以举一个如何尝试矩阵乘法的例子，那将会有所帮助。也许可以避免换位。您是如何在Python中做到这一点的？一切都只是BLAS，所以Julia在这里应该能和Python一样好。

使用-p 8标志运行Julia会启动Julia的九个副本，并关闭BLAS中的多线程，因此这可能会使您的情况更糟。

编辑1:
我在单线程和多线程模式下使用OpenBLAS和MKL尝试了您的示例，如果有任何区别，单线程会更快。我不确定OpenBLAS是否为您的示例使用了多个线程。对于细小的问题，可能很难实现多线程加速，因此，对于较严重的问题，它们只能打开多线程。如果您的示例在同一台计算机上运行时，使用Python的速度明显更快？如果是这样，请提供Python示例，以便可以与我的计算机上的示例进行比较。

编辑2:
我可以确认MKL中ddot的加速。确切地说，Edit 1中的注释是关于dgemm的。正在使Julia成为多线程的积极工作，当这种情况发生时，Julia的线程实现可能比OpenBLAS的ddot更快。您现在的选择是以下之一

如果可能，请重写您的问题，以使OpenBLAS的dgemm使用多个线程，即使用较胖的矩阵。在我的机器上，n=16似乎加入了多线程。如果可能的话，我认为这不会比MKL慢。通常，OpenBLAS的dgemm与MKL的比较好。

要求OpenBLAS开发人员使ddot多线程。如果您可以仅向项目提供一个版本，那会更好，但是高性能的BLAS代码很难编写。

购买MKL并开始销售Revolution Julia(包括MKL)以支付许可证费用。

编辑3:比较都是在具有Nahalem体系结构和80个内核的Ubuntu计算机上的Julia中进行的(但我最多使用16个)。我使用了Julia相同开发版本的两种不同的构建:一种与从源代码构建的OpenBLAS链接，另一种与MKL的链接。我还在最新的Haswell系统上使用OpenBLAS尝试了Julia，以确保未在该体系结构上为ddot实现线程。