含continue分支的循环向量化：Mandelbrot集计算性能优化求助

基础信息 #

• 硬件：搭载Intel(R) Xeon(R) Gold 6240 CPU的超算
• 编译器：icc，固定编译参数：-O3 -mavx2 -xHost -g -qopenmp-simd -qopt-report=1 -qopt-report-phase=vec
• 约束规则：仅允许向量化优化，禁止并行化；无法修改编译参数、图像宽高（宽高均为4096的倍数，内存对齐无问题）
• 已做优化：仅计算Mandelbrot集上半部分再复制到下半部分、为所有循环添加OpenMP SIMD指令、全循环向量化
• 性能差距：当前代码耗时约900ms，远慢于参考方案的500ms；Advixe Profiler分析显示内循环耗时主要集中在pointsRemaining、i2和r2变量

代码实现 #

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <cstring>
#include <stdlib.h>


#include "LineMandelCalculator.h"


LineMandelCalculator::LineMandelCalculator(unsigned matrixBaseSize, unsigned limit) :
    BaseMandelCalculator(matrixBaseSize, limit, "LineMandelCalculator")
{
    data = (int *)(_mm_malloc(height * width * sizeof(int), 64));
    zReal = (float *)(_mm_malloc(width * sizeof(float), 64));
    zImag = (float *)(_mm_malloc(width * sizeof(float), 64));
    iteration = (int *)(_mm_malloc(width * sizeof(int), 64));
    calculated_reals = (float *)(_mm_malloc(width * sizeof(float), 64));


    for (int i = 0; i < width; ++i) {
        calculated_reals[i] = x_start + i * dx;
        zReal[i] = 0.0f;
        zImag[i] = 0.0f;
        iteration[i] = limit;
    }
}

LineMandelCalculator::~LineMandelCalculator() {
    _mm_free(data);
    _mm_free(zReal);
    _mm_free(zImag);
    _mm_free(iteration);
    _mm_free(calculated_reals);

    data = NULL;
    zReal = NULL;
    zImag = NULL;
    iteration = NULL;
    calculated_reals = NULL;
}

int * LineMandelCalculator::calculateMandelbrot() {
    int *pdata = data;
    int half_height = height / 2;
    float *zReal = this -> zReal;
    float *zImag = this -> zImag;
    float *calculated_reals = this -> calculated_reals;
    int * iteration = this -> iteration;
    int limit = this -> limit;
    int width = this -> width;
    float dx = (float)this -> dx;
    float dy = (float)this -> dy;

    for (int i = 0; i < half_height; i++) {
        float y = y_start + i * dy;

        #pragma omp simd simdlen(64) aligned(zReal, zImag, iteration:64)
        for (int j = 0; j < width; j++) {
            zReal[j] = calculated_reals[j];
            zImag[j] = y;
            iteration[j] = limit;
        }

        int pointsRemaining = width; 

        for (int k = 0; k < limit && pointsRemaining!=0; ++k) {
            pointsRemaining = 0; 

            #pragma omp simd simdlen(64) reduction(+:pointsRemaining) aligned(zReal, zImag, calculated_reals, iteration:64)
            for (int j = 0; j < width; j++) {
                //if(iteration[j]!=limit){
                //    continue;
                //}
                float r2 = zReal[j] * zReal[j];
                float i2 = zImag[j] * zImag[j];
                bool inside = (r2 + i2 <= 4.0f);

                if (inside)pointsRemaining ++;

                zImag[j] = 2.0f * zReal[j] * zImag[j] + y;
                zReal[j] = r2 - i2 + calculated_reals[j];

                if (!inside && iteration[j] == limit){
                    iteration[j] = k;
                }
            }
        }

        memcpy(pdata + i * width, iteration, width * sizeof(int));

        memcpy(pdata + (height - i - 1) * width, pdata + i * width, width * sizeof(int));
    }

    return data;
}

代码逻辑与核心痛点 #

• 逻辑：对图像每行计算迭代次数（运行脚本指定，通常100-500次），只要行内存在未发散点，就对整行所有点执行迭代计算
• 痛点：曾尝试添加注释中的if分支跳过已收敛点，但导致代码变慢；若能实现该分支的向量化，预计性能可大幅提升

内容的提问来源于stack exchange，提问作者ampersander