Hi@all.

Ich habe ein großes "Array of Array of Array of double";

In diesem möchte ich per Finite-Differenzen-Methode die Lösung einer Poission Gleichung iterieren. Um eine Verbesserte Lösung an der Stelle (x,y,z) zu berechnen, muss auf alle Nachbarpunkte lesend, und auf die Stelle x,y,z schreibend zugegriffen werden.

Um das Verfahren zu beschleunigen, habe ich n Threads eingeführt.
Jeder Thread iteriert zunächst eine x-y-Ebene in diesem Array, erhöht dann z um 1 und iteriert die nächste x-y Ebene. Ist ein Thread beim Maximalwert von z angekommen, wird z wieder auf 0 gesetzt.
Damit ein Thread den anderen nicht "überholen" kann (Konsistenz der Lösung), fragt der i+1 Thread, ob die z Position des i-ten Thread mindestens um 1 größer ist, als die Position des i+1-ten Threads:

zusammenfalten · markieren

Delphi-Quellcode:

       while (not Iterrationsgenauigkeit_erreicht) do
       begin
        inc(j);
        for z := 0 to NodeCount_z-1 do
          begin
            current_z:=z;
            i:=Thread_id-1;
            if i=-1 then i:=Num_Thread-1;
            if assigned(Threads[i]) then
               while (Threads[i].Current_z>Current_z) and (Threads[i].Current_z-Current_z<=1) do
                  sleep(10);
            for x := 0 to NodeCount_x-1 do
               for y := 0 to NodeCount_y-1 do
                    begin
                        //Führe FDM-Iteration aus
                    end;
            end;
        end;

Zur kurzen und knappen Frage: Benötige ich irgendwo eine Critical-Section?
Da ich ja nicht auf den selben Array-Eintrag mit zwei verschiedenen Threads zugreife, gibts hier eig. kein Problem oder?
Benötige ich vielleicht für die variable current_z eine Critical-Section?

Besten Dank für eure Hilfe!

Euer Michael

Hat niemand eine Idee?

Zitat von alleinherrscher:

Da ich ja nicht auf den selben Array-Eintrag mit zwei verschiedenen Threads zugreife, gibts hier eig. kein Problem oder?

Wenn du das sicherstellen kannst dann brauchst du keine Critical Section.
Du könntest auch jedem Thread eine Kopie der Daten geben (ein Teilarray) und am Ende die Ergebnisse wieder zusammensetzen.
Das ist das Prinzip "Teile und Herrsche".
Wenn jeder Thread nur auf seinen eigenen Daten arbeitet dann kann er keinem anderen Thread in die Quere kommen.
Natürlich geht beim Kopieren von Arrays etwas Performance verloren; dann könnte aber im Vergleich zur Gesamtaufgabe unter einem Promille liegen.
Das Kopieren liese sich aber auch vermeiden, indem jeder Thread die Indexgrenzen (oder Zeiger auf Anfang und Ende) bekommen innerhalb dessen er arbeiten soll.
Das ist ja die Strategie die du z.Zt. verfolgst.
Bei harten Rechenaufgaben wie z.B. FFT sollte man übrigens nur soviele parallele Threads einsetzen wie der Prozessor an Kernen hat.

Hey, danke für die Infos. In der Tat hatte ich diese Methode (Teile und Herrsche) schon einmal implementiert. Allerdings geht aus folgendem Grund die performance sehr in den Keller:

1. An der Array-Grenze von zwei Threads, muss Thread 1 auf den Bereich von Thread 2 zugreifen und umgekehrt, da ich für die Berechnung am Punkt x,y,z jeweils die Punkte x+1,y,z ; x-1,y,z ; x,y+1,z ; x,y-1,z; x,y,z+1 ; x,y,z+1 benötige.

2. Ich muss sicherstellen, dass Thread 1 sein Teilgebiet nicht schon 200 mal durchlaufen hat, während Thread 2 erst 180 mal durchgelaufen ist.

3. Aus physikalischen Gründen benötigen die verschiedenen Teilgebiete unterschiedlich lange zur Berechnung. D.h.: Wegen dem vorher genannten Grund warten sich n-1 Thread die Füße in den Bauch, während Thread n noch rechnet....

Man muss die Teilstücke ja auch nicht starr einem Thread zuordnen, sondern kann einen Worker-Pool nutzen. Der enhält so viele Worker, wie es Prozessoren gibt.
Zu den Teilbereiche werden die angrenzenden Punkte mitgespeichert (sog. Geisterpunkte). Es sollte mehr Teilstücke geben als Worker.

Der Hauptthread macht nun folgendes:

1. die Teilstücke initialisieren
2. alle Teilstücke in Auftrag geben
3. blockieren, bis alle Teilstücke berechnet wurden
4. die Geisterpunkte aktualisieren
5. Abbruchkriterium prüfen, wenn nicht erfüllt: gehe zu 2.

Ein schöner Nebeneffekt ist, dass du nach Schritt 4. auch mal einen Zwischenstand speichern kannst.

Eine andere Sache: Wenn du weder Synchronisationsmittel nutzt, noch dich selbst darum kümmerst, können "merkwürdige" Sachen mit den Caches passieren. In den meisten Varianten sollte die Implementierung des Workerpool mit klassischen Synchronisationsmitteln schon ausreichen, um solche Effekte zu vermeiden.

Ah verstehe. Ja die Methode hört sich auch gut an

Zitat:

Eine andere Sache: Wenn du weder Synchronisationsmittel nutzt, noch dich selbst darum kümmerst, können "merkwürdige" Sachen mit den Caches passieren. In den meisten Varianten sollte die Implementierung des Workerpool mit klassischen Synchronisationsmitteln schon ausreichen, um solche Effekte zu vermeiden.

Ich sorge ja dafür dass nicht gleichzeitig auf gleiche Array-Einträge zugegriffen wird (mittels "current_z", s.o.)...hatte gehofft, dass das ausreicht?

Wenn du die Teilbereiche etwas intelligenter machst (kennt seine Nachbarn, kann sich und die Nachbarn sperren), dann erstellst du die Teilbereiche, packst diese in eine Queue und aus dieser Queue werden die Worker bedient.

markieren

Code:

repeat
  TB = Queue.Dequeue
  if not TB.TryLock then
    Queue.Enqueue( TB )
    TB = nil
until TB <> nil
WorkerThread.WorkOn( TB )

Dadurch organisieren sich die Teilbereiche und WorkerThreads selber und du kannst eine allgemein gültige Klasse schreiben, die solche Arten von Aufgaben durch Threads erledigen lässt.

Soll bei dieser Berechnung wirklich auf bereits verarbeitete Daten basiert werden? Ich spare mir gerade die Theorie dazu anzulesen, aber meistens sollte ein Durchgang auf einer Matrix auch komplett auf den Ursprungswerten erfolgen. Mal am Beispiel einer Weichzeichnung mit 3x3 Kernel á 1/9 (an den Rändern geclipped):

markieren

Code:

1.0  1.0  3.0
2.0  4.0  1.0
4.0  4.0  1.0
5.0  2.0  1.0

Spalte 1 von oben nach unten berechnen:

markieren

Code:

2.0  1.0  3.0
2.8  4.0  1.0
3.6  4.0  1.0
3.7  2.0  1.0

Spalte 2:

markieren

Code:

2.0  2.3  3.0
2.8  2.6  1.0
3.6  2.4  1.0
3.7  2.3  1.0

Spalte 3 spare ich mir mal. Eigentlich sollten sich hier aber alle Operationen auf die erste Matrix, und nur die beziehen:
Richtige Gesamtlösung:

markieren

Code:

2.0  2.0  2.3
2.6  2.3  2.3
3.5  2.7  2.2
3.8  2.8  2.0

Man würde also 2 Arrays brauchen: Eines aus dem gelesen wird, und eines in das die Ergebnisse kommen. Wenn mehrere Iterationen nötig sind, macht man im Folgeschritt einfach das Schreib-Array zum neuen Lese-Array. Dieser Weg hätte vor allem den riesen Vorteil, dass ruhig beliebig viele Threads auf dem Lese-Array wahlfrei lesen können, es muss nur sicher gestellt werden dass jeder Thread einen ganz eigenen Bereich im Schreib-Array bekommt. Keinerlei Synchronisierung nötig, abgesehen vom Fertig-Zeitpunkt aller Threads nach einem Durchgang. Dadurch wäre das ganze sogar so fein parallelisierbar, dass man eine GPU mit der Aufgabe betrauen könnte, und erschreckend große Performancegewinne hätte. (Defacto könnte für jeden Wert im Zielarray je ein Thread arbeiten, wenn dies für die Hardware günstig wäre.)

So kenne ich zumindest die meisten iterativen Verfahren auf Matrizen, sonst würden einem ja die Teilergebnisse einer Spalte/Zeile in die Basisdaten für die nächste "hinein bluten". (Was gewünscht sein könnte, aber eher selten ist. Daher frage ich so genau nach.) In 3D-Matrizen verhält sich das natürlich alles genau so.

Es sollten beide Möglichkeiten konvergieren (und afaik auch gleich schnell), also egal ob man nur von den alten Daten die neuen berechnet, oder die neuen auch dazu nimmt.
Je nach Dimension und Feinheit der Diskretisierung wird allerdings der Speicher schnell knapp, weswegen ein Kopie nicht gerade zu empfehlen ist:
Bsp: 10.000 Elemente in x und y Richtung -> 10^8 Elemente ~ 1 GB, interessant wirds dann in 3d (zumal die Matrix auch dicht besetzt ist/sein wird)

Edit: wobei ich mir jetzt nicht mehr ganz so sicher bin, ob es nicht doch einen (großen) Unterschied gibt. Wenn zB im Quadrat [0,1]x[0,1] oben und unten RB angegeben sind, so werden Störungen schneller in die eine als in die andere Richtung verbreitet (je nachdem wie die Elemente durchiteriert werden).

Was Notxor schreibt ist meiner Meinung nach richtig: Ich schreibe die neu berechneten Werte direkt wieder in mein Array zurück und berechne dann den nächsten Punkt mit all seinen Nachbarpunkten, unabhängig ob hier der Wert der n.ten oder n-1 -ten Berechnung gespeichert ist. Ein anderer Physiker, ein Professor in den Staaten verwendet die Methode mit 2 Arrays, hat aber Aufgrund von rotationssymmetrie eine Dimension weniger zu berechnen. Mein Problem ist allerdings nicht unbedingt rotationssymmetrisch. Wenn ich allerdings die selben Startbedingungen in meinem Programm setzt, wie er in seinem, kommt die selbe Lösung heraus. Daher denke ich, dass es für den Fall des Lösens einer Poisson gleichung egal ist, ob man immer das selbe Array benutzt, oder ein Lese- und ein Schreibarray. - Zumal man ja hier so lange iteriert, bis eine Steady-State Lösung herauskommt. Wenn ich mir also das Beispiel von Medium ansehe und vorstelle, dass ich unendlich oft über seine Matrix iteriere, wird dort auch das selbe Ergebnis herauskommen, egal ob ich 2 oder 1 Array verwende.

Um das ganze zu stabilisieren, habe ich noch eine Selbstkonsistenzprüfung eingebaut, welche nach der Iteration das Ergebnis so so lange zwischen dem neuen und dem vorhergehenden Wert hin und her schiebt, bis die Poisson Gleichung erfüllt (bzw. die Abweichung miniert) ist. (Verfahren des goldenen Schnitts).

In der Tat habe ich nicht genug Ram und zwei Arrays im Speicher zu halten.