使用OpenMP线程和std：:( experimentalal：:) simd计算Mandelbrot集

问题描述

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

由于该问题很难并行化，并且现代CPU可以使用矢量指令获得大量性能，因此我也想将OpenMP和SIMD结合起来。不幸的是，简单地添加#pragma omp simd并不能产生令人满意的结果，并且使用内在函数不是非常用户友好或未来的证明。或pretty。

幸运的是，正在按照C ++标准进行工作，以便更容易地通用实现矢量指令，如TS："Extensions for parallelism,version 2"中所述，尤其是关于数据并行类型的第9节。可以在here的基础上找到here的WIP实现。

假设我有以下课程（已更改为使其更简单一些）

#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>是单双精度还是是在其上执行了矢量指令的双精度集合？然后这个集合中有多少双打？

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

解决方法

这是一个非常基本的实现，有效（据我所知）。矢量的哪些元素的绝对值大于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,32
  jne .L6