Können Sie eine inexakte Jacobi-Berechnung in Julia ohne spezielle Arrays parallelisieren?

Question

In Julia wollten Sie eine ungenaue Jacobi basierend auf einer Vektorfunktion berechnen, f (x), die viel Rechenleistung benötigt, um zu bewerten. Die Bewertung der Jacobi ist offensichtlich ziemlich klar im Konzept parallelisierbar. Meine Frage ist, kann dies in Julia getan werden, ohne DistributedArray, SharedArray, etc?

Beispiel: Angenommen, Sie den Code haben:

function Jacob(f::Function,x) 
    eps=1e-7 
    delta=eps*eye(length(x)) 
    J=zeros(length(x),length(x)) 
    for i=1:length(x) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    end 
    J 
end 

Ist es möglich, dies auf die gleiche Weise parallelisieren, die Sie vielleicht die Summe von 200 Millionen Zufall Münzwürfe, gemäß dem Handbuch parallelisieren? Das heißt, etwas zu äquivalente

nheads = @parallel (+) for i=1:200000000 
    int(randbool()) 
end 

Ich habe versucht, dies:

function Jacob(f::Function,x) 
    require("testfunc.jl"); 
    eps=1e-7 
    delta=eps*eye(length(x)) 
    J=zeros(length(x),length(x)) 
    J=@parallel (+) for i=1:length(x) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    J 
    end 
    J 
end 

wo „testfunc.jl“ ist der Name der Datei, in diesem Code und die Definition von f selbst gefunden werden . Wenn ich das versuchte, mit f einfach auswerten zu x.^2 + cos (x), war ich in der Lage, eine richtige (diagonale) Matrix herauszuholen, aber die Werte stimmten nicht mit denen überein, die durch den nichtparallelen Code gegeben wurden kann bestätigen, dass sie die richtigen Werte haben). Weitere Untersuchungen deuten darauf hin, dass die resultierende Jacobi einige Werte mit 2 oder 3 multipliziert hat, wenn julia -p 4.

Ist der Ansatz, den ich beschrieben habe, eine plausible (und erfordert nur Tweaking, um Doppelungen von Auswertungen zu verhindern)? Wenn nicht, gibt es eine andere Methode, mit der ich die Jacobi ohne die Verwendung der komplizierteren speziellen Array-Typen auswerten kann?

Es scheint, dass das Hinzufügen von "J = Nullen (n, n)" als erste Operation innerhalb der parallelen for-Schleife dieses Duplikationsproblem korrigiert. Kann dasselbe getan werden, ohne auf solch brutales Ausräumen des J-Arrays zurückzugreifen?

Answer 1

Was ich von oben Code zu verstehen ist, dass, wenn Sie schreiben:

J=zeros(length(x),length(x)) 
    J=@parallel (+) for i=1:length(x) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    J 
    end 

Julia sendet eine Kopie von J zu neuen Prozess wertet dann f(x) und Summe zusammen führt. ich denke, die bessere und effizientere Art und Weise zu senden J zwischen Threads zu verhindern ist, und wie folgt vorgehen:

@parallel (+) for i=1:length(x) 
    J=zeros(length(x),length(x)) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    J 
    end 

Mit dem obigen Code jeden Thread auf einer frischen J arbeitet und so Summe gibt die richtige Antwort.