SIMD optimization of Mandelbrot AVX2 OpenMP CUDA

# Optimizing Mandelbrot Set Computation: From Scalar Code to SIMD, OpenMP, and CUDA

The Mandelbrot set is computed independently for each point in the complex plane: we check whether the sequence z_{n+1} = z_n² + c remains bounded, starting with z_0 = 0. If |z_n| > 2, the sequence diverges. Capping iterations at 256 ensures the computation finishes. The task is perfectly parallelizable—each pixel is processed independently.

Pixel color depends on the iteration count at which |z| exceeds radius 2. Points that stay within bounds after 256 steps are colored black.

Scalar Implementation

The basic algorithm in C++ maps screen coordinates [0, WIDTH] × [0, HEIGHT] to the complex plane region [-1, 1] × [-1, 1]:

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

The iteration loop is optimized by storing z_x² and z_y² separately:

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

On an AMD Ryzen 5 5600H at 1920×1080 with 256 iterations, the scalar code with -O2 delivers 7.0 ± 0.1 FPS (GCC/Clang).

AVX2 Vectorization

SIMD on 256-bit ymm registers processes 8 floats per instruction. Intrinsics from <x86intrin.h> provide explicit control over vectorization.

Key vectors:

• _c_x: base coordinate + offsets [0..7] × c_step_x
• _z_x, _z_y, _z_x2, _z_y2, _z_xy: states for 8 pixels
• _iterations: iteration counters (int32)

Offset initialization:

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 the iteration loop:

• Compute radius: _radius2 = _mm256_add_ps(_z_x2, _z_y2)
• Continuation mask: _mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OS) → all -1.0 or 0.0
• Update states under the active mask
• Shift along X after processing 8 pixels: _c_x = _mm256_add_ps(_c_x, _8_c_steps_x)

The mask allows breaking the loop for all 8 pixels simultaneously if they've all escaped.

OpenMP Multithreading

OpenMP distributes rows across threads using the #pragma omp parallel for directive. Each thread independently computes its strip of pixels.

#pragma omp parallel for schedule(dynamic)

schedule(dynamic) for load balancing—different fractal regions require varying numbers of iterations.

Combining AVX2 + OpenMP on the 6-core Ryzen 5 5600H yields up to 40-50 FPS.

GPU Acceleration with CUDA

Porting to CUDA is straightforward: each thread computes one 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;
}

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

On a laptop GPU (RTX 3050), it achieves 500+ FPS at Full HD.

| Implementation | FPS (Ryzen 5 5600H) |

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

| Scalar | 7.0 |

| AVX2 | 45.2 |

| AVX2+OpenMP | 48.1 |

| CUDA | 520+ |

Key Takeaways

• The Mandelbrot set is a benchmark embarrassingly parallel task for testing SIMD, OpenMP, and CUDA
• AVX2 provides ~6.5x speedup over scalar code by processing 8 pixels per loop
• OpenMP adds 6-7% gain on 6 cores due to synchronization overhead
• CUDA delivers order-of-magnitude speedup thanks to thousands of parallel threads
• Performance measurements exclude rendering (hyperfine, N runs ~5s)

— Editorial Team