OpenMPスレッドとstd::( experimental::) simdを使用してマンデルブロ集合を計算する

Aug 19 2020

さまざまな種類のHPCパラダイムを使用して、単純なマンデルブロ集合プロッターを実装することを検討しています。その長所と短所、および実装の容易さや難しさを示しています。GPGPU(CUDA / OpenACC / OpenMP4.5)、スレッド化/ OpenMPおよびMPIについて考えてみてください。そして、これらの例を使用して、HPCを初めて使用するプログラマーに手がかりを与え、可能性が何であるかを確認します。コードの明確さは、ハードウェアから絶対的な最高のパフォーマンスを引き出すことよりも重要です。それが2番目のステップです;)

問題は並列化するのは簡単であり、最近のCPUはベクトル命令を使用して膨大な量のパフォーマンスを得ることができるため、OpenMPとSIMDも組み合わせたいと思います。残念ながら、単にaを追加し#pragma omp simdても満足のいく結果は得られず、組み込み関数を使用することはあまりユーザーフレンドリーではなく、将来も保証されません。またはかなり。

幸いなことに、TS:「並列処理の拡張バージョン2」、特にデータ並列型に関するセクション9で説明されているように、ベクトル命令を一般的に実装するのが簡単になるように、C ++標準の作業が行われています。WIPの実装はここにあります。これはVCに基づいており、ここにあります。

私が次のクラスを持っていると仮定します(少し簡単にするために変更されました)

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

使用するSIMD命令に応じて、x方向とy方向の両方の解像度が2/4/8 / ..の倍数であると想定できます。

残念ながら、オンラインで入手できる情報はほとんどありませんstd::experimental::simd。私が見つけることができる限り、重要な例もありません。

Vc gitリポジトリには、マンデルブロ集合計算機の実装がありますが、それは非常に複雑であり、コメントが不足しているため、フォローするのはかなり困難です。

関数内のdoubleのデータ型を変更する必要があることは明らかcomputeMandelbrot()ですが、何がわからないのです。TSは、いくつかのタイプTの2つの主要な新しいデータタイプについて言及しています。

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を使用してベクトル命令を実装する方法についてのポインタを教えてくれたら、とてもありがたいです。

回答

NigelOvermars Aug 24 2020 at 19:45

これは非常に基本的な実装であり、(私が知る限り)機能します。ベクトルのどの要素の絶対値が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