SIMD optimalizace Mandelbrotovy AVX2 OpenMP CUDA

# Optimalizace výpočtu Mandelbrotovy množiny: od skalárního kódu k SIMD, OpenMP a CUDA

Mandelbrotova množina se počítá pro každý bod komplexní roviny nezávisle: ověřuje se omezenost posloupnosti z_{n+1} = z_n² + c při z_0 = 0. Pokud |z_n| > 2, posloupnost se rozbíhá. Omezením iterací na 256 lze výpočet dokončit. Tato úloha je ideálně paralelizovatelná – každý pixel se zpracovává samostatně.

Barva pixelu závisí na čísle iterace, při které |z| překročil poloměr 2. Body, které za 256 kroků neopustily oblast, se vybarví černě.

Skalární implementace

Základní algoritmus v C++ převádí souřadnice obrazovky [0, WIDTH] × [0, HEIGHT] do oblasti [-1, 1] × [-1, 1] komplexní roviny:

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

Smyčka iterací je optimalizována ukládáním z_x² a z_y² odděleně:

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++;
}

Na AMD Ryzen 5 5600H při 1920×1080 a 256 iteracích dává skalární kód s -O2 7.0 ± 0.1 FPS (GCC/Clang).

Vektorizace AVX2

SIMD na 256bitových registrech ymm zpracovává 8 float za instrukci. Intrinsiky <x86intrin.h> umožňují explicitní řízení vektorizace.

Klíčové vektory:

• _c_x: základní souřadnice + posuny [0..7] × c_step_x
• _z_x, _z_y, _z_x2, _z_y2, _z_xy: stavy pro 8 pixelů
• _iterations: čítače iterací (int32)

Inicializace posunů:

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));

Ve smyčce iterací:

• Výpočet poloměru: _radius2 = _mm256_add_ps(_z_x2, _z_y2)
• Maska pokračování: _mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OS) → vše -1.0 nebo 0.0
• Aktualizace stavů podle aktivní masky
• Posun podle X po zpracování 8 pixelů: _c_x = _mm256_add_ps(_c_x, _8_c_steps_x)

Maska umožňuje ukončit smyčku pro všech 8 pixelů najednou, pokud všichni opustili orbitu.

Mnohopotocnost OpenMP

OpenMP rozděluje řádky mezi vlákna direktivou #pragma omp parallel for. Každé vlákno samostatně počítá svůj pruh pixelů.

#pragma omp parallel for schedule(dynamic)
schedule(dynamic) pro vyvážení zátěže — různé oblasti fraktálu vyžadují různý počet iterací.

Kombinace AVX2 + OpenMP na 6 jádrech Ryzen 5 5600H dává nárůst na 40-50 FPS.

Zrychlení na GPU s CUDA

Přenos na CUDA je triviální: každý vlákno počítá jeden pixel. Jádro CUDA:

__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;
}

Spouštění: dim3 block(16,16); dim3 grid((width+15)/16, (height+15)/16); mandelbrot_cuda<<<grid,block>>>(...)

Na GPU jednoho notebooku (RTX 3050) se dosahuje 500+ FPS při Full HD.

| Implementace | FPS (Ryzen 5 5600H) |

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

| Skalární | 7.0 |

| AVX2 | 45.2 |

| AVX2+OpenMP | 48.1 |

| CUDA | 520+ |

Co je důležité

• Mandelbrotova množina – referenční embarrassingly parallel úloha pro testování SIMD, OpenMP, CUDA
• AVX2 dává ~6.5násobné zrychlení oproti skalárnímu kódu díky zpracování 8 pixelů za cyklus
• OpenMP přidává 6-7% nárůst na 6 jádrech kvůli overheadu synchronizace
• CUDA zajišťuje řádový nárůst díky tisícům paralelních vláken
• Měření výkonu probíhá bez vykreslování (hyperfine, N spuštění ~5s)

— Editorial Team