Cette optimisation peut donner un gain de vitesse par rapport à ce que l'on peut faire de mieux en réévaluant le polynôme séparément pour chaque i . Mais seulement si vous le généralisez à un pas plus grand, pour avoir encore assez de parallélisme dans la boucle. Ma version peut stocker 1 vecteur (4 doubles) par cycle d'horloge sur mon Skylake. . Ou en théorie 1 vecteur pour 2 horloges sur les Intel antérieurs, y compris Sandybridge où un stockage de 256 bits prend 2 cycles d'horloge dans le port de stockage.
Il est évident qu'il doit tenir dans le cache pour fonctionner aussi rapidement, donc répétez un petit LEN plusieurs fois.
Ce site optimisation de la réduction de la force (il suffit d'ajouter au lieu de commencer avec un nouveau i et en multipliant) introduit une dépendance sérielle entre les itérations de la boucle impliquant des calculs de type FP plutôt que des incréments entiers.
L'original est a parallélisme des données sur chaque élément de sortie chacun ne dépend que des constantes et de sa propre i valeur. Les compilateurs peuvent auto-vectoriser avec SIMD (SSE2, ou AVX si vous utilisez -O3 -march=native ), et les CPU peuvent chevaucher le travail entre les itérations de la boucle avec une exécution hors ordre. Malgré la quantité de travail supplémentaire, le CPU est capable d'appliquer une force brute suffisante, avec l'aide du compilateur.
Mais la version qui calcule poly(i+1) en termes de poly(i) a un parallélisme très limité ; pas de vectorisation SIMD, et votre CPU ne peut exécuter que deux additions scalaires par 4 cycles, par exemple, où 4 cycles est la latence de l'addition FP sur Intel Skylake à Tiger Lake. ( https://uops.info/ ).
La réponse de @huseyin tugrul buyukisik montre comment on peut se rapprocher du débit maximal de la version originale sur un CPU plus moderne, avec deux opérations FMA pour évaluer le polynôme (schéma de Horner), plus une conversion int->FP ou un incrément FP. (Ce dernier crée une chaîne FP dep qu'il faut dérouler pour la cacher).
Dans le meilleur des cas, vous avez donc 3 opérations mathématiques FP par vecteur SIMD de sortie. (Plus un stockage). Les CPU Intel actuels n'ont que deux unités d'exécution FP qui peuvent exécuter des opérations mathématiques FP. (Avec des vecteurs de 512 bits avec AVX-512, les CPU actuels ferment l'ALU vectorielle sur le port 1, il n'y a donc que 2 ports SIMD ALU, donc les opérations mathématiques non FP comme l'incrémentation d'un entier seront également en concurrence pour le débit).
Depuis Zen2, AMD a deux unités FMA/mul sur deux ports, et deux unités FP add/sub sur deux ports différents, donc si vous utilisez FMA pour faire l'addition, vous avez un maximum théorique de quatre additions SIMD par cycle d'horloge.
Haswell/Broadwell ont un FMA à 2/clock, mais seulement un FP add/sub à 1/clock (avec une latence plus faible). C'est une bonne chose pour le code naïf, pas génial pour le code qui a été optimisé pour avoir beaucoup de parallélisme. C'est probablement pour cela qu'Intel l'a modifié dans Skylake.
Vos processeurs Sandybridge (E5-1620) et Nehalem (W5580) ont un FP add/sub 1/clock, un FP mul 1/clock, sur des ports séparés. C'est ce sur quoi Haswell s'est basé. Et pourquoi l'ajout de multiplications supplémentaires n'est pas un gros problème : elles peuvent fonctionner en parallèle avec les additions existantes.
Trouver le parallélisme : réduction de la force avec une foulée arbitraire
Votre compute2 calcule le prochain Y et le prochain Z en fonction de la valeur immédiatement précédente, c'est-à-dire qu'avec un stride de 1, les valeurs dont vous avez besoin pour data[i+1] . Ainsi, chaque itération est dépendante de celle qui la précède immédiatement.
Si vous généralisez cela à d'autres foulées, vous pouvez avancer 4, 8, 16 ou plus de valeurs Y et Z distinctes, de sorte qu'elles sautent toutes au même rythme, indépendamment les unes des autres. Cela permet de retrouver un certain parallélisme dont le compilateur et/ou le CPU peuvent tirer parti.
poly(i) = A i^2 + B i + C
poly(i+4) = A (i+4)^2 + B (i+4) + C
= A*i^2 + A*8i + A*16 + B i + 4B + C
= poly(i) + 4*(2A)*i + 16*A + 4B
(J'ai choisi un 4 en partie parce que l'écriture des exposants est désordonnée dans un texte ASCII. Sur le papier, j'aurais pu simplement faire (i+s)^2 = i^2 + 2is + s^2 et ainsi de suite. Une autre approche consisterait à prendre les différences entre les points successifs, puis les différences de ces différences. Du moins, c'est là que nous voulons aboutir, avec des valeurs que nous pouvons ajouter pour générer la séquence (de valeurs FP). En gros, cela revient à prendre la 1ère et la 2ème dérivée, et ensuite la boucle optimisée est effectivement en train d'intégrer pour récupérer la fonction originale).
- La partie de la 2ème dérivée (anciennement
Z += A*2 )
devient Z[j] += A * (stride * 2); de la A * 8i terme.
- La croissance linéaire (autrefois le
Z = A+B; )
devient (A*stride + B) * stride de la 16*A + 4*B sur l'avancement par 4.
Comme vous pouvez le constater, pour stride=1 (pas de déroulage), ceux-ci se simplifient aux valeurs originales. Mais avec des stride mais ils laissent quand même cette méthode fonctionner, sans multiplications supplémentaires ou autre chose dans la boucle.
// UNTESTED, but compiles to asm that looks right.
// Hopefully I got the math right.
#include <stdalign.h>
#include <stdlib.h>
#define LEN 1024
alignas(64) double data[LEN];
void compute2_unrolled(void)
{
const double A = 1.1, B = 2.2, C = 3.3;
const int stride = 16; // unroll factor. 1 reduces to the original
const double deriv2 = A * (stride * 2);
double Z[stride], Y[stride];
for (int j = 0 ; j<stride ; j++){ // this loop will fully unroll
Y[j] = j*j*A + j*B + C; // poly(j) starting values to increment
Z[j] = (A*stride + B) * stride; // 1st derivative = same initial Z for all elements
}
for(size_t i=0; i < LEN - (stride-1); i+=stride) {
#if 0
for (int j = 0 ; j<stride ; j++){
data[i+j] = Y[j];
Y[j] += Z[j];
Z[j] += deriv2;
}
#else
// loops that are easy(?) for a compiler to roll up into some SIMD vectors
for (int j=0 ; j<stride ; j++) data[i+j] = Y[j]; // store
for (int j=0 ; j<stride ; j++) Y[j] += Z[j]; // add
for (int j=0 ; j<stride ; j++) Z[j] += deriv2; // add
#endif
}
// cleanup for the last few i values
for (int j = 0 ; j < LEN % stride ; j++) {
// letting the compiler see the %stride can help it decide *not* to auto-vectorize this part for non-huge strides.
size_t i = LEN - (stride-1) + j;
//data[i] = poly(i);
}
}
( Explorateur de compilateur Godbolt ) En fait, cela permet d'auto-vectoriser joliment avec stride=16 (4x YMM vecteurs de 4 double s chacun) avec clang14 -O3 -march=skylake -ffast-math .
Il semble que clang ait encore déroulé par 2 (mais sans inventer de nouveaux accumulateurs, bien sûr). Donc chaque itération de la boucle asm 2x 8 vaddpd instructions et 2x 4 magasins. Et un add/jne macro-fusion de l'overhead de la boucle, puisque clang compte vers zéro, en indexant depuis la fin du tableau. Donc, avec les hits de cache L1d pour les magasins, cela devrait fonctionner (sur Skylake) à 1 magasin par cycle d'horloge, faisant également 2x vaddpd chaque cycle. C'est le débit maximum pour ces deux choses. Le frontal n'a besoin de suivre qu'un peu plus de 3 uops / cycle d'horloge, mais il est large de 4 depuis Core2, et le cache uop de la famille Sandybridge fait que cela ne pose aucun problème. (Sauf si vous vous heurtez à l'erratum JCC sur Skylake, alors j'ai utilisé -mbranches-within-32B-boundaries d'avoir des instructions clang pad pour éviter cela .)
Je ne l'ai pas vraiment testé sur mon bureau. vaddpd La latence est de 4 cycles, donc 4 chaînes de dépannage sont tout juste suffisantes pour maintenir en vol 4 opérations indépendantes. Et avec deux chaînes dep imbriquées, il est probablement facile pour un ordonnancement imparfait d'entraîner des cycles perdus. Donc un peu plus de déroulage serait approprié si nous pouvions obtenir du compilateur qu'il fasse encore un asm sympa pour stride=32 .
(Je ne suis pas sûr de savoir pourquoi -ffast-math est nécessaire pour que clang puisse vectoriser décemment avec cette version ; je ne pense pas que clang ré-associe quoi que ce soit, il fait juste toutes les opérations mathématiques dans l'ordre de la source, à moins qu'il y ait quelque chose qui m'échappe ; je n'ai pas regardé de près les numéros de registre. Le source a été soigneusement écrit d'une manière SIMD-friendly pour rendre possible la compilation vers un asm sympa. Mise à jour : une autre version qui prend un pointeur comme argument permet de vectoriser avec ou sans -ffast-math .)
L'autre façon d'écrire les boucles est avec une boucle interne au lieu de toutes les opérations Y, puis toutes les opérations Z. clang n'auto-vectorise pas ceci, cependant, même avec -ffast-math seulement le style SIMD ci-dessus. Je l'ai inclus dans un #if 0 / #else / #endif bloc sur Godbolt. Mise à jour : dans la version prenant un pointeur arg, ceci vectorise mieux dans certains cas.
for(int i=0; i < LEN - (stride-1); i+=stride) {
for (int j = 0 ; j<stride ; j++){ // doesn't auto-vectorize with clang or GCC
data[i+j] = Y[j];
Y[j] += Z[j];
Z[j] += deriv2;
}
}
Nous devons choisir manuellement un montant de déroulement approprié. . Un facteur d'unroll trop important peut même empêcher le compilateur de voir ce qui se passe et de garder les tableaux temp dans les registres. par exemple 32 est un problème pour clang, mais pas 24 . Il peut y avoir des options de réglage pour forcer le compilateur à dérouler les boucles jusqu'à un certain nombre ; il en existe pour GCC qui peuvent parfois être utilisées pour lui permettre de voir quelque chose au moment de la compilation.
(Une autre approche serait la vectorisation manuelle avec #include <immintrin.h> y __m256d Z[4] au lieu de double Z[16] . Mais cette version peut vectoriser pour d'autres ISA comme AArch64).
D'autres inconvénients d'un facteur de déroulement important sont de laisser plus de travail de nettoyage lorsque la taille du problème n'est pas un multiple du déroulement. (Vous pouvez utiliser le compute1 pour le nettoyage, laissant le compilateur vectoriser cela pendant une itération ou deux avant de faire du scalaire).
En théorie, un compilateur devrait être autorisé de le faire pour vous avec -ffast-math soit de compute1 en effectuant la réduction de force sur le polynôme d'origine, ou de compute2 voir comment la foulée s'accumule. Mais en pratique, c'est très compliqué et les humains doivent le faire eux-mêmes. A moins que/jusqu'à ce que quelqu'un apprenne aux compilateurs comment rechercher des modèles comme celui-ci et appliquer le calcul différentiel (ou les différences successives avec un choix de stride).
Résultats expérimentaux :
Sur mon ordinateur de bureau (i7-6700k) avec clang13.0.0, j'ai testé et j'ai été agréablement surpris de constater qu'il fonctionne en fait à un stockage SIMD par cycle d'horloge avec la combinaison d'options de compilation (fast-math ou non) et de #if 0 vs. #if 1 sur la stratégie de la boucle intérieure, avec une version du cadre de référence de @huseyin tugrul buyukisik améliorée pour répéter une quantité plus mesurable entre rdtsc des instructions. Et pour compenser la différence entre la fréquence de l'horloge centrale et celle de l'horloge de l'ordinateur. fréquence "de référence" de l'horloge rdtsc lit Dans mon cas, 3.9GHz contre 4008 MHz. (Le turbo max nominal est de 4.2GHz, mais avec EPP = 1,5GHz). balance_performance sous Linux, il ne veut tourner que jusqu'à 3,9 GHz).
Code source que j'ai testé sur Godbolt : utiliser une seule boucle interne, plutôt que 3 boucles séparées. j<16 les boucles, et pas en utilisant -ffast-math . Et en utilisant __attribute__((noinline)) pour éviter qu'il ne s'insère dans la boucle de répétition. D'autres variations des options et de la source ont conduit à quelques vpermpd se déplace à l'intérieur de la boucle.
$ clang++ -std=gnu++17 -O3 -march=native -mbranches-within-32B-boundaries poly-eval.cpp -Wall
# warning about noipa, only GCC knows that attribute
$ perf stat --all-user -etask-clock,context-switches,cpu-migrations,page-faults,cycles,instructions,uops_issued.any,uops_executed.thread,fp_arith_inst_retired.256b_packed_double -r10 ./a.out
... (10 runs of the whole program, ending with)
0.252072 cycles per data element (corrected from ref cycles to core clocks)
0.252236 cycles per data element (corrected from ref cycles to core clocks)
0.252408 cycles per data element (corrected from ref cycles to core clocks)
0.251945 cycles per data element (corrected from ref cycles to core clocks)
0.252442 cycles per data element (corrected from ref cycles to core clocks)
0.252295 cycles per data element (corrected from ref cycles to core clocks)
0.252109 cycles per data element (corrected from ref cycles to core clocks)
xor=4303
min cycles per data = 0.251868
Performance counter stats for './a.out' (10 runs):
298.92 msec task-clock # 0.989 CPUs utilized ( +- 0.49% )
0 context-switches # 0.000 /sec
0 cpu-migrations # 0.000 /sec
129 page-faults # 427.583 /sec ( +- 0.56% )
1,162,430,637 cycles # 3.853 GHz ( +- 0.49% ) # time spent in the kernel for system calls and interrupts isn't counted, that's why it's not 3.90 GHz
3,772,516,605 instructions # 3.22 insn per cycle ( +- 0.00% )
3,683,072,459 uops_issued.any # 12.208 G/sec ( +- 0.00% )
4,824,064,881 uops_executed.thread # 15.990 G/sec ( +- 0.00% )
2,304,000,000 fp_arith_inst_retired.256b_packed_double # 7.637 G/sec
0.30210 +- 0.00152 seconds time elapsed ( +- 0.50% )
fp_arith_inst_retired.256b_packed_double compte 1 pour chaque instruction FP add ou mul (2 pour FMA), donc nous obtenons 1,98 vaddpd instructions par cycle d'horloge pour l'ensemble du programme, y compris l'impression et ainsi de suite. C'est très proche du maximum théorique de 2/clock, ne souffrant apparemment pas d'un ordonnancement sous-optimal des instructions d'addition. (J'ai augmenté la boucle de répétition pour qu'elle y passe plus de temps, il est donc plus utile d'utiliser perf stat sur l'ensemble du programme).
Le but de cette optimisation était d'effectuer le même travail avec moins de FLOPS, mais cela signifie également que nous maximisons essentiellement la limite de 8 FLOPS/clock pour Skylake sans utiliser la FMA. (Et nous obtenons 30,58 GFLOP/s à 3,9 GHz sur un seul cœur).
Asm de la fonction non-inline ( objdump -drwC -Mintel ) ; clang a utilisé 4 ensembles de 2 paires de vecteurs YMM, et a déroulé la boucle 3 fois de plus. Notez la add rax,0x30 faire 3 * 0x10 doubles par itération. Cela donne donc un multiple exact de 24KiB sans nettoyage.
0000000000001440 <void compute2<3072>(double*)>:
# just loading constants; the setup loop did flatten out and optimize into vector constants
1440: c5 fd 28 0d 18 0c 00 00 vmovapd ymm1,YMMWORD PTR [rip+0xc18] # 2060 <_IO_stdin_used+0x60>
1448: c5 fd 28 15 30 0c 00 00 vmovapd ymm2,YMMWORD PTR [rip+0xc30] # 2080 <_IO_stdin_used+0x80>
1450: c5 fd 28 1d 48 0c 00 00 vmovapd ymm3,YMMWORD PTR [rip+0xc48] # 20a0 <_IO_stdin_used+0xa0>
1458: c4 e2 7d 19 25 bf 0b 00 00 vbroadcastsd ymm4,QWORD PTR [rip+0xbbf] # 2020 <_IO_stdin_used+0x20>
1461: c5 fd 28 2d 57 0c 00 00 vmovapd ymm5,YMMWORD PTR [rip+0xc57] # 20c0 <_IO_stdin_used+0xc0>
1469: 48 c7 c0 d0 ff ff ff mov rax,0xffffffffffffffd0
1470: c4 e2 7d 19 05 af 0b 00 00 vbroadcastsd ymm0,QWORD PTR [rip+0xbaf] # 2028 <_IO_stdin_used+0x28>
1479: c5 fd 28 f4 vmovapd ymm6,ymm4
147d: c5 fd 28 fc vmovapd ymm7,ymm4
1481: c5 7d 28 c4 vmovapd ymm8,ymm4
1485: 66 66 2e 0f 1f 84 00 00 00 00 00 data16 cs nop WORD PTR [rax+rax*1+0x0]
# top of outer loop. The NOP before this is to align it.
1490: c5 fd 11 ac c7 80 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x180],ymm5
1499: c5 d5 58 ec vaddpd ymm5,ymm5,ymm4
149d: c5 dd 58 e0 vaddpd ymm4,ymm4,ymm0
14a1: c5 fd 11 9c c7 a0 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x1a0],ymm3
14aa: c5 e5 58 de vaddpd ymm3,ymm3,ymm6
14ae: c5 cd 58 f0 vaddpd ymm6,ymm6,ymm0
14b2: c5 fd 11 94 c7 c0 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x1c0],ymm2
14bb: c5 ed 58 d7 vaddpd ymm2,ymm2,ymm7
14bf: c5 c5 58 f8 vaddpd ymm7,ymm7,ymm0
14c3: c5 fd 11 8c c7 e0 01 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x1e0],ymm1
14cc: c5 bd 58 c9 vaddpd ymm1,ymm8,ymm1
14d0: c5 3d 58 c0 vaddpd ymm8,ymm8,ymm0
14d4: c5 fd 11 ac c7 00 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x200],ymm5
14dd: c5 d5 58 ec vaddpd ymm5,ymm5,ymm4
14e1: c5 dd 58 e0 vaddpd ymm4,ymm4,ymm0
14e5: c5 fd 11 9c c7 20 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x220],ymm3
14ee: c5 e5 58 de vaddpd ymm3,ymm3,ymm6
14f2: c5 cd 58 f0 vaddpd ymm6,ymm6,ymm0
14f6: c5 fd 11 94 c7 40 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x240],ymm2
14ff: c5 ed 58 d7 vaddpd ymm2,ymm2,ymm7
1503: c5 c5 58 f8 vaddpd ymm7,ymm7,ymm0
1507: c5 fd 11 8c c7 60 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x260],ymm1
1510: c5 bd 58 c9 vaddpd ymm1,ymm8,ymm1
1514: c5 3d 58 c0 vaddpd ymm8,ymm8,ymm0
1518: c5 fd 11 ac c7 80 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x280],ymm5
1521: c5 d5 58 ec vaddpd ymm5,ymm5,ymm4
1525: c5 dd 58 e0 vaddpd ymm4,ymm4,ymm0
1529: c5 fd 11 9c c7 a0 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x2a0],ymm3
1532: c5 e5 58 de vaddpd ymm3,ymm3,ymm6
1536: c5 cd 58 f0 vaddpd ymm6,ymm6,ymm0
153a: c5 fd 11 94 c7 c0 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x2c0],ymm2
1543: c5 ed 58 d7 vaddpd ymm2,ymm2,ymm7
1547: c5 c5 58 f8 vaddpd ymm7,ymm7,ymm0
154b: c5 fd 11 8c c7 e0 02 00 00 vmovupd YMMWORD PTR [rdi+rax*8+0x2e0],ymm1
1554: c5 bd 58 c9 vaddpd ymm1,ymm8,ymm1
1558: c5 3d 58 c0 vaddpd ymm8,ymm8,ymm0
155c: 48 83 c0 30 add rax,0x30
1560: 48 3d c1 0b 00 00 cmp rax,0xbc1
1566: 0f 82 24 ff ff ff jb 1490 <void compute2<3072>(double*)+0x50>
156c: c5 f8 77 vzeroupper
156f: c3 ret
En rapport :
-
Limites de latence et limites de débit des processeurs pour les opérations qui doivent se produire en séquence. - analyse du code avec deux chaînes de dépouillement, l'une lisant l'autre et antérieure à elle-même. Même modèle de dépendance que la boucle réduite en force, sauf que l'une de ses chaînes est un multiplicateur FP. (C'est aussi un schéma d'évaluation polynomiale, mais pour un grand polynôme).
-
Optimisation SIMD d'une courbe calculée à partir de la dérivée seconde Un autre cas où l'on peut avancer à grands pas dans la dépendance sérielle.
-
Est-il possible d'utiliser SIMD sur une dépendance série dans un calcul, comme un filtre de moyenne mobile exponentielle ? - S'il existe une formule fermée pour les n étapes à venir, vous pouvez l'utiliser pour contourner les dépendances en série.
-
Exécution en désordre, comment résoudre la dépendance réelle ? - Les processeurs doivent attendre lorsqu'une instruction dépend d'une autre qui n'a pas encore été exécutée.
-
Analyse de la chaîne de dépendance une analyse de la chaîne de dépendance sans boucle, à partir d'un des exemples d'Agner Fog.
-
Microprocesseurs modernes Un guide de 90 minutes ! - contexte général sur l'exécution hors ordre et les pipelines. La technologie moderne SIMD à vecteur court de type CPU existe sous cette forme pour faire passer plus de travail dans le pipeline d'un seul CPU sans élargir le pipeline. En revanche, les GPU disposent de nombreux pipelines simples.
-
Pourquoi mulss ne prend que 3 cycles sur Haswell, différent des tables d'instructions d'Agner (dérouler des boucles FP avec accumulateurs multiples) ? - Quelques chiffres expérimentaux avec le déroulage pour masquer la latence des chaînes de dépendance de la FP, et quelques informations sur l'architecture du CPU concernant le renommage des registres.
