SIMD Mandelbrot计算

Question

前几天我在玩GPU电脑着色器，创建了一个Mandelbrot着色器。不幸的是，金属不支持双精度的计算机着色器，所以超过一定的缩放水平，我需要切换回CPU。在这样做时，我决定尝试为计算编写SIMD代码，以使其更快。

在下面的代码中，我使用了AVX512指令，并且在标量代码上得到了一个加速比。我将图像分解成64x64像素的瓷砖，并将它们输出到可用的核心。对于一个特定测试图像上的标量代码，计算一个瓷砖的平均时间是0.757288秒。下面的SIMD版本是0.466437。这大约是33%的增长，这是可以的。考虑到我一次计算8倍的像素，我希望得到更多。

这些是我在代码中使用的一些有用的类型。

#include <immintrin.h>
typedef struct RGBA8Pixel {
    uint8_t red;
    uint8_t green;
    uint8_t blue;
    uint8_t alpha;
} RGBA8Pixel;

typedef union intVec8 {
    __m512i ivec;
    int64_t vec[8];
} intVec8;

typedef union doubleVec8 {
    __m512d dvec;
    double vec[8];
} doubleVec8;

这是我计算1 64x64瓷砖的函数：

- (void)calculateSIMDFromRow:(int)startPixelRow
                       toRow:(int)endPixelRow
                     fromCol:(int)startPixelCol
                       toCol:(int)endPixelCol;
{
    if (!_keepRendering)
    {
        return;
    }

    const doubleVec8 k0s = {
        .vec[0] = 0.0,
        .vec[1] = 0.0,
        .vec[2] = 0.0,
        .vec[3] = 0.0,
        .vec[4] = 0.0,
        .vec[5] = 0.0,
        .vec[6] = 0.0,
        .vec[7] = 0.0,
    };

    const intVec8   k1s = {
        .vec[0] = 1,
        .vec[1] = 1,
        .vec[2] = 1,
        .vec[3] = 1,
        .vec[4] = 1,
        .vec[5] = 1,
        .vec[6] = 1,
        .vec[7] = 1,
    };

    const doubleVec8    k2s = {
        .vec[0] = 2.0,
        .vec[1] = 2.0,
        .vec[2] = 2.0,
        .vec[3] = 2.0,
        .vec[4] = 2.0,
        .vec[5] = 2.0,
        .vec[6] = 2.0,
        .vec[7] = 2.0,
    };

    const doubleVec8    k4s = {
        .vec[0] = 4.0,
        .vec[1] = 4.0,
        .vec[2] = 4.0,
        .vec[3] = 4.0,
        .vec[4] = 4.0,
        .vec[5] = 4.0,
        .vec[6] = 4.0,
        .vec[7] = 4.0,
    };

    UInt64      maxIterations   = [self maxIterations];
    NSSize      viewportSize    = [self viewportSize];
    for (int row = startPixelRow; (row < endPixelRow) && (_keepRendering); ++row)
    {
        RGBA8Pixel* nextPixel   = _outputBitmap + (row * (int)viewportSize.width) + startPixelCol;
        double      yCoord      = _yCoords [ row ];
        doubleVec8  yCoords;
        for (int i = 0; i < 8; i++)
        {
            yCoords.vec [ i ] = yCoord;
        }
        double*     nextXCoord  = &_xCoords [ startPixelCol ];
        for (int col = startPixelCol; (col < endPixelCol) && (_keepRendering); col += 8)
        {
            __m512d as = _mm512_load_pd(nextXCoord);
            nextXCoord += 8;
            __m512d bs = yCoords.dvec;
            __m512d cs = as;
            __m512d ds = bs;

            UInt64 scalarIters = 1;
            __m512i iterations  = k1s.ivec;
            __m512d dists       = k0s.dvec;
            __mmask8 allDone    = 0;
            while ((allDone != 0xFF) && (_keepRendering))
            {
                // newA = a * a - b * b + c
                __m512d newA;
                __m512d newB;
                newA = _mm512_mul_pd(as, as);
                newA = _mm512_sub_pd(newA, _mm512_mul_pd(bs, bs));
                newA = _mm512_add_pd(newA, cs);

                //double    newB    = 2 * a * b + d;
                newB = _mm512_mul_pd(_mm512_mul_pd(k2s.dvec, as), bs);
                newB = _mm512_add_pd(newB, ds);

                as = newA;
                bs = newB;

                dists = _mm512_mul_pd(newB, newB);
                dists = _mm512_add_pd(_mm512_mul_pd(newA, newA), dists);
                __mmask8 escaped = _mm512_cmplt_pd_mask(dists, k4s.dvec);

                iterations = _mm512_mask_add_epi64(iterations, escaped, iterations, k1s.ivec);
                scalarIters++;
                __mmask8 hitMaxIterations = (scalarIters == maxIterations) ? 0xFF : 0;

                allDone = ~escaped | hitMaxIterations;
            }

            intVec8 iters = { .ivec = iterations };
            for (int i = 0; i < 8; i++)
            {
                UInt64  nextIteration = iters.vec [ i ];
                if (nextIteration == maxIterations)
                {
                    *nextPixel = kBlack;
                }
                else
                {
                    *nextPixel = kPalette [ nextIteration % kPaletteSize ];
                }
                nextPixel++;
            }
        }
    }
}

我是英特尔SIMD指令的新手，坦率地说，我发现它们非常令人困惑。如果有更好的方法来做上述任何一件事，请告诉我。我试着使用融合的乘法-加法和乘法加-负指令，不幸的是，他们使代码比使用2或3个单独的指令慢得多。

我正在使用Xcode 10.2.1开发macOS，使用在<immintrin.h>中找到的Intel数据类型和本质。

Answer 1

有几个想法：

1. 你说，考虑到我一次计算8倍的像素，我希望得到更多。是的，simd在做向量/矩阵运算时提供了一些非常出色的性能改进，但是对于Mandelbrot，您所做的只是元素相加和乘法，所以您应该看到改进，但是对于这些简单的元素级计算，没有接近8×。在我的测试中，simd的性能几乎是i9 Mac和iPhone Xs标量渲染性能的两倍。
2. 两个算法观察：
   • 我注意到您对as和bs进行了两次平方，一次是在算法期间，另一次是在转义测试中。我建议对此进行重构，这样您就可以在算法和转义测试中使用这个平方运算的结果。
   • 我注意到您正在为2做2×a×b部分的向量乘法。我使用的是向量×标量积，而不是向量元素积。那可能会更快一些。

3. 如果您想要消除那些不直观的_mm512_xxx调用，可以考虑simd库，它是加速框架的一部分。这是一个更高层次的抽象，特别是在Swift中，您最终得到了非常自然的代码，这是IMHO，更容易阅读。
4. 在概念上，应该注意到simd性能的提高将被一些像素已经转义的边界情况所抵消，而其他的则没有，因为您将计算向量中所有八个像素的迭代，包括那些已经转义的像素。这可能对性能没有实质性的影响，但值得注意，尤其是在处理Mandelbrot集中复杂的部分(这是最有趣的部分)的情况下。

不管它的价值是什么，这就是Swift simd版本看起来的样子：

import simd

func calculate(real: simd_double8, imaginary: simd_double8) -> simd_double8 {
    var zReal = real // simd_double8.zero
    var zImaginary = imaginary // simd_double8.zero

    let thresholds = simd_double8(repeating: 4)
    let maxIterations = 10_000.0

    var notEscaped = SIMDMask<SIMD8<Double.SIMDMaskScalar>>(repeating: true)
    let isDone = SIMDMask<SIMD8<Double.SIMDMaskScalar>>(repeating: false)

    var currentIterations = 0.0
    var iterations = simd_double8.zero

    repeat {                                                    // z = z^2 + c
        currentIterations += 1.0
        iterations.replace(with: currentIterations, where: notEscaped)

        let zRealSquared = zReal * zReal
        let zImaginarySquared = zImaginary * zImaginary

        zImaginary = 2.0 * zReal * zImaginary + imaginary       // 2 × zr × zi + ci
        zReal = zRealSquared - zImaginarySquared + real         // zr^2 - zi^2 + cr

        notEscaped = zRealSquared + zImaginarySquared .< thresholds
    } while notEscaped != isDone && currentIterations < maxIterations

    iterations.replace(with: 0, where: notEscaped)

    return iterations
}

值得注意的是，它非常类似于标量呈现，没有神秘的方法引用。例如，这里是标量版本：

func calculate(real: Double, imaginary: Double) -> Int {
    var zReal = real
    var zImaginary = imaginary

    let thresholds = 4.0
    let maxIterations = 10_000

    var notEscaped = false

    var currentIterations = 0

    repeat {                                                    // z = z^2 + c
        currentIterations += 1

        let zRealSquared = zReal * zReal
        let zImaginarySquared = zImaginary * zImaginary

        zImaginary = 2.0 * zReal * zImaginary + imaginary       // 2 × zr × zi + ci
        zReal = zRealSquared - zImaginarySquared + real         // zr^2 - zi^2 + cr

        notEscaped = zRealSquared + zImaginarySquared < thresholds
    } while notEscaped && currentIterations < maxIterations

    return currentIterations >= maxIterations ? 0 : currentIterations
}

Answer 2

我试着使用融合的乘法-加法和乘法加-负指令，不幸的是，他们使代码比使用2或3个单独的指令慢得多。

这通常表示延迟瓶颈，

这大约是33%的增长，这是可以的。考虑到我一次计算8倍的像素，我希望得到更多。

这个也是。

实际上，这对于Mandelbrot来说是预期的，因为它基于迭代函数应用程序，因此它本身就有一个非平凡的循环携带依赖项。Intel上的浮点操作具有较高的吞吐量，但与其吞吐量相比，它们仍然是缓慢的操作，具有较高的延迟时间。在Skylake上(我猜您使用的是AVX512)，一个FMA需要4个周期，但是处理器可以每个周期启动两个周期。因此，如果它们太“捆绑”(而且FMA的事情变得更加紧密，因为每个FMA都在等待3而不是2个输入准备好)，可能总是会执行一些浮点操作，但实际上在Skylake X上，我们希望8个操作在任何时候都“繁忙”。在哈斯韦尔，这是更糟的，采取10个“重叠”的FMA饱和浮点单位。

这种情况可以通过交错计算几个(4? 8?)独立的8像素行，虽然这的一个不幸的副作用是，它也会“将”循环计数到最大计数块中的所有像素。这种情况现在已经在较小的范围内发生了，但情况会变得更糟，并抑制了这样做的潜在收益。