Как мне индексировать codistributed массивы в блоке 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и документация для 2013a пятнистый на том, как повернуть Composite данные в матрицу (cell2mat не работает). Это также плохо масштабируется, потому что чем больше у меня ядер, тем больше матрица, и эту большую матрицу приходится копировать каждому работнику, который затем игнорирует большую часть данных. Вопрос 1: как превратить составной тип данных в пригодный для использования массив?

Второе, что я попробовал, было использовать codistributed массив.

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

Это говорит мне о том, что у каждого работника есть один вектор размера [K 1], что, как я считаю, является тем, что я хочу, но когда я пытаюсь объединить вышеупомянутые методы

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: как правильно индексировать в и из массива codistributed в spmd блок?

Answer 1

Что касается вопроса № 1, то Composite переменная в клиенте в основном относится к нераспределенному массиву вариантов, хранящемуся на рабочих. Вы можете получить доступ к массиву от каждого работника: {}-индексация с использованием соответствующего лабиндекса (например: 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