Intuition sur la disposition de la mémoire pour une conception rapide orientée SIMD / données

Une des choses qui n'a pas encore été mentionnée est que l'accès à la mémoire a un peu de latence. Et bien sûr, lors de la lecture à partir de 4 pointeurs, le résultat est disponible lorsque la dernière valeur arrive. Ainsi, même si 3 valeurs sur 4 sont dans le cache, la dernière valeur devra peut-être provenir de la mémoire et bloquer toute votre opération.

C'est pourquoi SSE ne prend même pas en charge ce mode. Toutes vos valeurs doivent être contiguës en mémoire, et pendant un certain temps, elles ont dû être alignées (afin qu'elles ne puissent pas franchir une limite de ligne de cache).

Fait important, cela signifie que votre exemple (Structure of Arrays) ne fonctionne pas dans le matériel SSE . Vous ne pouvez pas utiliser l'élément [1] de 4 vecteurs différents en une seule opération. Vous pouvez utiliser des éléments [0] à [3] à partir d'un seul vecteur.

J'ai implémenté un benchmark simple pour les deux méthodes.

Résultat: La mise en page par bandes est au mieux 10% plus rapide que la mise en page standard *. Mais avec SSE4.1, nous pouvons faire beaucoup mieux.

* Lorsque compilé avec gcc -Ofast sur un processeur i5-7200U .

La structure est em> un peu plus facile à travailler, mais beaucoup moins polyvalent. Cela pourrait cependant avoir un petit avantage dans un scénario réel, une fois que l'allocateur est suffisamment occupé.

Disposition par bandes

Durée 4624 ms

/* -----------------------------------------------------------------------------
        AoS layout [{X,y,Z,w},{X,y,Z,w},{X,y,Z,w},{X,y,Z,w},...]
----------------------------------------------------------------------------- */

using AoS_scene = std::vector<__m128>;

void print_scene(AoS_scene const &scene)
{
        // This is likely undefined behavior. Data might need to be stored
        // differently, but this is simpler to index.
        auto &&punned_data = reinterpret_cast<float const *>(scene.data());
        auto scene_size = std::size(scene);

        // Limit to 8 lines
        for(size_t j = 0lu; j < std::min(scene_size, 8lu); ++j) {
                for(size_t i = 0lu; i < 4lu; ++i) {
                        printf("%10.3e ", punned_data[j * 4lu + i]);
                }
                printf("\n");
        }
        if(scene_size > 8lu) {
                printf("(%lu more)...\n", scene_size - 8lu);
        }
        printf("\n");
}

void normalize(AoS_scene &scene)
{
        // Euclidean norm, SIMD 4 x 4D-vectors at a time.
        for(size_t i = 0lu; i < scene.size(); i += 4lu) {
                __m128 vec = scene[i];
                __m128 dot = _mm_dp_ps(vec, vec, 255);
                __m128 norms = _mm_sqrt_ps(dot);
                scene[i] = _mm_div_ps(vec, norms);
        }
}

float randf()
{
        std::random_device random_device;
        std::default_random_engine random_engine{random_device()};
        std::uniform_real_distribution<float> distribution(-10.0f, 10.0f);
        return distribution(random_engine);
}

int main()
{
        // Scene description, e.g. cameras, or particles, or boids etc.
        std::vector<__m128> scene(40'000);

        for(size_t i = 0lu; i < std::size(scene); ++i) {
                scene[i] = _mm_set_ps(randf(), randf(), randf(), randf());
        }

        // Print, normalize 100'000 times, print again

        // Compiler is hopefully not smart enough to realize
        // idempotence of normalization
        using std::chrono::steady_clock;
        using std::chrono::duration_cast;
        using std::chrono::milliseconds;
        // >:(

        print_scene(scene);
        printf("Working...\n");

        auto begin = steady_clock::now();
        for(int j = 0; j < 100'000; ++j) {
                normalize(scene);
                //break;
        }
        auto end = steady_clock::now();
        auto duration = duration_cast<milliseconds>(end - begin);

        printf("Time %lu ms\n", duration.count());
        print_scene(scene);

        return 0;
}

Memory usage summary: heap total: 746552, heap peak: 713736, stack peak: 2720
         total calls   total memory   failed calls
 malloc|          5         746552              0
realloc|          0              0              0  (nomove:0, dec:0, free:0)
 calloc|          0              0              0
   free|          2         672816
Histogram for block sizes:
    0-15              1  20% =========================
 1024-1039            1  20% =========================
32816-32831           1  20% =========================
   large              2  40% ==================================================

Disposition SoA

Durée 4982 ms

#include <chrono>
#include <cstdio>
#include <random>
#include <vector>
#include <xmmintrin.h>

/* -----------------------------------------------------------------------------
        SoA layout [X,X,X,X,...], [y,y,y,y,...], [Z,Z,Z,Z,...], ...
----------------------------------------------------------------------------- */

struct SoA_scene {
        size_t size;
        float *xs;
        float *ys;
        float *zs;
        float *ws;
};

void print_scene(SoA_scene const &scene)
{
        // This is likely undefined behavior. Data might need to be stored
        // differently, but this is simpler to index.

        // Limit to 8 lines
        for(size_t j = 0lu; j < std::min(scene.size, 8lu); ++j) {
                printf("%10.3e ", scene.xs[j]);
                printf("%10.3e ", scene.ys[j]);
                printf("%10.3e ", scene.zs[j]);
                printf("%10.3e ", scene.ws[j]);
                printf("\n");
        }
        if(scene.size > 8lu) {
                printf("(%lu more)...\n", scene.size - 8lu);
        }
        printf("\n");
}

void normalize(SoA_scene &scene)
{
        // Euclidean norm, SIMD 4 x 4D-vectors at a time.
        for(size_t i = 0lu; i < scene.size; i += 4lu) {
                __m128 xs = _mm_load_ps(&scene.xs[i]);
                __m128 ys = _mm_load_ps(&scene.ys[i]);
                __m128 zs = _mm_load_ps(&scene.zs[i]);
                __m128 ws = _mm_load_ps(&scene.ws[i]);

                __m128 xxs = _mm_mul_ps(xs, xs);
                __m128 yys = _mm_mul_ps(ys, ys);
                __m128 zzs = _mm_mul_ps(zs, zs);
                __m128 wws = _mm_mul_ps(ws, ws);

                __m128 xx_yys = _mm_add_ps(xxs, yys);
                __m128 zz_wws = _mm_add_ps(zzs, wws);

                __m128 xx_yy_zz_wws = _mm_add_ps(xx_yys, zz_wws);

                __m128 norms = _mm_sqrt_ps(xx_yy_zz_wws);

                __m128 normed_xs = _mm_div_ps(xs, norms);
                __m128 normed_ys = _mm_div_ps(ys, norms);
                __m128 normed_zs = _mm_div_ps(zs, norms);
                __m128 normed_ws = _mm_div_ps(ws, norms);

                _mm_store_ps(&scene.xs[i], normed_xs);
                _mm_store_ps(&scene.ys[i], normed_ys);
                _mm_store_ps(&scene.zs[i], normed_zs);
                _mm_store_ps(&scene.ws[i], normed_ws);
        }
}

float randf()
{
        std::random_device random_device;
        std::default_random_engine random_engine{random_device()};
        std::uniform_real_distribution<float> distribution(-10.0f, 10.0f);
        return distribution(random_engine);
}

int main()
{
        // Scene description, e.g. cameras, or particles, or boids etc.
        // Has to be a multiple of 4!   -- No edge case handling.
        auto scene_size = 40'000lu;
        std::vector<float> xs(scene_size);
        std::vector<float> ys(scene_size);
        std::vector<float> zs(scene_size);
        std::vector<float> ws(scene_size);

        for(size_t i = 0lu; i < scene_size; ++i) {
                xs[i] = randf();
                ys[i] = randf();
                zs[i] = randf();
                ws[i] = randf();
        }

        SoA_scene scene{
                scene_size,
                std::data(xs),
                std::data(ys),
                std::data(zs),
                std::data(ws)
        };
        // Print, normalize 100'000 times, print again

        // Compiler is hopefully not smart enough to realize
        // idempotence of normalization
        using std::chrono::steady_clock;
        using std::chrono::duration_cast;
        using std::chrono::milliseconds;
        // >:(

        print_scene(scene);
        printf("Working...\n");

        auto begin = steady_clock::now();
        for(int j = 0; j < 100'000; ++j) {
                normalize(scene);
        }
        auto end = steady_clock::now();
        auto duration = duration_cast<milliseconds>(end - begin);

        printf("Time %lu ms\n", duration.count());
        print_scene(scene);

        return 0;
}

Un inconvénient majeur de l'entrelacement sur la largeur du vecteur est que vous devez modifier la disposition pour tirer parti de vecteurs plus larges. (AVX, AVX512).

Mais oui, lorsque vous vectorisez uniquement manuellement (sans boucles qu'un compilateur pourrait vectoriser automatiquement avec son choix de largeur de vecteur), cela pourrait peut-être en valoir la peine si tout vos boucles (importantes) utilisent toujours tous les membres de la structure.

Sinon, le point de Max s'applique: une boucle qui ne touche que x et y gaspillera bande passante sur les membres z et w .

Ce ne sera pas moyen em> plus rapide, cependant ; avec une quantité raisonnable de déroulement de boucle, l'indexation de 4 tableaux ou l'incrémentation de 4 pointeurs est à peine pire que 1. La prélecture matérielle sur les processeurs Intel peut suivre un flux avant + 1 flux arrière par page 4k, donc 4 flux d'entrée est fondamentalement bien.

(Mais L2 est associatif à 4 voies dans Skylake, contre 8 auparavant, donc plus de 4 flux d'entrée, tous avec le même alignement par rapport à une page de 4k, provoqueraient des erreurs de conflit / une défaite de la prélecture. Donc, avec plus de 4 grands / les tableaux alignés sur les pages, un format entrelacé pourrait éviter ce problème.)

Pour les petits tableaux, où l'ensemble de la chose entrelacée tient dans une page 4k, oui, c'est un avantage potentiel. Sinon, c'est à peu près le même nombre total de pages touchées et d'éventuels ratés TLB, avec un 4x plus souvent au lieu de venir par groupe de 4. Cela pourrait bien être mieux pour la prélecture TLB, si elle peut faire une page à l'avance à la place d'être submergé par plusieurs ratés TLB venant en même temps.

Modification de la structure SoA:

Cela peut aider à informer le compilateur que la mémoire pointée par chacun des pointeurs ne se chevauche pas. La plupart des compilateurs C ++ (y compris les 4 compilateurs x86 majeurs, gcc / clang / MSVC / ICC) prennent en charge __restrict comme mot-clé avec la même sémantique que C99 restrict . Ou pour la portabilité, utilisez #ifdef / #define pour définir un mot-clé restrict comme vide ou __restrict ou autre, comme approprié pour le compilateur.

struct SoA_scene {
        size_t size;
        float *__restrict xs;
        float *__restrict ys;
        float *__restrict zs;
        float *__restrict ws;
};

Cela peut certainement aider avec la vectorisation automatique, sinon le compilateur ne sait pas que xs [i] = foo; code> ne change pas la valeur de ys [i + 1] pour la prochaine itération.

Si vous lisez ces variables dans des variables locales (le compilateur est donc sûr que ce pointeur les affectations ne modifient pas le pointeur lui-même dans la structure), vous pouvez les déclarer comme float * __ restrict xs = soa.xs; et ainsi de suite.

Le format entrelacé évite par nature cette possibilité d'alias.