【问题标题】：Matlab CUDA basic experimentMatlab CUDA基础实验
【发布时间】：2013-04-06 00:03:40
【问题描述】：

（正确且有指导意义的回答，见下文）

我开始用 matlab 和 gpu (nvidia gtx660) 做实验。

现在，我编写了这个简单的蒙特卡罗算法来计算 PI。以下是CPU版本：

function pig = mc1vecnocuda(n)
countr=0;
A=rand(n,2);
 for i=1:n

   if norm(A(i,:))<1
    countr=countr+1;
   end
 end
pig=(countr/n)*4;
end

这需要很少的时间在 CPU 上执行，将 100000 个点“扔”到单位圆中：

   >> tic; mc1vecnocuda(100000);toc;

      Elapsed time is 0.092473 seconds.

相反，看看 gpu 化版本的算法会发生什么：

   function pig = mc1veccuda(n)
   countr=0;
   gpucountr=gpuArray(countr);
   A=gpuArray.rand(n,2);
   parfor (i=1:n,1024)
    if norm(A(i,:))<1
        gpucountr=gpucountr+1;
    end
   end

   pig=(gpucountr/n)*4;
   end

现在，这需要很长时间才能执行：

>> tic; mc1veccuda(100000);toc;
Elapsed time is 21.137954 seconds.

我不明白为什么。我对 1024 个工作人员使用 parfor 循环，因为使用 gpuDevice 查询我的 nvidia 卡，1024 是 gtx660 上允许的最大并发线程数。

有人可以帮助我吗？谢谢。

编辑：这是避免 IF: 的更新版本：

function pig = mc2veccuda(n)
countr=0;
gpucountr=gpuArray(countr);
A=gpuArray.rand(n,2);
parfor (i=1:n,1024)

    gpucountr = gpucountr+nnz(norm(A(i,:))<1);

end

pig=(gpucountr/n)*4;
end

这是按照 Bichoy 的指导方针编写的代码（正确代码实现结果）：

function pig = mc3veccuda(n)
countr=0;
gpucountr=gpuArray(countr);
A=gpuArray.rand(n,2);
Asq = A.^2;
Asqsum_big_column = Asq(:,1)+Asq(:,2);
Anorms=Asqsum_big_column.^(1/2);
gpucountr=gpucountr+nnz(Anorms<1);

pig=(gpucountr/n)*4;
end

请注意 n=1000 万的执行时间：

>> tic; mc3veccuda(10000000); toc;
Elapsed time is 0.131348 seconds.
>> tic; mc1vecnocuda(10000000); toc;
Elapsed time is 8.108907 seconds.

我没有测试我的原始 cuda 版本 (for/parfor)，因为它的执行需要几个小时，n=10000000。

伟大的比丘！ ;)

【问题讨论】：

我很确定在 for 循环的执行过程中，某些东西会从设备移动到主机并再次返回。有任何方法可以观察是否（以及何时）传输设备主机出现？？
这可能不是最好的方法。整个事情可以向量化为：pig = nnz(sum(X.^2,2)<=1) * 4 / size(X,1) where X=rand(N,2) or similar array on the gpu
跑题了，但这里有一个很好的模拟动画：pastebin.com/w0N46vJE :) 类似于：en.wikipedia.org/wiki/File:Pi_30K.gif
欢迎您的链接！

标签： matlab cuda

【解决方案1】：

我猜问题出在parfor！

parfor 应该在 MATLAB worker 上运行，那是你的主机而不是 GPU！我猜实际发生的情况是您在主机上（而不是在 GPU 上）启动了 1024 个线程，并且每个线程都在尝试调用 GPU。这会导致您的代码花费大量时间。

尝试重新编写代码以使用矩阵和数组运算，而不是 for 循环！这将显示一些加速。另外，请记住，您应该在 GPU 中进行更多计算，否则，内存传输只会支配您的代码。

代码：

这是包含几个人的所有更正和建议后的最终代码：

function pig = mc2veccuda(n)
  A=gpuArray.rand(n,2); % An nx2 random matrix
  Asq = A.^2; % Get the square value of each element
  Anormsq = Asq(:,1)+Asq(:,2); % Get the norm squared of each point
  gpucountr = nnz(Anorm<1); % Check the number of elements < 1
  pig=(gpucountr/n)*4;

【讨论】：

你太棒了，给我上了一课。我想您出于教育目的而使代码有些混乱，请参阅我的原始帖子编辑以获取更多详细信息！非常感谢，通过您的回答，我开始学习如何以正确的方式矢量化代码；）
我很高兴它帮助了你。这就是这个表格的全部目的。抱歉我的代码缺少教育性的 cmets，我会尝试对其进行编辑以使其更清晰，我很匆忙地回复了。谢谢你的好话。
顺便说一句，如果您认为合适，请随时建议对我的帖子进行任何编辑以使其更清晰。如果您无法进行编辑，只需在此处将它们添加到 cmets 中，我会对其进行调整。
我更正了您在最后一个问题中发布的答案。再次感谢:)
如您所见，其他人无法建议正确的路径。这并不是因为他们是糟糕的程序员，而是因为我们都受过从传统循环的角度进行思考的教育。正如您正确建议的那样，并行编程是另一回事，必须重新塑造他的思维方式；）

【解决方案2】：

很多原因，例如：

主机和设备之间的数据移动
每个循环内的计算量非常小
在 GPU 上调用 rand 可能不是并行的
if 循环内的条件可能会导致分歧
对公共变量的累加可能会以串行方式运行，但会产生开销

很难分析 Matlab+CUDA 代码。您可能应该尝试使用原生 C++/CUDA 并使用并行 Nsight 来查找瓶颈。

【讨论】：

感谢您的回复。根据这篇文章：stackoverflow.com/questions/9054949/… 你是对的。现在，我使用 nnz 函数而不是 IF 修改了我的代码，并且 r2013a 在 gpuArray 上支持 nnz。我得到了同样可怕的结果，所以瓶颈必须在 FOR 循环中。我很难理解的原因，似乎每次迭代的数据都会来回移动，这与我的信念相反。
你能用新号码更新你的帖子吗？我很想知道。
循环体非常小，你在同一个变量中累积（很难并行化）。检查是否有任何内置的归约操作。
当然：>> tic; mc2veccuda(100000);目录；经过的时间是 19.631224 秒。
注意我正在累积一个驻留在 gpu 内存中的变量，但我知道并行化是多么困难。你是不是建议分块积累，然后总结？？请更具体，如果可能，请提供代码 sn-ps.. 非常感谢！ ;)

【解决方案3】：

正如 Bichoy 所说，CUDA 代码应始终进行矢量化处理。在 MATLAB 中，除非您正在编写 CUDA 内核，否则您获得的唯一大幅加速是在具有数千个（慢速）内核的 GPU 上调用矢量化操作。如果您没有大向量和向量化代码，那将无济于事。

另一件没有提到的事情是，对于像 GPU 这样的高度并行架构，您希望使用与“标准”算法不同的随机数生成算法。因此，添加到 Bichoy 的答案中，添加参数“Threefry4x64”（64 位）或“Philox4x32-10”（32 位，速度更快！超快！）可以导致 CUDA 代码的大幅加速。 MATLAB 在这里解释了这一点：http://www.mathworks.com/help/distcomp/examples/generating-random-numbers-on-a-gpu.html

【讨论】：