Mandelbrot 的 SIMD 优化 AVX2 OpenMP CUDA

优化曼德勃罗集计算：从标量代码到 SIMD、OpenMP 和 CUDA

曼德勃罗集的每个复平面点都可以独立计算：从 z_0 = 0 开始，检查序列 z_{n+1} = z_n² + c 是否保持有界。如果 |z_n| > 2，则序列发散。将迭代上限设为 256 可确保计算完成。该任务完美并行化——每个像素独立处理。

像素颜色取决于 |z| 超过半径 2 时的迭代次数。经过 256 步仍保持在界内的点着黑色。

标量实现

基本的 C++ 算法将屏幕坐标 [0, WIDTH] × [0, HEIGHT] 映射到复平面区域 [-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);

迭代循环通过分别存储 z_x² 和 z_y² 来优化：

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

在 AMD Ryzen 5 5600H 上，以 1920×1080 分辨率和 256 次迭代，标量代码使用 -O2 优化可实现 7.0 ± 0.1 帧/秒（GCC/Clang）。

AVX2 向量化

256 位 ymm 寄存器的 SIMD 每条指令可处理 8 个 float。<x86intrin.h> 中的 intrinsics 提供对向量化的明确控制。

关键向量：

• _c_x：基础坐标 + 偏移 [0..7] × c_step_x
• _z_x、_z_y、_z_x2、_z_y2、_z_xy：8 个像素的状态
• _iterations：迭代计数器（int32）

偏移初始化：

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

在迭代循环中：

• 计算半径：_radius2 = _mm256_add_ps(_z_x2, _z_y2)
• 继续掩码：_mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OS) → 全为 -1.0 或 0.0
• 在活动掩码下更新状态
• 处理 8 个像素后沿 X 方向移位：_c_x = _mm256_add_ps(_c_x, _8_c_steps_x)

掩码允许如果所有 8 个像素均已逃逸，则同时中断循环。

OpenMP 多线程

OpenMP 使用 #pragma omp parallel for 指令将行分布到各线程。每个线程独立计算其像素条带。

#pragma omp parallel for schedule(dynamic)

schedule(dynamic) 用于负载均衡——不同分形区域需要不同迭代次数。

在 6 核 Ryzen 5 5600H 上结合 AVX2 + OpenMP 可达 40-50 帧/秒。

使用 CUDA 的 GPU 加速

移植到 CUDA 非常简单：每个线程计算一个像素。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;
}

启动：dim3 block(16,16); dim3 grid((width+15)/16, (height+15)/16); mandelbrot_cuda<<<grid,block>>>(...)

在笔记本 GPU (RTX 3050) 上，全高清分辨率下可达 500+ 帧/秒。

| 实现 | 帧/秒 (Ryzen 5 5600H) |

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

| 标量 | 7.0 |

| AVX2 | 45.2 |

| AVX2+OpenMP | 48.1 |

| CUDA | 520+ |

关键要点

• 曼德勃罗集是测试 SIMD、OpenMP 和 CUDA 的尴尬并行基准任务
• AVX2 通过每循环处理 8 个像素，比标量代码提速约 6.5 倍
• OpenMP 在 6 核上额外提升 6-7%，因同步开销所限
• CUDA 借助数千并行线程，实现数量级提速
• 性能测量排除渲染（hyperfine，N 次运行约 5s）

— Editorial Team