2015-06-01 7 views
9

Quando ho avuto un processore Haswell per la prima volta ho provato a implementare FMA per determinare il set Mandelbrot. L'algoritmo principale è questo:Ottimizza per moltiplicazione rapida ma aggiunta lenta: FMA e doppio raddoppia

intn = 0; 
for(int32_t i=0; i<maxiter; i++) { 
    floatn x2 = square(x), y2 = square(y); //square(x) = x*x 
    floatn r2 = x2 + y2; 
    booln mask = r2<cut; //booln is in the float domain non integer domain 
    if(!horizontal_or(mask)) break; //_mm256_testz_pd(mask) 
    n -= mask 
    floatn t = x*y; mul2(t); //mul2(t): t*=2 
    x = x2 - y2 + cx; 
    y = t + cy; 
} 

Questo determina se n pixel sono Mandelbrot. Quindi per il doppio punto mobile viene eseguito su 4 pixel (floatn = __m256d, intn = __m256i). Ciò richiede 4 moltiplicazioni in virgola mobile SIMD e quattro aggiunte in virgola mobile SIMD.

Poi ho modificato questo per lavorare con FMA come questo

intn n = 0; 
for(int32_t i=0; i<maxiter; i++) { 
    floatn r2 = mul_add(x,x,y*y); 
    booln mask = r2<cut; 
    if(!horizontal_or(mask)) break; 
    add_mask(n,mask); 
    floatn t = x*y; 
    x = mul_sub(x,x, mul_sub(y,y,cx)); 
    y = mul_add(2.0f,t,cy); 
} 

dove mul_add chiama _mm256_fmad_pd e mul_sub chiamate _mm256_fmsub_pd. Questo metodo utilizza 4 operazioni SIMD FMA e due moltiplicazioni SIMD che sono due operazioni aritmetiche in meno quindi senza FMA. Inoltre, FMA e moltiplicazione possono utilizzare due porte e aggiungere solo una.

Per rendere i miei test meno precisi ho ingrandito una regione che è interamente nel set di Mandelbrot, quindi tutti i valori sono maxiter. In questo caso il metodo che utilizza FMA è circa il 27% più veloce. Questo è sicuramente un miglioramento, ma passare da SSE a AVX ha raddoppiato le mie prestazioni, quindi speravo di ottenere un altro fattore due con FMA.

Ma poi ho trovato this risposta per quanto riguarda FMA in cui si dice

L'aspetto importante della fusa multiply-add l'istruzione è il (quasi) infinita precisione del risultato intermedio. Questo aiuta con le prestazioni, ma non tanto perché due operazioni sono codificate in una singola istruzione - Aiuta con le prestazioni perché la precisione virtualmente infinita del risultato intermedio è talvolta importante e molto costosa da recuperare con la moltiplicazione e l'aggiunta ordinarie quando questo livello di la precisione è davvero ciò che il programmatore sta cercando.

e successivamente riportato un esempio di doppia * doppio per double-double moltiplicazione

high = a * b; /* double-precision approximation of the real product */ 
low = fma(a, b, -high); /* remainder of the real product */ 

Da questo, conclusi che ero implementando FMA non ottimale e così decisi di implementare SIMD doppia-doppia. Ho implementato il doppio doppio basato sulla carta Extended-Precision Floating-Point Numbers for GPU Computation. La carta è per il doppio float quindi l'ho modificata per il doppio doppio. Inoltre, invece di impacchettare un valore doppio doppio in un registro SIMD, impongo 4 valori double-double in un registro alto AVX e un registro basso AVX.

Per il set Mandelbrot ciò di cui ho veramente bisogno è la doppia moltiplicazione e l'aggiunta. In questo documento sono le funzioni df64_add e df64_mult. L'immagine seguente mostra l'assieme per la mia funzione df64_mult per software FMA (sinistra) e hardware FMA (destra). Questo dimostra chiaramente che l'hardware FMA è un grande miglioramento per la doppia moltiplicazione.

fma software vs hardware

Così come l'hardware FMA eseguire in doppia-doppia di Mandelbrot set calcolo? La risposta è che è solo circa il 15% più veloce rispetto al software FMA. È molto meno di quanto speravo. Il calcolo doppio-doppio di Mandelbrot richiede 4 aggiunte doppio doppio e quattro moltiplicazioni doppio doppio (x*x, y*y, x*y e 2*(x*y)). Tuttavia, il 2*(x*y) multiplication is trivial for double-double così questa moltiplicazione può essere ignorata nel costo. Pertanto, la ragione per cui ritengo che il miglioramento usando l'hardware FMA sia così piccola è che il calcolo è dominato dalla lenta aggiunta doppia doppia (vedi l'assemblaggio sotto).

Una volta la moltiplicazione era più lenta dell'aggiunta (ei programmatori utilizzavano diversi trucchi per evitare la moltiplicazione) ma con Haswell sembra che sia il contrario. Non solo a causa di FMA ma anche perché la moltiplicazione può utilizzare due porte ma aggiunta solo una.

Così le mie domande (finalmente) sono:

  1. Come si fa a ottimizzare quando aggiunta è lento rispetto alla moltiplicazione?
  2. Esiste un modo algebrico per modificare il mio algoritmo per utilizzare più moltiplicazioni e meno aggiunte? So che ci sono metodi per fare il contrario, ad es. (x+y)*(x+y) - (x*x+y*y) = 2*x*y che utilizzano altre due aggiunte per una moltiplicazione in meno.
  3. C'è un modo per semplicemente la funzione df64_add (ad esempio utilizzando FMA)?

Nel caso in cui qualcuno si chieda, il metodo del doppio doppio è circa dieci volte più lento del doppio. Non è così male, penso che se ci fosse un tipo di hardware di precisione quadrupla probabilmente sarebbe almeno due volte più lento del doppio, quindi il mio metodo software è circa cinque volte più lento di quello che mi aspetterei per l'hardware, se esistesse.

df64_add assemblaggio

vmovapd 8(%rsp), %ymm0 
movq %rdi, %rax 
vmovapd 72(%rsp), %ymm1 
vmovapd 40(%rsp), %ymm3 
vaddpd %ymm1, %ymm0, %ymm4 
vmovapd 104(%rsp), %ymm5 
vsubpd %ymm0, %ymm4, %ymm2 
vsubpd %ymm2, %ymm1, %ymm1 
vsubpd %ymm2, %ymm4, %ymm2 
vsubpd %ymm2, %ymm0, %ymm0 
vaddpd %ymm1, %ymm0, %ymm2 
vaddpd %ymm5, %ymm3, %ymm1 
vsubpd %ymm3, %ymm1, %ymm6 
vsubpd %ymm6, %ymm5, %ymm5 
vsubpd %ymm6, %ymm1, %ymm6 
vaddpd %ymm1, %ymm2, %ymm1 
vsubpd %ymm6, %ymm3, %ymm3 
vaddpd %ymm1, %ymm4, %ymm2 
vaddpd %ymm5, %ymm3, %ymm3 
vsubpd %ymm4, %ymm2, %ymm4 
vsubpd %ymm4, %ymm1, %ymm1 
vaddpd %ymm3, %ymm1, %ymm0 
vaddpd %ymm0, %ymm2, %ymm1 
vsubpd %ymm2, %ymm1, %ymm2 
vmovapd %ymm1, (%rdi) 
vsubpd %ymm2, %ymm0, %ymm0 
vmovapd %ymm0, 32(%rdi) 
vzeroupper 
ret 
+1

Come Agner Fog rileva in alcune delle sue risorse di ottimizzazione (http://www.agner.org/optimize/ ma non ricordo l'esatto file o pagina), a volte FMA rende le cose più lente.Penso di ricordare che l'esempio era 'x * x + y * y', dove, per alcuni modelli di processori Intel, la latenza, se implementata come mul-fma, rende il timing peggiore di mul-mul-add, che contiene più parallelismo. –

+1

@PascalCuoq, finora ho utilizzato [IACA] (https://stackoverflow.com/questions/26021337/what-is-iaca-and-how-do-i-use-it/26021338#26021338) e quindi guardando i tempi. Per il doppio è stato per lo più un gioco a indovinare poiché ci sono così tante permutazioni con FMA per il set Mandelbrot. Con il doppio doppio è abbastanza chiaro dove usarlo. Ma il problema principale è l'aggiunta con doppio doppio. È così lento rispetto alla moltiplicazione. Cosa si può fare? –

+2

Un'aggiunta doppia doppia programmata correttamente comprende esattamente 20 operazioni di base in virgola mobile. Potresti incontrare versioni che usano meno istruzioni, tuttavia queste perdono la precisione proprio quando ne hai bisogno e la vuoi, che è la sottrazione effettiva di operandi quasi identici in grandezza. L'accelerazione dell'aggiunta doppia doppia richiede miglioramenti HW. Le operazioni di grandezza massima e di magnitudo minima aiutano un po 'e l'aggiunta a tre ingressi con arrotondamento singolo aiuta molto. Alcune architetture di processori offrono il primo, non ne conosco nessuna che offra il secondo, – njuffa

risposta

4

di rispondere alla mia terza domanda ho trovato una soluzione più veloce per la doppia-doppia aggiunta. Ho trovato una definizione alternativa nel documento Implementation of float-float operators on graphics hardware.

Theorem 5 (Add22 theorem) Let be ah+al and bh+bl the float-float arguments of the following 
algorithm: 
Add22 (ah ,al ,bh ,bl) 
1 r = ah ⊕ bh 
2 if | ah | ≥ | bh | then 
3  s = (((ah ⊖ r) ⊕ bh) ⊕ b l) ⊕ a l 
4 e l s e 
5  s = (((bh ⊖ r) ⊕ ah) ⊕ a l) ⊕ b l 
6 (rh , r l) = add12 (r , s) 
7 return (rh , r l) 

Ecco come ho implementato questa (pseudo-codice):

static inline doubledoublen add22(doubledoublen const &a, doubledouble const &b) { 
    doublen aa,ab,ah,bh,al,bl; 
    booln mask; 
    aa = abs(a.hi);    //_mm256_and_pd 
    ab = abs(b.hi); 
    mask = aa >= ab;    //_mm256_cmple_pd 
    // z = select(cut,x,y) is a SIMD version of z = cut ? x : y; 
    ah = select(mask,a.hi,b.hi); //_mm256_blendv_pd 
    bh = select(mask,b.hi,a.hi); 
    al = select(mask,a.lo,b.lo); 
    bl = select(mask,b.lo,a.lo); 

    doublen r, s; 
    r = ah + bh; 
    s = (((ah - r) + bh) + bl) + al; 
    return two_sum(r,s); 
} 

Questa definizione di Add22 usa 11 aggiunte anziché 20 ma richiede codice aggiuntivo per determinare se |ah| >= |bh|. Here is a discussion on how to implement SIMD minmag and maxmag functions. Fortunatamente, la maggior parte del codice addizionale non utilizza la porta 1. Ora solo 12 istruzioni vanno alla porta 1, invece di 20.

Ecco una forma di analisi rendimento IACA per il nuovo Add22

Throughput Analysis Report 
-------------------------- 
Block Throughput: 12.05 Cycles  Throughput Bottleneck: Port1 

Port Binding In Cycles Per Iteration: 
--------------------------------------------------------------------------------------- 
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | 
--------------------------------------------------------------------------------------- 
| Cycles | 0.0 0.0 | 12.0 | 2.5 2.5 | 2.5 2.5 | 2.0 | 10.0 | 0.0 | 2.0 | 
--------------------------------------------------------------------------------------- 


| Num Of |     Ports pressure in cycles      | | 
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | | 
--------------------------------------------------------------------------------- 
| 1 |   |  | 0.5 0.5 | 0.5 0.5 |  |  |  |  | | vmovapd ymm3, ymmword ptr [rip] 
| 1 |   |  | 0.5 0.5 | 0.5 0.5 |  |  |  |  | | vmovapd ymm0, ymmword ptr [rdx] 
| 1 |   |  | 0.5 0.5 | 0.5 0.5 |  |  |  |  | | vmovapd ymm4, ymmword ptr [rsi] 
| 1 |   |  |   |   |  | 1.0 |  |  | | vandpd ymm2, ymm4, ymm3 
| 1 |   |  |   |   |  | 1.0 |  |  | | vandpd ymm3, ymm0, ymm3 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vcmppd ymm2, ymm3, ymm2, 0x2 
| 1 |   |  | 0.5 0.5 | 0.5 0.5 |  |  |  |  | | vmovapd ymm3, ymmword ptr [rsi+0x20] 
| 2 |   |  |   |   |  | 2.0 |  |  | | vblendvpd ymm1, ymm0, ymm4, ymm2 
| 2 |   |  |   |   |  | 2.0 |  |  | | vblendvpd ymm4, ymm4, ymm0, ymm2 
| 1 |   |  | 0.5 0.5 | 0.5 0.5 |  |  |  |  | | vmovapd ymm0, ymmword ptr [rdx+0x20] 
| 2 |   |  |   |   |  | 2.0 |  |  | | vblendvpd ymm5, ymm0, ymm3, ymm2 
| 2 |   |  |   |   |  | 2.0 |  |  | | vblendvpd ymm0, ymm3, ymm0, ymm2 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm3, ymm1, ymm4 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm2, ymm1, ymm3 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm1, ymm2, ymm4 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm1, ymm1, ymm0 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm0, ymm1, ymm5 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm2, ymm3, ymm0 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm1, ymm2, ymm3 
| 2^ |   |  |   |   | 1.0 |  |  | 1.0 | | vmovapd ymmword ptr [rdi], ymm2 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm0, ymm0, ymm1 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm1, ymm2, ymm1 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm3, ymm3, ymm1 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm0, ymm3, ymm0 
| 2^ |   |  |   |   | 1.0 |  |  | 1.0 | | vmovapd ymmword ptr [rdi+0x20], ymm0 

e qui è l'analisi di throughput dal vecchio

Throughput Analysis Report 
-------------------------- 
Block Throughput: 20.00 Cycles  Throughput Bottleneck: Port1 

Port Binding In Cycles Per Iteration: 
--------------------------------------------------------------------------------------- 
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | 
--------------------------------------------------------------------------------------- 
| Cycles | 0.0 0.0 | 20.0 | 2.0 2.0 | 2.0 2.0 | 2.0 | 0.0 | 0.0 | 2.0 | 
--------------------------------------------------------------------------------------- 

| Num Of |     Ports pressure in cycles      | | 
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | | 
--------------------------------------------------------------------------------- 
| 1 |   |  | 1.0 1.0 |   |  |  |  |  | | vmovapd ymm0, ymmword ptr [rsi] 
| 1 |   |  |   | 1.0 1.0 |  |  |  |  | | vmovapd ymm1, ymmword ptr [rdx] 
| 1 |   |  | 1.0 1.0 |   |  |  |  |  | | vmovapd ymm3, ymmword ptr [rsi+0x20] 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm4, ymm0, ymm1 
| 1 |   |  |   | 1.0 1.0 |  |  |  |  | | vmovapd ymm5, ymmword ptr [rdx+0x20] 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm2, ymm4, ymm0 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm1, ymm1, ymm2 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm2, ymm4, ymm2 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm0, ymm0, ymm2 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm2, ymm0, ymm1 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm1, ymm3, ymm5 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm6, ymm1, ymm3 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm5, ymm5, ymm6 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm6, ymm1, ymm6 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm1, ymm2, ymm1 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm3, ymm3, ymm6 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm2, ymm4, ymm1 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm3, ymm3, ymm5 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm4, ymm2, ymm4 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm1, ymm1, ymm4 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm0, ymm1, ymm3 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vaddpd ymm1, ymm2, ymm0 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm2, ymm1, ymm2 
| 2^ |   |  |   |   | 1.0 |  |  | 1.0 | | vmovapd ymmword ptr [rdi], ymm1 
| 1 |   | 1.0 |   |   |  |  |  |  | CP | vsubpd ymm0, ymm0, ymm2 
| 2^ |   |  |   |   | 1.0 |  |  | 1.0 | | vmovapd ymmword ptr [rdi+0x20], ymm0 

Una soluzione migliore sarebbe se ci fossero tre operandi istruzioni per la modalità singola di arrotondamento oltre FMA. Mi sembra che ci dovrebbero essere istruzioni in modalità rounding singolo per

a + b + c 
a * b + c //FMA - this is the only one in x86 so far 
a * b * c 
+1

Cosa ti dà 'a * b * c'? La sum-of-three rende banale la doppia addizione doppia e FMA rende banale la doppia moltiplicazione doppia, ma 'a * b * c' non sembra avere uno scopo ovvio come quello. – tmyklebu

+0

@tmyklebu, non lo so. Forse hai ragione. Per il doppio doppio sono d'accordo con te sul fatto che 'a * b * c' non sembra avere un punto. Mi chiedo quali altre applicazioni siano utili per queste tre istruzioni per la modalità rounding dell'operando? L'ho usato solo per il doppio doppio. –

+0

@tmyklebu, in realtà vieni a pensarci solo quelli che sono necessari sono 'a + b - c' e' a * b - c'. –

Problemi correlati