# Optymalizacja obliczania zbioru Mandelbrota: od kodu skalarnego do SIMD, OpenMP i CUDA
Zbiór Mandelbrota jest obliczany dla każdej tochki płaszczyzny zespolonej niezależnie: sprawdza się ograniczoność ciągu z_{n+1} = z_n² + c przy z_0 = 0. Jeśli |z_n| > 2, ciąg rozbiega się. Ograniczenie iteracji do 256 pozwala zakończyć obliczenia. Zadanie idealnie nadaje się do paralelizowania — każdy piksel jest przetwarzany autonomicznie.
Kolor piksela zależy od numeru iteracji, w której |z| przekroczył promień 2. Punkty, które nie wyszły poza granice w 256 krokach, są malowane na czarno.
Implementacja skalarna
Podstawowy algorytm w C++ przekształca współrzędne ekranu [0, WIDTH] × [0, HEIGHT] w obszar [-1, 1] × [-1, 1] płaszczyzny zespolonej:
float c_y = -1.0f + screen_y * (2.0f / WINDOW_HEIGHT);
float c_x = -1.0f + screen_x * (2.0f / WINDOW_WIDTH);Pętla iteracji jest zoptymalizowana poprzez przechowywanie z_x² i z_y² oddzielnie:
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 przy 1920×1080 i 256 iteracjach kod skalarny z -O2 daje 7.0 ± 0.1 FPS (GCC/Clang).
Wektoryzacja AVX2
SIMD na 256-bitowych rejestrach ymm przetwarza 8 float na instrukcję. Instrynystyki <x86intrin.h> pozwalają na jawne sterowanie wektoryzacją.
Kluczowe wektory:
_c_x: bazowa współrzędna + przesunięcia [0..7] × c_step_x_z_x,_z_y,_z_x2,_z_y2,_z_xy: stany dla 8 pikseli_iterations: liczniki iteracji (int32)
Inicjalizacja przesunięć:
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));W pętli iteracji:
- Obliczenie promienia:
_radius2 = _mm256_add_ps(_z_x2, _z_y2) - Maska kontynuacji:
_mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OS)→ wszystkie -1.0 lub 0.0 - Aktualizacja stanów przy aktywnej masce
- Przesunięcie po X po przetworzeniu 8 pikseli:
_c_x = _mm256_add_ps(_c_x, _8_c_steps_x)
Maska pozwala przerwać pętlę dla wszystkich 8 pikseli jednocześnie, jeśli wszystkie opuściły orbitę.
Wielowątkowość OpenMP
OpenMP rozdziela wiersze między wątki za pomocą dyrektywy #pragma omp parallel for. Każdy wątek niezależnie oblicza swój pasek pikseli.
#pragma omp parallel for schedule(dynamic)schedule(dynamic) do równoważenia obciążenia — różne obszary fraktala wymagają różnej liczby iteracji.
Kombinacja AVX2 + OpenMP na 6 rdzeniach Ryzen 5 5600H daje przyrost do 40-50 FPS.
Przyspieszenie GPU CUDA
Przeniesienie na CUDA jest trywialne: każdy wątek oblicza jeden piksel. 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;
}Uruchomienie: dim3 block(16,16); dim3 grid((width+15)/16, (height+15)/16); mandelbrot_cuda<<<grid,block>>>(...)
Na GPU jednego laptopa (RTX 3050) osiąga się 500+ FPS przy Full HD.
| Implementacja | FPS (Ryzen 5 5600H) |
|---------------|----------------------|
| Skalarna | 7.0 |
| AVX2 | 45.2 |
| AVX2+OpenMP | 48.1 |
| CUDA | 520+ |
Co ważne
- Zbiór Mandelbrota — wzorcowe zadanie embarrassingly parallel do testowania SIMD, OpenMP, CUDA
- AVX2 daje ~6.5x przyspieszenie w porównaniu do kodu skalarnego dzięki przetwarzaniu 8 pikseli na cykl
- OpenMP dodaje 6-7% przyrostu na 6 rdzeniach z powodu narzutu synchronizacji
- CUDA zapewnia przyspieszenie rzędu wielkości dzięki tysiącom równoległych wątków
- Pomiar wydajności przeprowadzono bez renderowania (hyperfine, N uruchomień ~5s)
— Editorial Team
Brak komentarzy.