c++ - 使用 OpenMP 线程和 std::(experimental::)simd 计算 Mandelbrot 集

Question

我希望使用不同类型的 HPC 范例来实现一个简单的 Mandelbrot 集绘图器，展示它们的优点和缺点以及它们的实现的难易程度。想想 GPGPU (CUDA/OpenACC/OpenMP4.5)、线程/OpenMP 和 MPI。并使用这些示例为刚接触 HPC 的程序员提供帮助并了解可能性。代码的清晰度比从硬件中获得绝对的顶级性能更重要，这是第二步；)

因为这个问题对于并行化来说是微不足道的，而且现代 CPU 可以使用 vector 指令获得巨大的性能，所以我还想结合 OpenMP 和 SIMD。不幸的是，只需添加一个 #pragma omp simd不会产生令人满意的结果，并且使用内在函数不是非常用户友好或面向 future 。或者 pretty .

幸运的是，正在对 C++ 标准进行改进，这样应该更容易通用地实现 vector 指令，如 TS 中所述:"Extensions for parallelism, version 2" ，特别是关于数据并行类型的第 9 节。可以找到 WIP 实现 here ，基于VC可以找到here .

假设我有以下类(已对其进行更改以使其更简单一些)

#include <stddef.h>

using Range = std::pair<double, double>;
using Resolution = std::pair<std::size_t, std::size_t>;

class Mandelbrot
{
    double* d_iters;
    Range d_xrange;
    Range d_yrange;
    Resolution d_res;
    std::size_t d_maxIter;
    
public:
    Mandelbrot(Range xrange, Range yrange, Resolution res, std::size_t maxIter);
    ~Mandelbrot();

    void writeImage(std::string const& fileName);
    void computeMandelbrot();
private:
    void calculateColors();
}; 

以及以下 computeMandelbrot() 的实现使用 OpenMP

void Mandelbrot::computeMandelbrot()
{
    double dx = (d_xrange.second - d_xrange.first) / d_res.first;
    double dy = (d_yrange.second - d_yrange.first) / d_res.second;

    #pragma omp parallel for schedule(dynamic)
    for (std::size_t row = 0; row != d_res.second; ++row)
    {
        double c_imag = d_yrange.first + row * dy;
        for (std::size_t col = 0; col != d_res.first; ++col)
        {
            double real = 0.0;
            double imag = 0.0;
            double realSquared = 0.0;
            double imagSquared = 0.0;
            double c_real = d_xrange.first + col * dx;

            std::size_t iter = 0;
            while (iter < d_maxIter && realSquared + imagSquared < 4.0)
            {
                realSquared = real * real;
                imagSquared = imag * imag;
                imag = 2 * real * imag + c_imag;
                real = realSquared - imagSquared + c_real;
                ++iter;
            }
            d_iters[row * d_res.first + col] = iter;
        }   
    }
}

我们可以假设 x 和 y 方向的分辨率都是 2/4/8/.. 的倍数，具体取决于我们使用的 SIMD 指令。

不幸的是，关于std::experimental::simd 的在线信息很少。 .据我所知，也没有任何重要的例子。

在 Vc git 存储库中，有一个 Mandelbrot 集合计算器的实现，但它非常复杂，并且由于缺少注释而很难理解。

很明显，我应该在函数 computeMandelbrot() 中更改 double 的数据类型，但我不确定是什么。 TS 提到了一些类型 T 的两种主要新数据类型，

native_simd = std::experimental::simd<T, std::experimental::simd_abi::native>;

和

fixed_size_simd = std::experimental::simd<T, std::experimental::simd_abi::fixed_size<N>>;

使用 native_simd最有意义，因为我在编译时不知道我的界限。但是我不清楚这些类型代表什么，是 native_simd<double>单个 double 还是执行 vector 指令的 double 集合？那么这个系列中有多少 double ？

如果有人能给我指出使用这些概念的示例，或者给我一些关于如何使用 std::experimental::simd 实现 vector 指令的指示，我将不胜感激。

Answer 1

这是一个非常基本的实现，它可以工作(据我所知)。测试 vector 中哪些元素的绝对值大于 2 的方法非常繁琐且效率低下。一定有更好的方法来做到这一点，但我还没有找到。

我在 AMD Ryzen 5 3600 上获得了大约 72% 的性能提升，并为 g++ 提供了 -march=znver2 选项，这低于预期。

template <class T>
void mandelbrot(T xstart, T xend,
                T ystart, T yend)
{
    namespace stdx = std::experimental;

    constexpr auto simdSize = stdx::native_simd<T>().size();
    constexpr unsigned size = 4096;
    constexpr unsigned maxIter = 250;

    assert(size % simdSize == 0);
    unsigned* res = new unsigned[size * size];

    T dx = (xend - xstart) / size;
    T dy = (yend - ystart) / size;

    for (std::size_t row = 0; row != size; ++row)
    {
        T c_imag = ystart + row * dy;
        for (std::size_t col = 0; col != size; col += simdSize)
        {
            stdx::native_simd<T> real{0};
            stdx::native_simd<T> imag{0};
            stdx::native_simd<T> realSquared{0};
            stdx::native_simd<T> imagSquared{0};
            stdx::fixed_size_simd<unsigned, simdSize> iters{0};

            stdx::native_simd<T> c_real;
            for (int idx = 0; idx != simdSize; ++idx)
            {
                c_real[idx] = xstart + (col + idx) * dx;
            }

            for (unsigned iter = 0; iter != maxIter; ++iter)
            {
                realSquared = real * real;
                imagSquared = imag * imag;
                auto isInside = realSquared + imagSquared > stdx::native_simd<T>{4};
                for (int idx = 0; idx != simdSize; ++idx)
                {
                    // if not bigger than 4, increase iters
                    if (!isInside[idx])
                    {
                        iters[idx] += 1;
                    }
                    else
                    {
                        // prevent that they become inf/nan
                        real[idx] = static_cast<T>(4);
                        imag[idx] = static_cast<T>(4);
                    }
                }

                if (stdx::all_of(isInside) )
                {
                    break;
                }        

                imag = static_cast<T>(2.0) * real * imag + c_imag;
                real = realSquared - imagSquared + c_real;
            }
            iters.copy_to(res + row * size + col, stdx::element_aligned);
        }

    }
    delete[] res;
}

整个测试代码(从auto test = (...)开始)编译为

  .L9:
  vmulps ymm1, ymm1, ymm1
  vmulps ymm13, ymm2, ymm2
  xor eax, eax
  vaddps ymm2, ymm13, ymm1
  vcmpltps ymm2, ymm5, ymm2
  vmovaps YMMWORD PTR [rsp+160], ymm2
  jmp .L6
.L3:
  vmovss DWORD PTR [rsp+32+rax], xmm0
  vmovss DWORD PTR [rsp+64+rax], xmm0
  add rax, 4
  cmp rax, 32
  je .L22
.L6:
  vucomiss xmm3, DWORD PTR [rsp+160+rax]
  jp .L3
  jne .L3
  inc DWORD PTR [rsp+96+rax]
  add rax, 4
  cmp rax, 32
  jne .L6