Matlab Parfor Scheibe korrekt

Question

Ich habe zwei verschachtelte Schleifen, die ich parallelisieren möchte.

n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = zeros(1,m); 
for i=1:n 
    q = ones(1,m); 
    for j=1:n 
     q = q .* (xx-x(j))/(x(i)-x(j)); 
    end 
    r = r + q; 
end 

Um diese Funktion für die Palatalisierung vorzubereiten, habe ich lokale Variablen in globale geändert.

n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = ones(n,m); 
for i=1:n 
    for j=1:n 
     r(i,:) = r(i,:) .* (xx-x(j))/x(i)-x(j)) 
    end 
end 
r = sum(r,1); 

Anstatt einen ganzen Vektor auf einmal zu transformieren, versuchen wir es mit nur einem Skalar. Verwende auch das einfachste Element von x, das von i und j abhängt. Ich habe auch die sum am Ende entfernt. Wir können es später wieder hinzufügen.

n=100; 
x=rand(1,n); 

r = ones(n,1); 
for i=1:n 
    for j=1:n 
     y = x(i)+x(j); 
     r(i) = r(i) * y; 
    end 
end 

Der obige Code ist die Beispielfunktion, ich parallelisieren wollen.

Die innere Schleife muss immer auf denselben Vektor r(i) für eine Iteration der äußeren Schleife zugreifen i. Dieser Zugriff ist ein Schreiben Betrieb (*=), aber die Bestellung ist nicht wichtig für diesen Vorgang.

Da verschachtelte parfor Schleifen in Matlab nicht erlaubt sind, habe ich versucht, alles in eine parfor Schleife zu packen.

n=100; 
x=rand(1,n); 

r = ones(n,1); 
parfor k=1:(n*n) 
    %i = floor((k-1)/n)+1; % outer loop 
    %j = mod(k-1,n)+1;  % inner loop 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(i) = r(i) * y;  % ERROR here 
end 

Da Indies berechnet werden, weiß Matlab noch nicht heiß, um es zu schneiden. Also entschied ich, die Multiplikation nach außen zu verschieben und lineare Indizes zu verwenden.

n=100; 
x=rand(1,n); 

r = ones(n,n); 
parfor k=1:(n*n) 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(k) = y; 
end 
r = prod(r,1); 
r = squeeze(r); % remove singleton dimensions 

Während dies nicht für skalare Werte in der inneren Schleife arbeiten, es funktioniert nicht für Vektoren, die in der inneren Schleife, da Indizes müssen neu berechnet werden.

Obwohl es funktioniert, wenn ich das Array umformen.

n=100; 
x=rand(1,n); 
m=5; 

r = ones(n*n,m); 
parfor k=1:(n*n) 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(k,:) = y.*(1:m); % ERROR here 
end 
r = reshape(r,n,n,m); 
r = prod(r,2); 
r = squeeze(r); % remove singleton dimensions 

Auf diese Weise kann ich r einen Vektor xx zu einem anderen Vektor zu transformieren.

n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = ones(n*n,m); 
parfor k=1:(n*n) 
    [j,i] = ind2sub([n,n],k); 
    y = x(i)+x(j); 
    r(k,:) = y.*xx; % ERROR here 
end 
r = reshape(r,n,n,m); 
r = prod(r,2); 
r = sum(r,1); 
r = reshape(r,size(xx)); % reshape output vector to input vector 

Für meine parallele Lösung, ich brauche eine n*n*m Array anstelle eines n*m Array, das ziemlich ineffizient zu sein scheint. Gibt es eine bessere Art zu tun, was ich will? Was sind die Vorteile anderer Möglichkeiten (schönerer Code, weniger CPU, weniger RAM, ...)?

UPDATE

In der Reihenfolge weggelassen zu versuchen, die Aufgabe und reduzieren sie auf das minimalen Arbeits Beispiel für das Problem, ich habe die Kontrolle von i~=j zu vereinfachen es einfacher zu machen, obwohl NaN Ergebnis in einem alle resultierenden . Außerdem führt die Art des Codes beim Hinzufügen dieser Überprüfung zu einem Ergebnis von 1. Damit der Code Sinn ergibt, sind die Faktoren nur Gewichte für einen anderen Vektor z.

Das aufwändigere Problem sieht wie folgt aus:

n=100; 
x=rand(1,n); 
z=rand(1,n); 
m=5; 
xx=rand(1,m); 

r = zeros(1,m); 
for i=1:n 
    q = ones(1,m); 
    for j=1:n 
     if i~=j 
      q = q .* (xx-x(j))/(x(i)-x(j)); 
     end 
    end 
    r = r + z(i) .* q; 
end 

Answer 1

Dieses Problem muss keine parallel zur Schleife auszuführen. Ein Problem ist, dass x(i)-x(j) viele Male redundant berechnet wird. Dies ist ineffizient. Der vorgeschlagene Ansatz berechnet jede Zahl genau einmal und vektorisiert die Operationen für jedes Element in xx. Da xx bei weitem der kürzeste Vektor ist, ist er fast vollständig vektorisiert. Wenn Sie auch die letzte Schleife vektorisieren wollen, wird dies wahrscheinlich nur eine versteckte for-Schleife sein, es wird viel mehr Speicher und der Code wäre komplizierter (wie 3D-Matrizen und so). Ich habe mir die Freiheit genommen, im Nenner Minus zu Plus zu wechseln, nur zum Testen. Minus würde NaN für alle Zahlen erzeugen. Der letzte Ansatz ist etwas schneller. Etwa 10 mal für n = 10000. Ich schlage vor, Sie versuchen ein bisschen mehr Benchmark.

function test() 
% Initiate variables 
n=100; 
x=rand(1,n); 
m=5; 
xx=rand(1,m); 

tic; 
% Alternative 1 
r = zeros(1,m); 
for i=1:n 
    q = ones(1,m); 
    for j=1:n 
     q = q .* (xx-x(j))/(x(i)+x(j)); 
    end 
    r = r + q; 
end 
toc; 

tic; 
% Alternative 2 
xden = bsxfun(@plus, x, x.'); % Calculate denominator 
xnom = repmat(x,n,1); % Calculate nominator 
xfull = (xnom./xden).'; % calculate right term on rhs. 

for (k = 1:m) 
    tmp= prod(xx(k)./xden - xfull); % Split in 2 calculations 
    r2(k) = sum(tmp); % "r = r + xx(k)" 
end 
toc; 

disp(r); 
disp(r2); 

Nur eine Anmerkung am Ende. Alternative 2 ist schneller, aber auch speicherintensiv. Im Falle von Speicherproblemen ist eine Schleife vorzuziehen. Außerdem sind im Falle der Parallelisierung keine globalen Variablen erforderlich. Für den Fall, dass Sie dies benötigen, müssen Sie wahrscheinlich über Ihr Design schauen (aber für den Fall, dass der Code kurz ist, gibt es keine kritischen, also sollten Sie sich nicht so viel Mühe machen).