矢量化带有“continue”分支的循环

问：

我目前正在接受一门课程，在该课程中，我们可以使用超级计算机。CPU 是 Intel（R） Xeon（R） Gold 6240 CPU。我们的任务是矢量化（但不允许使用并行化）曼德布洛特集的计算。目标是使代码尽可能快。我尝试了一些优化。仅计算曼德布洛特集的上半部分，然后向下复制值，利用所有可能的 OpenMP 编译指示，并尝试对每个循环进行矢量化。然而，我的代码比它应该的要慢得多（~900ms，而参考解决方案在~500ms）。我需要帮助使此代码运行得更快。我们使用这些标志在编译器上编译。我运行了一个Advixe Profiler，并在内部循环上收到了以下时间。但是，在分析器的多次运行中，每行花费的时间差异很大，尽管大多数时间总是花在 and or 变量上：icc-O3 -mavx2 -xHost -g -qopenmp-simd -qopt-report=1 -qopt-report-phase=vecpointsRemainingi2r2

这是我的代码（它将从自定义脚本运行 - 我无法更改编译，也无法更改集合的宽度或高度等任何内容 - 但这两者将是 4096 的倍数，所有对齐都可以）。对于图像中的每一行，它都会计算迭代次数（其数量由运行计算器的脚本给出，通常为 ~100-500）。只要至少有一个点没有偏离曼德布洛特集，程序就会计算一行中所有点的迭代次数（我尝试跳过每个点的计算，但添加注释掉的“if”子句只会减慢代码速度 - i我认为如果我设法矢量化注释的“if”， 它会运行得更快）：

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