如何在spmd块中索引共分布数组

Question

我正在做一个非常大的计算(大气吸收)，它有许多个别的窄峰，它们都是在最后加起来的。对于每个峰值，我都预先计算了峰值形状函数的值超过我所选择的阈值的范围，然后我将逐行进行，并将这些峰值添加到我的光谱中。下面给出了一个最低限度的例子：

X = 1:1e7;
K = numel(a); % count the number of peaks I have.
spectrum = zeros(size(X));
for k = 1:K
    grid = X >= rng(1,k) & X <= rng(2,k);
    spectrum(grid) = spectrum(grid) + peakfn(X(grid),a(k),b(k),c(k)]);
end

在这里，每个峰值都有一些参数来定义位置和形状(a、b、c)，以及计算范围(rng)。这是很好的工作，在我的机器上它的基准测试大约220秒来完成一个完整的数据集。然而，我有一台4核心机器，我最终想要在集群上运行它，所以我想并行化它并使它可扩展。

因为每个循环都依赖于上一次迭代的结果，所以我不能使用parfor，所以我开始学习如何使用spmd块。我的第一次尝试是这样的：

X = 1:1e7;
cores = matlabpool('size');
K = numel(a);
spectrum = zeros(size(X),cores);
spmd
    n = labindex:cores:K
    N = numel(n);
    for k = 1:N
        grid = X >= rng(1,n(k)) & X <= rng(2,n(k));
        spectrum(grid,labindex) = spectrum(grid,labindex) + peakfn(X(grid),a(n(k)),b(n(k)),c(n(k))]);
    end
end
finalSpectrum = sum(spectrum,2);

这几乎成功了。程序在最后一行崩溃，因为spectrum是Composite类型，关于如何将Composite数据转换为矩阵(cell2mat不起作用)，2013a的文档不清楚。这也不能很好地扩展，因为我拥有的核越多，矩阵就越大，而这个大型矩阵必须被复制到每个工作人员，然后忽略大部分数据。问题1:如何将复合数据类型转换为可使用的数组？

我尝试的第二件事是使用codistributed数组。

spmd
    spectrum = codistributed.zeros(K,cores);
    disp(size(getLocalPart(spectrum)))
end

这告诉我，每个工人都有一个K1大小的向量，我相信这就是我想要的，但是当我尝试将上述方法结合起来时，

spmd
    spectrum = codistributed.zeros(K,cores);
    n = labindex:cores:K
    N = numel(n);
    for k = 1:N
        grid = X >= rng(1,n(k)) & X <= rng(2,n(k));
        spectrum(grid) = spectrum(grid) + peakfn(X(grid),a(n(k)),b(n(k)),c(n(k))]);        end
    finalSpectrum = gather(spectrum);
end
finalSpectrum = sum(finalSpectrum,2);

我得到了Matrix dimensions must agree错误。由于它是在一个并行块中，所以我不能使用我的常规调试拐杖，通过遍历循环，查看每个块在每个点的大小来查看正在发生的事情。问题2:，什么是在spmd块中对一个共分布数组进行索引的适当方法？

Answer 1

对于question#1，客户机中的Composite变量基本上是指存储在工作人员上的非分布式变体数组。您可以通过{}-indexing使用相应的labindex (例如：spectrum{1}、spectrum{2}、.)从每个工作人员访问数组。

对于您的代码是：finalSpectrum = sum(cat(2,spectrum{:}), 2);

现在我用随机数据自己尝试了这个问题。下面是三个需要比较的实现(请参阅这里了解分布式和非分布式阵列之间的区别)。首先，我们从共同的数据开始：

len = 100;    % spectrum length
K = 10;       % number of peaks
X = 1:len;

% random position and shape parameters
a = rand(1,K); b = rand(1,K); c = rand(1,K);

% random peak ranges (lower/upper thresholds)
ranges = sort(randi([1 len], [2 K]));

% dummy peakfn() function
fcn = @(x,a,b,c) x+a+b+c;

% prepare a pool of MATLAB workers
matlabpool open

( 1)串行循环：

spectrum = zeros(size(X));
for i=1:size(ranges,2)
    r = ranges(:,i);
    idx = (r(1) <= X & X <= r(2));
    spectrum(idx) = spectrum(idx) + fcn(X(idx), a(i), b(i), c(i));
end
s1 = spectrum;

clear spectrum i r idx

2)具有组合阵列的SPMD

spmd
    spectrum = zeros(1,len);
    ind = labindex:numlabs:K;
    for i=1:numel(ind)
        r = ranges(:,ind(i));
        idx = (r(1) <= X & X <= r(2));
        spectrum(idx) = spectrum(idx) + ...
            feval(fcn, X(idx), a(ind(i)), b(ind(i)), c(ind(i)));
    end
end
s2 = sum(vertcat(spectrum{:}));

clear spectrum i r idx ind

3)共分布阵列的SPMD。

spmd
    spectrum = zeros(numlabs, len, codistributor('1d',1));
    ind = labindex:numlabs:K;
    for i=1:numel(ind)
        r = ranges(:,ind(i));
        idx = (r(1) <= X & X <= r(2));
        spectrum(labindex,idx) = spectrum(labindex,idx) + ...
            feval(fcn, X(idx), a(ind(i)), b(ind(i)), c(ind(i)));
    end
end
s3 = sum(gather(spectrum));

clear spectrum i r idx ind

这三个结果都应该相等(在可接受的小误差范围内)。

>> max([max(s1-s2), max(s1-s3), max(s2-s3)])
ans =
   2.8422e-14