Warning: Declaration of syntax_plugin_nbsp::handle($aMatch, $aState, $aPos, &$aHandler) should be compatible with DokuWiki_Syntax_Plugin::handle($match, $state, $pos, Doku_Handler $handler) in /homepages/13/d110133045/htdocs/foerstemann/dokuwiki/lib/plugins/nbsp/syntax.php on line 40

Warning: Declaration of syntax_plugin_nbsp::render($aFormat, &$aRenderer, $aData) should be compatible with DokuWiki_Syntax_Plugin::render($format, Doku_Renderer $renderer, $data) in /homepages/13/d110133045/htdocs/foerstemann/dokuwiki/lib/plugins/nbsp/syntax.php on line 40

Warning: Cannot modify header information - headers already sent by (output started at /homepages/13/d110133045/htdocs/foerstemann/dokuwiki/lib/plugins/nbsp/syntax.php:40) in /homepages/13/d110133045/htdocs/foerstemann/dokuwiki/inc/actions.php on line 207

Warning: Cannot modify header information - headers already sent by (output started at /homepages/13/d110133045/htdocs/foerstemann/dokuwiki/lib/plugins/nbsp/syntax.php:40) in /homepages/13/d110133045/htdocs/foerstemann/dokuwiki/lib/tpl/dokuwiki/main.php on line 12
algorithmus_zur_numerischen_flaechenabschaetzung_der_mandelbrotmenge [name://Thorsten.Förstemann/Labor]

Benutzer-Werkzeuge

Webseiten-Werkzeuge


algorithmus_zur_numerischen_flaechenabschaetzung_der_mandelbrotmenge

Algorithmus zur Flächenabschätzung der Mandelbrotmenge

Die hier vorgenommene Abschätzung des Flächeninhalts der Mandelbrotmenge erfolgt mit Hilfe von Quasi-Monte-Carlo-Methoden (Pixel Counting). Im Prinzip wird ein hochaufgelöstes Schwarz/Weiss-Bild der Mandelbrotmenge gezeichnet und dann werden die Pixel ausgezählt.

Komplexität wächst dabei etwa quadratisch: eine Dezimalstelle mehr bedeutet etwa 100-facher Aufwand. Dies kann durch Verwendung von CUDA oder Open/CL kompensiert werden: Monte-Carlo-Methoden lassen sich in der Regel gut Parallelisieren und man kann Geschwindigkeitszuwächse um ca. Faktor 100 erreichen.

Umsetzung erfolgte mit Mathematica 8 auf einem Intel Core i7 mit Radeon HD 5970.

Erzeugung des Rasters

Die hier verwendete Monte-Carlo-Methode ist auf eine effiziente Erzeugung sehr vieler Pixel bzw. Koordinaten angewiesen. Dazu werden die Koordinaten durch ein ij-Raster erzeugt, das immer wieder mit einem Block-Offset Onm und einem Random-Offset RN verschoben wird.

Die Größe des ij-Rasters wird durch den Grafikspeicher begrenzt. Der Block-Offset sorgt für handliche Arbeitspakete (1 Stunde Rechenzeit) und der Random-Offset sorgt für einen nicht versiegenden Strom an Arbeitspaketen.

Die Berechnung der Koordinaten erfolgt in der GPU und kostet kaum Rechenzeit.

Für die komplexen Koordinaten, die bei der Berechnung der Mandelbrotmenge auftauchen, gilt im Folgenden z = x + iy.

Das ij-Raster

Figure 2: The ij-Raster (marked gray). The block offset O<sub>nm</sub> (marked red) is drawn on a lager scale.

Das ij-Raster hat eine Größe von 4096×4096 Pixel. Größere Raster lassen sich nicht zuverlässig in den Speicher der Grafikkarte laden.

Das ij-Raster hat eine Seitenlänge von 2.6, die Rasterweite beträgt also w = 2.6/4096. In der Abbildung 2 ist ein ij-Raster durch den grauen Kasten symbolisiert.

Es gilt für die Koordidinaten des ( i , j )-ten Rasterpunkts mit Rasterweite w:

Xij = ( xij , yij ) = w * ( i , j ) + (-2.1 , -1.3 ).

Der kleine rote Kasten symbolisiert den Block-Offset Onm. Der rote Kasten ist nicht maßstabsgetreu gezeichnet: Die Seitenlänge beträgt nur die Rasterweite w des ij-Rasters.

Der Block-Offset

Figure 3: The block offsets O<sub>nm</sub> (marked red). The random offset is marked green.

Die Größe des ij-Rasters wurde durch den Speicher der Grafikkarte bestimmt. Ein handliches Arbeitspaket besteht aus 32×32 ij-Rastern – hierfür braucht die eingesetzte Hardware ca. eine Stunde Rechenzeit.

Um den Arbeitsaufwand auf die ij-Raster gleich zu verteilen, werden diese um einen Bruchteil ihrer Rasterweite w verschoben – die verschobenen ij-Raster zeigen dann im Prinzip den gleichen Bildinhalt. In der Abbildung sind die verschiedenen Block-Offsets Onm rot dargestellt. Der Ausschnitt eines ij-Raster an Position (n,m)=(20,25) ist beispielhaft eingezeichnet. Die Verschiebung beträgt Vielfache von w/32 und es gilt:

Onm = w/32 * ( n , m )

Der kleine güne Kasten symbolisiert den Random-Offset RN. Der grüne Random-Offset ist maßstabsgetreu gezeichnet.

Der Random-Offset

Figur 4: The random offsets R<sub>N</sub> (marked green)

Ein Arbeitspaket besteht aus 32×32 ij-Rastern. Die Rasterweite des Arbeitspakets beträgt damit w/32. Jedes Arbeitspaket wird um einen zufälligen Random-Offset RN innerhalb der Rasterweite w/32 des Arbeitspakets verschoben. Die einzelnen Random-Offsets RN werden durch eine 2-dimensionale Niederreiter-Sequenz1) erzeugt. Man erkennt, dass die Seitenlänge des Random-Offset-Bereichs 0.000 019 8 = 2.6/4096/32 beträgt. Es gilt für das N-te Arbeitspaket:

RN = w/32 Niedereiter( N , dim=2 ) - w/64

In der Abbildung sind die verwendeten Random-Offsets grün dargestellt. Die Arbeitspakete sind also nur leicht verschoben und besitzen den gleichen zu erwartenden Rechenaufwand.

Jeder grüne Punkt entspricht einem blauen Punkt aus der ersten Abbildung 1.

Insgesamt gilt also für die Koordinaten des ( i , j )-ten Rasterpunkts mit ( n , m )-Offset im N-ten Arbeitspaket:

XNnmij = RN + Onm + Xij

mit w = 2.6/4096.

Berechnung der Iterationen

Wie oben beschrieben wird die Aufgabe in sich wiederholende Arbeitspakete zerlegt. Beim Abarbeiten eines Arbeitspaketes wird in einem ersten Schritt schwerpunktmäßig die GPU auslastet. Dann wird die restliche Arbeit in einem zweiten Schritt von der CPU erledigt, während die GPU schon den ersten Schritt des nächsten Arbeitspaketes bearbeitet. Um die Hardware optimal auszulasten, ist der Zeitpunkt der Übergabe so zu wählen, dass beide Teile etwa gleiche Laufzeit besitzen.

Der GPU-Teil besteht dabei aus nochmals zwei Unterabschnitten mit unterschiedlichen Methoden. Im ersten GPU-Unterabschnitt wird eine grobe Vorselektion der ij-Raster vorgenommen. Im zweiten Unterabschnitt werden dann Koordinatenlisten nach und nach ausgedünnt. Der erste Unterabschnitt arbeitet mit Bitmaps während der zweite mit Listen arbeitet.

Kernproblem der folgenden Überlegungen ist, die zahlreichen GPUs in Gegenwart des sehr inhomogenen (chaotischen) Verhaltens der Mandelbrotmenge gleichmäßig hoch auszulasten. Mit einem Pixel mit 109 Iterationen muss man anders umgehen als mit 109 Pixel mit einer Iteration, wenn man effizient sein will.

Hilfsfunktionen

Mit den folgenden Funktionen kann für einen Pixel geprüft werden, ob er zur Mandelbrotmenge gehört. Am Ende der Prüfung besitzt ein Pixel eine der folgenden Eigenschaften:

  • member: Der Pixel gehört zur Mandelbrotmenge.
  • refugee: Der Pixel gehört nicht zur Mandelbrotmenge.
  • unknown: Es kann nicht entschieden werden, ob der Pixel zur Mandelbrotmenge gehört.

Es werden folgende recht triviale Hilfsfunktionen benötigt, die auf der GPU laufen:

  • iter( n ) berechnet die Iteration eines Pixels bis zur maximalen Iterationstiefe n. Dies ist die typische Mandelbrot-Iteration2).
    Ergebnis: entweder refugee oder unknown
  • test( ) prüft, ob der Pixel Teil der Kardioide oder des Kreises ist3).
    Ergebnis: entweder member oder unknown
  • orbit( n ) berechnet die Iterationen eines Pixels mit Orbit-Detektion4). Die Funktion ist durch eine extra Prüfung etwas aufwändiger als iter( n ).
    Ergebnis: member, refugee oder unknown

Werden zum Beispiel die Hilfsfunktionen iter( n ) und iter( m ) hintereinander ausgeführt, so ist das identisch mit iter( n + m ).

Pixel-Vorselektion (Map Selection)

Mit Hilfe der Hilfsfunktionen und einer Funktion list( ) kann die Vorselektion geschrieben werden als:

  iter( n )
  test( )
  loop ( e = 1..m )
     iter ( n * pow(2,e) )
     orbit ( n * pow(2,e+1) )
  list( )

Die Funktion list( ) wandelt unknown Pixel aus der Bitmap in eine Koordinatenliste und läuft auf der CPU. Dabei gilt n = 64 und m = 2. Diese Werte sind Erfahrungswerte. Die maximale Iterationstiefe für diese Parameter beträgt 64 + 128 + 256 + 256 + 1024 = 1728.

Für ein ij-Raster mit 40962 Pixel benötigt die Vorselektion ca. 0.3 Sekunden auf der GPU. Die Funktion list( ) benötigt ca. 0.5 Sekunden auf der CPU. Die Koordinatenliste wird für die folgenden Berechnungen benötigt.

In ca. 0.8 Sekunden werden die 40962 Pixel in typischerweise 20.9 % member, 77.6 % refugee und 1.5 % unkown Pixel aufgeteilt. Der nächste Schritt behandelt nur noch die verbleibenden ca. 1.5 % * 40962 = 250 000 unknown Pixel.

Koordinatenlisten-Ausdünnung (List Dieout)

Mit Hilfe der Hilfsfunktionen und einer Funktion dieout( ) kann die Listen-Ausdünnung geschrieben werden als:

  loop ( e = 0..m )
     iter ( n * pow(2,e) )
     orbit ( n * pow(2,e+1) )
     dieout( )

Die Funktion dieout( ) entfernt alle member und refugee Pixel aus der Koordinatenliste und läuft auf der CPU. Dabei gilt n = 1024 und m = 11. Diese Werte sind Erfahrungswerte. Die maximale Iterationstiefe für diese Parameter beträgt 12 579 840. Dieser Wert ergibt sich als Summe in der Spalte maximum iteration in Tabelle 1.

In Tabelle 1 sind einige beispielhafte Parameter während des loops für e = 0..m angegeben. Gestartet wird mit den ca. 250 000 unknown Pixeln aus der Vorselektion.

Table 1: typical parameters while running the loop ( e = 0..m )
loop found member found refugee unknown maximum iteration total iteration time Giga it./sec.
start 0 0 240 846 0 0 0 sec. 0
e = 0 + 128 967 + 13 885 97 994 3 072739 878 912 0.5 sec. 1.4
e = 1 + 64 437 + 2 987 30 570 6 144602 075 136 0.3 sec. 2.0
e = 2 + 16 702 + 799 13 069 12 288375 644 160 0.2 sec. 1.9
e = 3 + 6 746 + 311 6 012 24 576321 183 744 0.1 sec. 3.2
e = 4 + 3 003 + 154 2 855 49 152295 501 824 0.1 sec. 2.9
e = 5 + 1 318 + 66 1 471 98 304280 657 920 0.1 sec. 2.8
e = 6 + 743 + 25 703 196 608289 210 368 0.1 sec. 2.8
e = 7 + 349 + 9 345 393 216276 430 848 0.1 sec. 2.7
e = 8 + 175 + 6 164 786 432271 319 040 0.2 sec. 1.3
e = 9 + 77 + 4 83 1 572 864257 949 696 0.2 sec. 1.2
e = 10 + 41 + 3 39 3 145 728261 095 424 0.4 sec. 0.7
e = 11 + 18 + 1 20 6 291 456245 366 784 0.9 sec. 0.3

Erläuterungen zu den Spalten in Tabelle 1:

  • loop: Parameter e
  • member: gefundene member Pixel
  • refugee: gefundene refugee Pixel
  • unknown: verbleibende unknown Pixel – mit diesen wird weitergerechnet.
  • maximum iteration: maximale Iterationen, n * pow(2,e) + n * pow(2,e+1)
  • total iteration: Anzahl berechneter Iterationen, unknown(e-1) x maximum iteration(e)
  • time: GPU/CPU-Zeit (ungenau)
  • Giga it./sec.: Giga Iterationen pro Sekunde

Man sieht, dass die Iterationen pro Sekunde für kleine und große e nicht optimal sind. Bei kleinen e ist die Koordinatenliste sehr lang (Spalte unknown) und die Funktion dieout( ) kostet viel Zeit. Für große e wird die Liste immer kürzer und die GPU kann nicht mehr mit Pixeln ausreichend befüllt werden. Hier sind mit Sicherheit noch Optimierungen möglich.

Zusammen mit der Vorselektion benötigt die Listen-Ausdünnung ca. 4 Sekunden, um die 40962 Pixel des ij-Rasters auf ca. 20 unknown Pixel zu reduzieren. Ein Arbeitspaket mit 32×32 ij-Rastern benötigt also ca. 1 Stunde. Nach einer Stunde werden die Koordinaten der ca. 32 * 32 * 20 = 20480 unknown Pixel in eine Datei geschrieben, die von der Nach-Iteration weiterverarbeitet werden.

Randpunkt-Nachiteration (Post Iteration)

Jede Stunde fallen ca. 20 000 unknown Pixel nach der Listen-Ausdünnung an. Diese Pixel sind

  • weder nach 12 579 840 Iterationen divergiert (refugee Pixel)
  • noch sind die Pixel nach spätestens 6 288 384 Iterationen in einen Orbit mit maximaler Länge 6 291 456 eingetreten (member Pixel).

Mit Hilfe der Hilfsfunktionen kann die Nachiteration geschrieben werden als:

  loop ( e = 0..m )
     iter ( n1 * pow(2,e) )
     orbit ( n2 * pow(2,e) )

Dabei gilt n1 = 3 200 000 und n2 = 4 194 304 = 222 und m = 9. Diese Werte sind Erfahrungswerte. Die maximale Iterationstiefe für diese Parameter beträgt 7 564 372 992. Ein Orbit wird erkannt, wenn seine Länge höchstens 2 147 483 648 beträgt und spätestens nach 5 416 889 344 Iterationen begonnen hat (siehe dazu auch Tabelle 3 im Abschnitt Fehlerbetrachtungen).

Von den ca. 20 000 unknown Pixel nach der Koordinatenlisten-Ausdünnung bleiben nach der Nachiteration noch ca. 15 unknown Pixel. Dafür wird ca. 1 Stunde CPU-Zeit benötigt im Parallelbetrieb mit der Vorselektion und der Listen-Ausdünnung. Die Anzahl der member und refugee Pixel wird mit den Koordinaten der unknown Pixel für eine spätere Auswertung in eine Datei geschrieben.

Typische Parameter für ein Arbeitspaket mit 22×17 Pixel während des gesamten Algorithmus sind in Tabelle 2 angegeben. Man erkennt, dass der Anteil der refugee Pixel an den verbliebenden unknown Pixel vermutlich 5 % nicht übersteigt.

Table 2: typical parameters of an raster with 22×17 pixels
task found members found refugees unknown pixels GPU time CPU time
start 0 0 16 400 000 000 0 min. 0 min.
map select (GPU) + 3 420 000 000 + 12 700 000 000 245 000 000 4 min. 7 min.
list dieout (GPU) + 223 000 + 18 000 20 000 50 min. 3 min.
post iteration (CPU) + 19 175 + 810 15 0 min. 50 min.
algorithmus_zur_numerischen_flaechenabschaetzung_der_mandelbrotmenge.txt · Zuletzt geändert: 2015/09/10 10:11 von vaumann