ottimizzata 2x2 moltiplicazione matriciale: assemblaggio lenta rispetto SIMD veloce

Question

Problemi

sto studiando algoritmi matrice moltiplicazione alte prestazioni come OpenBLAS o GotoBLAS e sto cercando di riprodurre alcuni dei risultati. Questa domanda riguarda il kernel interno di un algoritmo di moltiplicazione della matrice. Nello specifico, sto cercando di calcolare l'C += AB, dove A e B sono matrici 2x2 di tipo double alla velocità di picco della mia CPU. Ci sono due modi per farlo. Un modo è usare le istruzioni SIMD. Il secondo modo è codificare direttamente in assembly usando i registri SIMD.

Quello che ho guardato finora

OpenBLAS Tutti i documenti rilevanti, pagine web del corso, molti molti SO Q & Come affrontare il tema (troppi da elencare), ho compilato sul mio computer, ha esaminato i codici sorgente di OpenBLAS, GotoBLAS e BLIS, i manuali di Agner.

Hardware

mio CPU è un Intel i5 - 540M. È possibile trovare le informazioni CPUID rilevanti su cpu-world.com. La microarchitettura è Nehalem (westmere), quindi può teoricamente calcolare 4 doppi flop di precisione per core per ciclo. Userò solo un core (senza OpenMP), quindi con hyperthreading off e Intel Turbo Boost a 4 fasi, dovrei vedere un picco di (2.533 Ghz + 4*0.133 Ghz) * (4 DP flops/core/cycle) * (1 core) = 12.27 DP Gflops. Per riferimento, con entrambi i core in esecuzione al picco, Intel Turbo Boost offre una velocità in due passaggi e dovrei ottenere un picco teorico di 22.4 DP Gflops.

Setup

dichiaro le mie matrici 2x2 come double e inizializzo con le voci casuali come illustrato nel frammento di codice seguente.

srand(time(NULL)); 
const int n = 2; 
double A[n*n]; 
double B[n*n]; 
double C[n*n]; 
double T[n*n]; 
for(int i = 0; i < n*n; i++){ 
    A[i] = (double) rand()/RAND_MAX; 
    B[i] = (double) rand()/RAND_MAX; 
    C[i] = 0.0; 
} 

ho calcolare una vera risposta utilizzando naive multiplcation matrice-matrice (mostrato di seguito) che mi permette di controllare il risultato visivamente o calcolando la norma L2 di tutti gli elementi

// "true" answer 
for(int i = 0; i < n; i++) 
    for(int j = 0; j < n; j++) 
     for(int k = 0; k < n; k++) 
      T[i*n + j] += A[i*n + k]*B[k*n + j]; 

Per eseguire il codice e ottenere una stima del Gflops, chiamo una volta ogni funzione di moltiplicazione per il riscaldamento e quindi eseguirlo all'interno di un ciclo for per maxiter volte, assicurandosi di azzerare la matrice C ogni volta che sto calcolando C += AB. Il ciclo for si trova all'interno di due istruzioni clock() e viene utilizzato per stimare Gflops. Il frammento di frammento di codice illustra questa parte.

C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 
mult2by2(A,B,C); //warmup 
time1 = clock(); 
for(int i = 0; i < maxiter; i++){ 
     mult2by2(A,B,C); 
     C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 
} 
time2 = clock() - time1; 
time3 = (double)(time2)/CLOCKS_PER_SEC; 
gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; 
mult2by2(A,B,C); // to compute the norm against T 
norm = L2norm(n,C,T); 

codice SIMD

mia CPU supporta vettori a 128 bit, quindi posso ospitare 2 double s in ciascun vettore. Questo è il motivo principale per cui sto facendo una moltiplicazione di matrice 2x2 nel kernel interno. Il codice SIMD calcola un'intera riga di C alla volta.

inline void 
    __attribute__ ((gnu_inline))   
    __attribute__ ((aligned(16))) mult2by2B(  
      const double* restrict A, 
      const double* restrict B, 
      double* restrict C 
     ) 

    { 

    register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; 
    xmm0 = _mm_load_pd(C); 
    xmm1 = _mm_load1_pd(A); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 1); 
    xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C,xmm2); 

    xmm0 = _mm_load_pd(C + 2); 
    xmm1 = _mm_load1_pd(A + 2); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 3); 
    //xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C + 2,xmm2); 
} 

Assmebly (Intel sintassi)

Il mio primo tentativo è stato quello di creare una routine di assembly separato per questa parte e lo chiamano dal main routine. Tuttavia, è stato estremamente lento perché non riesco a incorporare le funzioni extern. Ho scritto l'assemblaggio come assemblaggio in linea come mostrato di seguito. È identico a a quello prodotto da gcc -S -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel. Da quello che ho capito degli schemi microarchitettura di Nehalem, questo processore può eseguire , SSE MUL e SSE MOV in parallelo, il che spiega l'interlacciamento delle istruzioni MUL, ADD, MOV. Noterai che le istruzioni SIMD di cui sopra sono in un ordine diverso perché ho avuto una comprensione diversa dai manuali di Agner Fog. Tuttavia, gcc è intelligente e il codice SIMD sopra riportato viene compilato nell'assieme mostrato nella versione in linea.

inline void 
__attribute__ ((gnu_inline))   
__attribute__ ((aligned(16))) mult2by2A 
    ( 
     const double* restrict A, 
     const double* restrict B, 
     double* restrict C 
    ) 
    { 
    __asm__ __volatile__ 
    (
    "mov  edx, %[A]     \n\t" 
    "mov  ecx, %[B]     \n\t" 
    "mov  eax, %[C]     \n\t" 
    "movapd  xmm3, XMMWORD PTR [ecx]  \n\t" 
    "movapd  xmm2, XMMWORD PTR [ecx+16] \n\t" 
    "movddup xmm1, QWORD PTR [edx]  \n\t" 
    "mulpd  xmm1, xmm3     \n\t" 
    "addpd  xmm1, XMMWORD PTR [eax]  \n\t" 
    "movddup xmm0, QWORD PTR [edx+8]  \n\t" 
    "mulpd  xmm0, xmm2     \n\t" 
    "addpd  xmm0, xmm1     \n\t" 
    "movapd  XMMWORD PTR [eax], xmm0  \n\t" 
    "movddup xmm4, QWORD PTR [edx+16] \n\t" 
    "mulpd  xmm4, xmm3     \n\t" 
    "addpd  xmm4, XMMWORD PTR [eax+16] \n\t" 
    "movddup xmm5, QWORD PTR [edx+24] \n\t" 
    "mulpd  xmm5, xmm2     \n\t" 
    "addpd  xmm5, xmm4     \n\t" 
    "movapd  XMMWORD PTR [eax+16], xmm5 \n\t" 
    : // no outputs 
    : // inputs 
    [A] "m" (A), 
    [B] "m" (B), 
    [C] "m" (C) 
    : //register clobber 
    "memory", 
    "edx","ecx","eax", 
    "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5" 
    ); 
} 

Risultati

compilo il mio codice con i seguenti flag:

gcc -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel 

I risultati per maxiter = 1000000000 sono al di sotto:

********** Inline ASM 
L2 norm: 0.000000e+000, Avg. CPU time: 9.563000, Avg. Gflops: 1.673115 

********** SIMD Version 
L2 norm: 0.000000e+000, Avg. CPU time: 0.359000, Avg. Gflops: 44.568245 

Se forzo la versione SIMD non essere in linea con __attribute__ ((noinline)), i risultati sono:

********** Inline ASM 
L2 norm: 0.000000e+000, Avg. CPU time: 11.155000, Avg. Gflops: 1.434334 

********** SIMD Version 
L2 norm: 0.000000e+000, Avg. CPU time: 11.264000, Avg. Gflops: 1.420455 

Domande

1. Se sia l'ASM in linea e le implementazioni SIMD producono uscita di montaggio identici, perché è la versione di montaggio in modo molto più lento? È come se l'assemblaggio in linea non fosse stato allineato, il che è reso evidente dal secondo set di risultati che mostra prestazioni identiche per ASM "in linea" rispetto a SIMD "noinline". L'unica spiegazione che posso trovare è in Agner Fog Volume 2 pagina 6:

   codice compilato può essere più veloce di codice assembly perché i compilatori possono fare ottimizzazione inter-procedurale e l'ottimizzazione dell'intero programma. Il programmatore dell'assieme di solito deve eseguire funzioni ben definite con un'interfaccia call ben definita che obbedisce a tutte le convenzioni di chiamata per rendere verificabile il codice e verificare . Ciò impedisce molti dei metodi di ottimizzazione utilizzati dai compilatori, come l'inlineamento, l'allocazione del registro, la propagazione costante, l'annullamento comune delle sottoespressioni tramite le funzioni, la pianificazione delle funzioni, ecc. Questi vantaggi possono essere ottenuti utilizzando il codice C++ con funzioni intrinseche del codice assembly .

   Ma l'output dell'assembler per entrambe le versioni è esattamente lo stesso.

2. Perché vedo 44 Gflops nel primo set di risultati?Questo è molto al di sopra del picco di 12 Gflops calcolato, ed è quello che mi aspetterei se eseguissi entrambi i core con calcoli a precisione singola.

EDIT 1 Il commento dice che ci potrebbe essere l'eliminazione del codice morto posso confermare che questo sta accadendo per le istruzioni SIMD. L'output -S mostra che il ciclo for per SIMD ha zeri solo per la matrice C. Posso disabilitarlo disattivando l'ottimizzazione del compilatore con -O0. In tal caso, SIMD esegue 3 volte più lentamente di ASM, ma ASM funziona sempre alla stessa velocità. Anche la norma è ora diversa da zero, ma è ancora OK a 10^-16. Vedo anche che la versione ASM in linea è in linea con i tag APP e NO_APP, ma viene anche srotolata 8 volte nel ciclo for. Penso che srotolare molte volte inciderà pesantemente sulle prestazioni, perché di solito srotolo i loop 4 volte. Nulla di più, secondo la mia esperienza, sembra degradare le prestazioni.

Answer 1

GCC sta ottimizzando via la vostra funzione inline utilizzando intrinseci, mult2by2B, a causa della linea

C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 

Senza di quella linea ci vogliono 2,9 secondi sul computer dal Coliru http://coliru.stacked-crooked.com/a/992304f5f672e257

E con la linea di esso richiede solo 0.000001 http://coliru.stacked-crooked.com/a/9722c39bb6b8590a

Si può anche vedere questo nel gruppo. Se si rilascia il codice di seguito in http://gcc.godbolt.org/, si vedrà che con quella riga di codice salta completamente la funzione.

Tuttavia, quando si esegue l'allineamento del gruppo, GCC NON sta ottimizzando la funzione, mult2by2A, (anche se è in linea). Puoi vederlo anche nell'assemblea.

#include <stdio.h> 
#include <emmintrin.h>     // SSE2 
#include <omp.h> 

inline void 
    __attribute__ ((gnu_inline))   
    __attribute__ ((aligned(16))) mult2by2B(  
      const double* __restrict A, 
      const double* __restrict B, 
      double* __restrict C 
     ) 

    { 

    register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; 
    xmm0 = _mm_load_pd(C); 
    xmm1 = _mm_load1_pd(A); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 1); 
    xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C,xmm2); 

    xmm0 = _mm_load_pd(C + 2); 
    xmm1 = _mm_load1_pd(A + 2); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 3); 
    //xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C + 2,xmm2); 
} 

int main() { 
    double A[4], B[4], C[4]; 
    int maxiter = 10000000; 
    //int maxiter = 1000000000; 
    double dtime; 
    dtime = omp_get_wtime(); 
    for(int i = 0; i < maxiter; i++){ 
     mult2by2B(A,B,C); 
     C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 
    } 
    dtime = omp_get_wtime() - dtime; 
    printf("%f %f %f %f\n", C[0], C[1], C[2], C[3]); 
    //gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; 
    printf("time %f\n", dtime); 
}