6

任何人都可以帮忙吗？我是一个相当有经验的 Matlab 用户，但是在加速下面的代码时遇到了麻烦。

使用 12 个核心，我能够在所有三个循环中一次运行的最快时间约为 200 秒。实际函数将被调用约 720 次，按此速率执行将需要 40 多个小时。根据 Matlab 分析器，大部分 CPU 时间都花在了指数函数调用上。我已经设法使用 gpuArray 大大加快了速度，然后在 Quadro 4000 显卡上运行 exp 调用，但这会阻止使用 parfor 循环，因为工作站只有一个显卡，这会抹杀任何收益。任何人都可以提供帮助，或者这段代码是否接近使用 Matlab 可以实现的最佳值？我用 openMP 编写了一个非常粗略的 c++ 实现，但收效甚微。

提前谢谢了

function SPEEDtest_CPU

% Variable setup:
% - For testing I'll use random variables. These will actually be fed into 
%   the function for the real version of this code.
sy    = 320;
sx    = 100;
sz    = 32;
A     = complex(rand(sy,sx,sz),rand(sy,sx,sz));
B     = complex(rand(sy,sx,sz),rand(sy,sx,sz));
C     = rand(sy,sx);
D     = rand(sy*sx,1);
F     = zeros(sy,sx,sz);
x     = rand(sy*sx,1);  
y     = rand(sy*sx,1);
x_ind = (1:sx) - (sx / 2) - 1;
y_ind = (1:sy) - (sy / 2) - 1;


% MAIN LOOPS 
%  - In the real code this set of three loops will be called ~720 times!
%  - Using 12 cores, the fastest I have managed is ~200 seconds for one
%    call of this function.
tic
for z = 1 : sz
    A_slice = A(:,:,z);
    A_slice = A_slice(:);
    parfor cx = 1 : sx       
        for cy = 1 : sy       
            E = ( x .* x_ind(cx) ) + ( y .* y_ind(cy) ) + ( C(cy,cx) .* D );                                                          

            F(cy,cx,z) = (B(cy,cx,z) .* exp(-1i .* E))' * A_slice; 
        end       
    end   
end
toc

end
4

8 回答 8

3

需要考虑的一些事情：

你有没有考虑过单身？

您能否对 cx、cy 部分进行矢量化，以便它们表示数组操作？

考虑更改浮点舍入或信号模式。

于 2013-10-04T10:13:30.070 回答
2

如果您的数据是真实的（不复杂），如您的示例所示，您可以节省替换时间

(B(cy,cx,z) .* exp(-1i .* E))'

经过

(B(cy,cx,z) .* (cos(E)+1i*sin(E))).'

具体来说，在我的机器上(cos(x)+1i*sin(x)).'花费的时间比exp(-1i .* x)'.

If Aand Bare complex:E仍然是真实的，因此您可以Bconj = conj(B)在循环之外预先计算（这需要大约 10 毫秒的数据大小，并且只完成一次）然后替换

(B(cy,cx,z) .* exp(-1i .* E))'

经过

(Bconj(cy,cx,z) .* (cos(E)+1i*sin(E))).'

以获得类似的收益。

于 2013-10-04T10:28:35.687 回答
1

加速 MATLAB 代码的主要方法有两种；预分配和向量化。

您已经很好地进行了预分配，但没有矢量化。为了最好地学习如何做到这一点，您需要很好地掌握线性代数以及repmat如何将向量扩展为多个维度。

矢量化可以导致多个数量级的加速，并将最佳地使用内核（假设标志向上）。

您正在计算的数学表达式是什么，我可以帮忙吗？

于 2013-10-04T10:21:56.790 回答
1

您可以移出x .* x_ind(cx)最里面的循环。我没有方便的 GPU 来测试时序，但您可以将代码分成三个部分，以允许您使用 GPU 和 parfor

for z = 1 : sz
    E = zeros(sy*sx,sx,sy);
    A_slice = A(:,:,z);
    A_slice = A_slice(:);
    parfor cx = 1 : sx
        temp = ( x .* x_ind(cx) );       
        for cy = 1 : sy       
            E(:, cx, cy) = temp + ( y .* y_ind(cy) ) + ( C(cy,cx) .* D );                                                          
        end
    end
    temp = zeros(zeros(sy*sx,sx,sy));
    for cx = 1 : sx
        for cy = 1 : sy       
             % Ideally use your GPU magic here
             temp(:, cx, cy) = exp(-1i .* E(:, cx, cy)));
        end
    end
    parfor cx = 1 : sx
        for cy = 1 : sy       
            F(cy,cx,z) = (B(cy,cx,z) .* temp(:, cx, cy)' * A_slice; 
        end       
    end   
end
于 2013-10-04T11:16:29.657 回答
0

为了允许适当的并行化，您需要确保循环是完全独立的，因此检查E在每次运行中不分配是否有帮助。

此外，尽量矢量化，一个简单的例子可能是：y.*y_ind(cy)

如果您只是一次为所有值创建适当的索引，则可以将其从最低循环中取出。

于 2013-10-04T10:22:40.463 回答
0

除了其他人在这里给出的其他好建议之外，乘法与循环A_slice无关，cx,cy可以在循环之外进行，F一旦两个循环都完成就相乘。

类似地， 的共轭B*exp(...)也可以在cx,cy循环之外进行，在乘以之前A_slice。

于 2013-10-10T08:14:48.427 回答
0

不确定它是否对速度有很大帮助 - 但由于 E 基本上是一个总和，也许你可以使用它exp (i cx(A+1)x) = exp(i cx(A) x) * exp(i x)并且exp(i x)可以预先计算。

这样你就不必在每次迭代时评估 exp ——而只需要乘法，这应该会更快。

于 2013-10-04T10:57:27.300 回答
0

这一行： ( x .* x_ind(cx) ) + ( y .* y_ind(cy) ) + ( C(cy,cx) .* D );

是某种类型的卷积，不是吗？循环卷积在频域中要快得多，并且使用 FTT 优化到/从频域的转换。

于 2017-02-04T17:03:22.587 回答