Optimización SIMD de Mandelbrot AVX2 OpenMP CUDA

Optimización del cómputo del conjunto de Mandelbrot: desde código escalar hasta SIMD, OpenMP y CUDA

El conjunto de Mandelbrot se calcula de forma independiente para cada punto en el plano complejo: verificamos si la secuencia z_{n+1} = z_n² + c permanece acotada, empezando con z_0 = 0. Si |z_n| > 2, la secuencia diverge. Limitar las iteraciones a 256 garantiza que el cómputo finalice. La tarea es perfectamente paralelizable: cada píxel se procesa de forma independiente.

El color del píxel depende del número de iteraciones en el que |z| excede el radio 2. Los puntos que permanecen dentro de los límites tras 256 pasos se colorean de negro.

Implementación escalar

El algoritmo básico en C++ mapea las coordenadas de pantalla [0, WIDTH] × [0, HEIGHT] a la región del plano complejo [-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);

El bucle de iteración se optimiza almacenando z_x² y z_y² por separado:

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

En un AMD Ryzen 5 5600H a 1920×1080 con 256 iteraciones, el código escalar con -O2 ofrece 7.0 ± 0.1 FPS (GCC/Clang).

Vectorización con AVX2

SIMD en registros ymm de 256 bits procesa 8 flotantes por instrucción. Las intrínsecas de <x86intrin.h> proporcionan un control explícito sobre la vectorización.

Vectores clave:

• _c_x: coordenada base + offsets [0..7] × c_step_x
• _z_x, _z_y, _z_x2, _z_y2, _z_xy: estados para 8 píxeles
• _iterations: contadores de iteraciones (int32)

Inicialización de offset:

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

En el bucle de iteración:

• Calcular radio: _radius2 = _mm256_add_ps(_z_x2, _z_y2)
• Máscara de continuación: _mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OS) → todos -1.0 o 0.0
• Actualizar estados bajo la máscara activa
• Desplazar a lo largo de X tras procesar 8 píxeles: _c_x = _mm256_add_ps(_c_x, _8_c_steps_x)

La máscara permite salir del bucle para los 8 píxeles simultáneamente si todos han escapado.

Multihilo con OpenMP

OpenMP distribuye las filas entre hilos mediante la directiva #pragma omp parallel for. Cada hilo calcula de forma independiente su tira de píxeles.

#pragma omp parallel for schedule(dynamic)

`schedule(dynamic) para equilibrado de carga: diferentes regiones del fractal requieren un número variable de iteraciones.

Combinando AVX2 + OpenMP en el Ryzen 5 5600H de 6 núcleos se obtienen hasta 40-50 FPS.

Aceleración con GPU mediante CUDA

La adaptación a CUDA es directa: cada hilo calcula un píxel. Kernel de 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;
}

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

En una GPU de portátil (RTX 3050), alcanza más de 500 FPS a Full HD.

| Implementación | FPS (Ryzen 5 5600H) |

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

| Escalar | 7.0 |

| AVX2 | 45.2 |

| AVX2+OpenMP | 48.1 |

| CUDA | 520+ |

Lecciones clave

• El conjunto de Mandelbrot es una tarea de referencia embarazosamente paralela para probar SIMD, OpenMP y CUDA
• AVX2 proporciona una aceleración de ~6.5x sobre el código escalar al procesar 8 píxeles por bucle
• OpenMP añade un 6-7% de mejora en 6 núcleos debido a la sobrecarga de sincronización
• CUDA ofrece una aceleración de orden de magnitud gracias a miles de hilos paralelos
• Las mediciones de rendimiento excluyen el renderizado (hyperfine, N ejecuciones ~5s)

— Editorial Team