Utilisation des threads OpenMP et std: :( ​​experimental: :) simd pour calculer l'ensemble de Mandelbrot

Aug 19 2020

Je cherche à implémenter un traceur de jeu Mandelbrot simple en utilisant différents types de paradigmes HPC, montrant leurs forces et leurs faiblesses et la facilité ou la difficulté de leurs implémentations. Pensez à GPGPU (CUDA / OpenACC / OpenMP4.5), threading / OpenMP et MPI. Et utilisez ces exemples pour donner aux programmeurs novices une prise en main et pour voir quelles sont les possibilités. La clarté du code est plus importante que l'obtention des meilleures performances absolues du matériel, c'est la deuxième étape;)

Comme le problème est simple à paralléliser et que les processeurs modernes peuvent gagner énormément de performances en utilisant des instructions vectorielles, je souhaite également combiner OpenMP et SIMD. Malheureusement, le simple fait d'ajouter a #pragma omp simdne donne pas de résultats satisfaisants et l'utilisation des intrinsèques n'est pas très conviviale ou à l'épreuve du temps. Ou jolie .

Heureusement, des travaux sont en cours sur le standard C ++ de sorte qu'il devrait être plus facile d'implémenter de manière générique des instructions vectorielles, comme mentionné dans le TS: "Extensions for parallelism, version 2" , en particulier la section 9 sur les types de parallélisme de données. Une implémentation WIP peut être trouvée ici , qui est basée sur VC qui peut être trouvée ici .

Supposons que j'ai la classe suivante (qui a été modifiée pour la rendre un peu plus simple)

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

Et l'implémentation suivante de l' computeMandelbrot()utilisation d'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;
        }   
    }
}

Nous pouvons supposer que les résolutions dans les directions x et y sont des multiples de 2/4/8 / .., en fonction des instructions SIMD que nous utilisons.

Malheureusement, il y a très peu d'informations disponibles en ligne sur std::experimental::simd. Ni d'exemples non triviaux pour autant que je puisse trouver.

Dans le dépôt Vc git, il y a une implémentation du calculateur d'ensemble de Mandelbrot, mais c'est assez compliqué et en raison du manque de commentaires plutôt difficile à suivre.

Il est clair que je devrais changer les types de données des doubles dans la fonction computeMandelbrot(), mais je ne sais pas quoi. Le TS mentionne deux nouveaux types de données principaux pour certains types T,

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

et

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

L'utilisation native_simdest la plus logique, car je ne connais pas mes limites au moment de la compilation. Mais alors, je ne sais pas ce que représentent ces types, est-ce native_simd<double>un simple double ou est-ce une collection de doubles sur laquelle une instruction vectorielle est exécutée? Et puis combien de doubles y a-t-il dans cette collection?

Si quelqu'un pouvait m'indiquer des exemples où ces concepts sont utilisés, ou me donner des conseils sur la façon d'implémenter des instructions vectorielles en utilisant std :: experimental :: simd, je serais très reconnaissant.

Réponses

NigelOvermars Aug 24 2020 at 19:45

Voici une implémentation très basique, qui fonctionne (pour autant que je sache). Le test des éléments du vecteur ayant une valeur absolue supérieure à 2 est effectué de manière très lourde et inefficace. Il doit y avoir une meilleure façon de faire cela, mais je ne l'ai pas encore trouvée.

J'obtiens une augmentation des performances d'environ 72% sur un AMD Ryzen 5 3600 et en donnant l'option g ++ -march=znver2, ce qui est moins que prévu.

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

L'ensemble du code de test (à partir de auto test = (...)) se compile jusqu'à

  .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