Pourquoi std::fill(0) est-il plus lent que std::fill(1) ?

Question

J'ai observé sur un système que std::fill sur un grand std::vector était significativement et constamment plus lent lors de la définition d'une valeur constante 0 par rapport à une valeur constante 1 ou une valeur dynamique :

5.8 Gio/s vs 7.5 Gio/s

Cependant, les résultats sont différents pour des tailles de données plus petites, où fill(0) est plus rapide :

Avec plus d'un thread, à une taille de données de 4 Gio, fill(1) montre une pente plus élevée, mais atteint un pic beaucoup plus bas que fill(0) (51 Gio/s vs 90 Gio/s) :

Cela soulève la question secondaire, pourquoi la bande passante maximale de fill(1) est beaucoup plus basse.

Le système de test pour cela était un CPU Intel Xeon E5-2680 v3 double socket réglé à 2,5 GHz (via /sys/cpufreq) avec 8x16 Gio DDR4-2133. J'ai testé avec GCC 6.1.0 (-O3) et le compilateur Intel 17.0.1 (-fast), les deux donnent des résultats identiques. GOMP_CPU_AFFINITY=0,12,1,13,2,14,3,15,4,16,5,17,6,18,7,19,8,20,9,21,10,22,11,23 était défini. Strem/add/24 threads obtient 85 Gio/s sur le système.

J'ai pu reproduire cet effet sur un système de serveur double socket Haswell différent, mais pas sur une autre architecture. Par exemple, sur Sandy Bridge EP, les performances de la mémoire sont identiques, tandis que dans le cache fill(0) est bien plus rapide.

Voici le code pour reproduire :

#include 
#include 
#include 
#include 
#include 

using value = int;
using vector = std::vector;

constexpr size_t write_size = 8ll * 1024 * 1024 * 1024;
constexpr size_t max_data_size = 4ll * 1024 * 1024 * 1024;

void __attribute__((noinline)) fill0(vector& v) {
    std::fill(v.begin(), v.end(), 0);
}

void __attribute__((noinline)) fill1(vector& v) {
    std::fill(v.begin(), v.end(), 1);
}

void bench(size_t data_size, int nthreads) {
#pragma omp parallel num_threads(nthreads)
    {
        vector v(data_size / (sizeof(value) * nthreads));
        auto repeat = write_size / data_size;
#pragma omp barrier
        auto t0 = omp_get_wtime();
        for (auto r = 0; r < repeat; r++)
            fill0(v);
#pragma omp barrier
        auto t1 = omp_get_wtime();
        for (auto r = 0; r < repeat; r++)
            fill1(v);
#pragma omp barrier
        auto t2 = omp_get_wtime();
#pragma omp master
        std::cout << data_size << ", " << nthreads << ", " << write_size / (t1 - t0) << ", "
                  << write_size / (t2 - t1) << "\n";
    }
}

int main(int argc, const char* argv[]) {
    std::cout << "size,nthreads,fill0,fill1\n";
    for (size_t bytes = 1024; bytes <= max_data_size; bytes *= 2) {
        bench(bytes, 1);
    }
    for (size_t bytes = 1024; bytes <= max_data_size; bytes *= 2) {
        bench(bytes, omp_get_max_threads());
    }
    for (int nthreads = 1; nthreads <= omp_get_max_threads(); nthreads++) {
        bench(max_data_size, nthreads);
    }
}

Résultats présentés compilés avec g++ fillbench.cpp -O3 -o fillbench_gcc -fopenmp.

Answer 1

De votre question + l'asm généré par le compilateur de votre réponse :

• fill(0) est un ERMSB rep stosb qui utilisera des stores de 256b dans une boucle microcodée optimisée. (Fonctionne mieux si le tampon est aligné, probablement au moins à 32B ou peut-être à 64B).

• fill(1) est une simple boucle de stockage de vecteur movaps de 128 bits. Seul un magasin peut s'exécuter par cycle d'horloge du cœur quelle que soit sa largeur, jusqu'à 256b AVX. Ainsi, les magasins de 128b ne peuvent remplir que la moitié de la largeur de bande d'écriture de cache L1D de Haswell. C'est pourquoi fill(0) est environ 2x plus rapide pour des tampons allant jusqu'à ~32kiB. Compilez avec -march=haswell ou -march=native pour corriger cela.

  Haswell peut à peine suivre les frais généraux de la boucle, mais il peut toujours exécuter 1 magasin par horloge même s'il n'est pas du tout déroulé. Mais avec 4 uops de domaine fusionné par horloge, c'est beaucoup de remplissage prenant de la place dans la fenêtre hors service. Un déroulement pourrait peut-être permettre aux fautes TLB de commencer à se résoudre plus loin de l'endroit où les magasins se produisent, puisqu'il y a plus de débit pour les uops d'adresse de magasin que pour les données de magasin. Le déroulement pourrait aider à combler le reste de la différence entre ERMSB et cette boucle de vecteur pour les tampons qui rentrent dans L1D. (Un commentaire sur la question dit que -march=native n'a aidé que fill(1) pour L1.)

Notez que rep movsd (qui pourrait être utilisé pour implémenter fill(1) pour des éléments int) performera probablement de la même manière que rep stosb sur Haswell. Bien que seule la documentation officielle garantisse que ERMSB offre un rep stosb rapide (mais pas rep stosd), les CPU réels prenant en charge ERMSB utilisent probablement un microcode tout aussi efficace pour rep stosd. Il y a un doute concernant IvyBridge, où peut-être seul b est rapide. Voir la réponse ERMSB de @BeeOnRope pour des mises à jour sur ce sujet.

gcc dispose de certaines options d'optimisation x86 pour les opérations de chaîne (comme -mstringop-strategy=alg et -mmemset-strategy=strategy), mais je ne sais pas s'il y en a qui permettent d'émettre réellement rep movsd pour fill(1). Probablement pas, car je suppose que le code commence initialement par une boucle, plutôt qu'un memset.

---

Avec plus d'un thread, à une taille de données de 4 GiB, fill(1) montre une pente plus élevée, mais atteint un pic beaucoup plus bas que fill(0) (51 Go/s contre 90 Go/s) :

Un stockage movaps normal sur une ligne de cache froide déclenche un Read For Ownership (RFO). Beaucoup de bande passante réelle de DRAM est dépensée pour lire des lignes de cache depuis la mémoire lorsque movaps écrit les premiers 16 octets. Les stores ERMSB utilisent un protocole sans RFO pour leurs stores, donc les contrôleurs mémoire n'écrivent que. (Sauf pour les lectures diverses, comme les tables de pages si des miss de page-walks se produisent même dans le cache L3, et peut-être quelques miss de chargement dans les gestionnaires d'interruption ou autre).

@BeeOnRope explique dans les commentaires que la différence entre les stores RFO réguliers et le protocole évitant les RFO utilisé par ERMSB a des inconvénients pour certaines plages de tailles de tampon sur les CPU serveur où il y a une latence élevée dans l'uncore/cache L3. Voir également la réponse ERMSB liée pour plus d'informations sur RFO vs non-RFO, et la haute latence de l'uncore (L3/mémoire) dans les CPU Intel à plusieurs cœurs posant problème pour la bande passante d'un seul cœur.

---

Les stores movntps (_mm_stream_ps()) sont faiblement ordonnés, donc ils peuvent contourner le cache et aller directement en mémoire une ligne de cache entière à la fois sans jamais lire la ligne de cache dans L1D. movntps évite les RFO, comme le fait rep stos. (rep stos stores peuvent se réorganiser les uns avec les autres, mais pas en dehors des limites de l'instruction.)

Vos résultats avec movntps dans votre réponse mise à jour sont surprenants.
Pour un seul thread avec de grands tampons, vos résultats sont movnt >> RFO régulier > ERMSB. C'est donc vraiment étrange que les deux méthodes sans RFO soient aux extrémités des simples anciens magasins, et que ERMSB soit si éloigné de l'optimal. Je n'ai actuellement pas d'explication pour cela. (modifications bienvenues avec une explication + des preuves convaincantes).

Comme nous le pensions, movnt permet à plusieurs threads d'atteindre une bande passante d'écriture agrégée élevée, comme ERMSB. movnt va toujours directement dans les buffers de remplissage de ligne puis en mémoire, il est donc beaucoup plus lent pour les tailles de tampon qui rentrent dans le cache. Un vecteur de 128b par horloge est suffisant pour saturer facilement la bande passante sans RFO d'un seul cœur vers DRAM. Probablement que vmovntps ymm (256b) est seulement un avantage mesurable par rapport à vmovntps xmm (128b) lorsqu'il stocke les résultats d'un calcul vectoriel AVX 256b lié au processeur (c'est-à-dire uniquement quand il évite la tâche de déballage en 128b).

La bande passante de movnti est faible parce que le stockage par tranches de 4 octets est un goulot d'étranglement sur 1 uop de stockage par horloge ajoutant des données aux buffers de remplissage de ligne, et non sur l'envoi de ces buffers remplis de lignes vers DRAM (jusqu'à ce que vous ayez suffisamment de threads pour saturer la bande passante mémoire).

---

@osgx a posté quelques liens intéressants dans les commentaires:

• Guide d'optimisation asm, tables d'instructions et guide microarch d'Agner Fog : http://agner.org/optimize/

• Guide d'optimisation Intel : http://www.intel.com/content/dam/www/public/us/en/documents/manuals/64-ia-32-architectures-optimization-manual.pdf.

• Snooping NUMA : http://frankdenneman.nl/2016/07/11/numa-deep-dive-part-3-cache-coherency/

• https://software.intel.com/en-us/articles/intelr-memory-latency-checker

• Protocole de cohérence de cache et performances mémoire de l'architecture Intel Haswell-EP

Voir également d'autres choses dans la balise x86 tag wiki.

Answer 2

Je vais partager mes conclusions préliminaires, dans l'espoir de encourager des réponses plus détaillées. J'ai juste pensé que ce serait trop comme partie de la question elle-même.

Le compilateur optimise fill(0) vers un memset interne. Il ne peut pas faire la même chose pour fill(1), car memset ne fonctionne que sur des octets.

Plus précisément, à la fois __memset_avx2 et __intel_avx_rep_memset de glibc sont implémentés avec une seule instruction chaude :

rep    stos %al,%es:(%rdi)

Alors que la boucle manuelle se compile en une instruction réelle de 128 bits :

add    $0x1,%rax                                                                                                       
add    $0x10,%rdx                                                                                                      
movaps %xmm0,-0x10(%rdx)                                                                                               
cmp    %rax,%r8                                                                                                        
ja     400f41

Intéressant de noter qu'il y a une optimisation de modèle/en-tête pour implémenter std::fill via memset pour les types d'octets, mais dans ce cas, c'est une optimisation du compilateur pour transformer la boucle réelle.

Étrangement, pour un std::vector, gcc commence également à optimiser fill(1). Le compilateur Intel ne le fait pas, malgré la spécification du modèle memset.

Étant donné que cela se produit uniquement lorsque le code fonctionne réellement en mémoire plutôt qu'en cache, il semble que l'architecture Haswell-EP échoue à consolider efficacement les écritures de bytes.

Je serais reconnaissance pour toute autre perspective sur le problème et les détails de la micro-architecture associée. Il n'est pas clair pour moi pourquoi cela se comporte si différemment pour quatre threads ou plus et pourquoi memset est tellement plus rapide en cache.

Mise à jour :

Voici un résultat en comparaison avec

• fill(1) qui utilise -march=native (avx2 vmovdq %ymm0) - il fonctionne mieux en L1, mais de manière similaire à la version movaps %xmm0 pour d'autres niveaux de mémoire.
• Variantes de 32, 128 et 256 bits pour les envois non-temporels. Ils ont des performances constantes avec la même performance quelle que soit la taille des données. Tous surpassent les autres variantes en mémoire, en particulier pour de petits nombres de threads. Les 128 bits et les 256 bits sont exactement similaires, pour un faible nombre de threads le 32 bits se comporte nettement moins bien.

Pour <= 6 threads, vmovnt a un avantage de 2x sur rep stos lorsqu'il opère en mémoire.

Bande passante mono-thread :

Bande passante agrégée en mémoire :

Voici le code utilisé pour les tests supplémentaires avec leurs boucles chaudes respectives :

void __attribute__ ((noinline)) fill1(vector& v) {
    std::fill(v.begin(), v.end(), 1);
}
add    $0x1,%rax
  vmovdq %ymm0,(%rdx)
  add    $0x20,%rdx
  cmp    %rdi,%rax
jb     e0

void __attribute__ ((noinline)) fill1_nt_si32(vector& v) {
    for (auto& elem : v) {
       _mm_stream_si32(&elem, 1);
    }
}
movnti %ecx,(%rax)
  add    $0x4,%rax
  cmp    %rdx,%rax
jne    18

void __attribute__ ((noinline)) fill1_nt_si128(vector& v) {
    assert((long)v.data() % 32 == 0); // alignement
    const __m128i buf = _mm_set1_epi32(1);
    size_t i;
    int* data;
    int* end4 = &v[v.size() - (v.size() % 4)];
    int* end = &v[v.size()];
    for (data = v.data(); data < end4; data += 4) {
        _mm_stream_si128((__m128i*)data, buf);
    }
    for (; data < end; data++) {
        *data = 1;
    }
}
vmovnt %xmm0,(%rdx)
  add    $0x10,%rdx
  cmp    %rcx,%rdx
jb     40

void __attribute__ ((noinline)) fill1_nt_si256(vector& v) {
    assert((long)v.data() % 32 == 0); // alignement
    const __m256i buf = _mm256_set1_epi32(1);
    size_t i;
    int* data;
    int* end8 = &v[v.size() - (v.size() % 8)];
    int* end = &v[v.size()];
    for (data = v.data(); data < end8; data += 8) {
        _mm256_stream_si256((__m256i*)data, buf);
    }
    for (; data < end; data++) {
        *data = 1;
    }
}
vmovnt %ymm0,(%rdx)
  add    $0x20,%rdx
  cmp    %rcx,%rdx
jb     40

Remarque : J'ai dû effectuer un calcul de pointeur manuel pour rendre les boucles aussi compactes. Sinon, il ferait de l'indexation vectorielle dans la boucle, probablement en raison de l'optimisation confondante de l'intrinsèque.