向量化一个带有“继续”分支的循环

Question

我目前正在学习一门课程，在课程中我们可以使用超级计算机。 CPU 是 Intel(R) Xeon(R) Gold 6240 CPU。我们的任务是向量化（但不允许使用并行化）曼德尔布罗集的计算。目标是使代码尽可能快。我尝试了很多优化。仅计算 mandelbrot 集的上半部分，然后向下复制值，利用所有可能的 OpenMP 编译指示，并尝试对每个循环进行矢量化。然而，我的代码比预期的要慢得多（约 900 毫秒，而参考解决方案约为 500 毫秒）。我需要帮助使该代码运行得更快。我们使用这些标志

icc

在

-O3 -mavx2 -xHost -g -qopenmp-simd -qopt-report=1 -qopt-report-phase=vec

编译器上进行编译。我运行了 Advixe Profiler，并在内循环上收到了以下时间。然而，在分析器的多次运行中，每行花费的时间差异很大，尽管大多数时间总是花费在

pointsRemaining

和

i2

或

r2

变量上：

这是我的代码（它将从自定义脚本运行 - 我无法更改编译，也无法更改集合的宽度或高度之类的任何内容 - 但这两个将是 4096 的倍数，所有对齐都可以）。对于图像中的每一行，它都会计算迭代次数（迭代次数由运行计算器的脚本给出，通常约为 100-500）。只要至少有一个点没有偏离曼德尔布罗特集，程序就会计算一行中所有点的迭代（我尝试跳过对每个点的计算，但添加注释掉的“if”子句只会减慢速度写下代码 - 我认为如果我设法对注释的“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;
}

Answer 1

这个问题肯定是由于

!inside && iteration[j] == limit

造成的，需要在 C 和 C++ 中进行 惰性求值 。编译器无法对此类条件进行向量化，因为获取

iteration[j]

具有可见行为：编译器无法评估

iteration[j]

，因为根据

j

的值，

iteration

可能超出

!inside

的范围。因此，ICC（和其他编译器）无法矢量化整个循环并生成条件跳转的缓慢级联。

如果您知道可以立即评估这样的条件，那么您可以将其替换为

(!inside) & (iteration[j] == limit)

。使用二元

&

运算符的表达式不会被延迟计算。这使得 ICC 能够生成完全矢量化的代码。

话虽这么说，生成的代码很大（大约 12 KiB），其中热循环有 1.5 KiB。太大的代码可能会降低性能。事实上，现代 x86-64 处理器将指令转换为微指令 (uop)。翻译通常可能成为此类热循环的瓶颈。因此，它们有一个微指令缓存，不能重新翻译指令，但它的大小有限。太大的循环会导致微指令缓存丢失，并且指令会一遍又一遍地重新解码，速度很慢。关于所使用的确切目标 CPU，这可能是一个问题，也可能不是一个问题。在 Intel(R) Xeon(R) Gold 6240（Cascade Lake CPU）上，我预计不会成为重大瓶颈。请注意，指令的对齐也很重要，尤其是在 Cascade Lake 上（请参阅这篇文章）。乍一看，

simdlen(64)

就是生成的代码如此之大的原因。这可能不合理。至少，对于 AVX-2 来说肯定不是，其中

simdlen(32)

肯定是首选。

您指定冲突的编译器标志。事实上，

-xHost

告诉编译器使用目标机器的可用指令集（AFAIK 目标平台上的 AVX-512），而

-mavx2

告诉编译器可以使用 AVX-2。因此，

-mavx2

对于

-xHost

来说毫无用处。指定是否要对编译器使用 AVX-2 或 AVX-512，而不是同时使用两者。您的目标 CPU 似乎支持 2 个 AVX-512 单元，因此 AVX-512 是最佳解决方案。

矢量化后，代码还不错，但效率也不是很高。一个非常好的 Mandelbrot 实现应该将大部分时间花在执行 FMA 指令上。这里的情况似乎并非如此。不过，看起来 AVX-512 ALU 是瓶颈，这很好。