SIMD-Optimierung von Mandelbrot AVX2 OpenMP CUDA

Optimierung der Mandelbrot-Menge-Berechnung: Vom skalaren Code zu SIMD, OpenMP und CUDA

Die Mandelbrot-Menge wird unabhängig für jeden Punkt in der komplexen Ebene berechnet: Wir prüfen, ob die Folge z_{n+1} = z_n² + c begrenzt bleibt, beginnend mit z_0 = 0. Wenn |z_n| > 2, divergiert die Folge. Eine Begrenzung der Iterationen auf 256 sorgt dafür, dass die Berechnung terminiert. Die Aufgabe ist perfekt parallelisierbar – jedes Pixel wird unabhängig verarbeitet.

Die Pixelfarbe hängt von der Iterationszahl ab, bei der |z| den Radius 2 überschreitet. Punkte, die nach 256 Schritten noch innerhalb der Grenzen sind, werden schwarz dargestellt.

Skalare Implementierung

Der grundlegende Algorithmus in C++ bildet Bildschirmkoordinaten [0, WIDTH] × [0, HEIGHT] auf den Bereich der komplexen Ebene [-1, 1] × [-1, 1] ab:

float c_y = -1.0f + screen_y * (2.0f / WINDOW_HEIGHT);
float c_x = -1.0f + screen_x * (2.0f / WINDOW_WIDTH);

Die Iterationsschleife ist optimiert, indem z_x² und z_y² separat gespeichert werden:

while (z_x2 + z_y2 < MAX_RADIUS_2 && iterations < MAX_ITERATION_DEPTH) {
    z_y = 2 * z_x * z_y + c_y;
    z_x = z_x2 - z_y2 + c_x;
    z_x2 = z_x * z_x;
    z_y2 = z_y * z_y;
    iterations++;
}

Auf einem AMD Ryzen 5 5600H bei 1920×1080 mit 256 Iterationen liefert der skalare Code mit -O2 7,0 ± 0,1 B/s (GCC/Clang).

AVX2-Vektorisierung

SIMD auf 256-Bit-ymm-Registern verarbeitet 8 Floats pro Instruktion. Intrinsics aus <x86intrin.h> bieten explizite Kontrolle über die Vektorisierung.

Wichtige Vektoren:

• _c_x: Basis-Koordinate + Offsets [0..7] × c_step_x
• _z_x, _z_y, _z_x2, _z_y2, _z_xy: Zustände für 8 Pixel
• _iterations: Iterationszähler (int32)

Offset-Initialisierung:

const __m256 _01234567 = _mm256_set_ps(7.0f, 6.0f, 5.0f, 4.0f, 3.0f, 2.0f, 1.0f, 0.0f);
_c_x = _mm256_add_ps(_c_x, _mm256_mul_ps(_c_step_x, _01234567));

In der Iterationsschleife:

• Radius berechnen: _radius2 = _mm256_add_ps(_z_x2, _z_y2)
• Fortsetzungs-Maske: _mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OS) → alles -1.0 oder 0.0
• Zustände unter der aktiven Maske aktualisieren
• Nach Verarbeitung von 8 Pixeln entlang X verschieben: _c_x = _mm256_add_ps(_c_x, _8_c_steps_x)

Die Maske ermöglicht es, die Schleife für alle 8 Pixel gleichzeitig zu beenden, sobald sie alle entkommen sind.

OpenMP-Multithreading

OpenMP verteilt Reihen auf Threads mithilfe der Direktive #pragma omp parallel for. Jeder Thread berechnet unabhängig seinen Streifen von Pixeln.

#pragma omp parallel for schedule(dynamic)

schedule(dynamic) für Lastverteilung – unterschiedliche Fraktalregionen benötigen variierende Iterationszahlen.

Die Kombination aus AVX2 + OpenMP auf dem 6-Kern-Ryzen 5 5600H liefert bis zu 40–50 B/s.

GPU-Beschleunigung mit CUDA

Die Portierung nach CUDA ist einfach: Jeder Thread berechnet ein Pixel. CUDA-Kernel:

__global__ void mandelbrot_cuda(float* output, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= width || y >= height) return;
    
    float c_x = -2.0f + (float)x * scale_x;
    float c_y = -1.5f + (float)y * scale_y;
    
    float z_x = 0.0f, z_y = 0.0f;
    int iter = 0;
    
    while (z_x*z_x + z_y*z_y < 4.0f && iter < 256) {
        float temp = z_x*z_x - z_y*z_y + c_x;
        z_y = 2.0f * z_x * z_y + c_y;
        z_x = temp;
        iter++;
    }
    output[y*width + x] = (float)iter;
}

Start: dim3 block(16,16); dim3 grid((width+15)/16, (height+15)/16); mandelbrot_cuda<<<grid,block>>>(...)

Auf einer Laptop-GPU (RTX 3050) erreicht es bei Full HD über 500 B/s.

| Implementierung | B/s (Ryzen 5 5600H) |

|---------------|----------------------|

| Skalar | 7,0 |

| AVX2 | 45,2 |

| AVX2+OpenMP | 48,1 |

| CUDA | 520+ |

Wichtige Erkenntnisse

• Die Mandelbrot-Menge ist ein Benchmark, der sich mühelos parallelisieren lässt, ideal zum Testen von SIMD, OpenMP und CUDA
• AVX2 bietet ~6,5-fache Beschleunigung gegenüber skalarem Code durch Verarbeitung von 8 Pixeln pro Schleifendurchlauf
• OpenMP bringt auf 6 Kernen nur 6–7 % Zuwachs wegen Synchronisationsoverhead
• CUDA liefert eine um Größenordnungen höhere Geschwindigkeit dank Tausender paralleler Threads
• Leistungsmessungen schließen das Rendering aus (hyperfine, N Läufe ~5 s)

— Editorial Team