使用內聯彙編來加速矩陣乘法

Question

我一直在試圖通過寄存器阻塞，SSE2矢量化和L1緩存阻塞來加速矩陣矩陣乘法C < - C + alpha * A * B（注意我特別選擇了轉置設置op（A）= A和op（B）= B）。經過一番努力，我的書寫代碼仍然是，比單線程模式下的GotoBLAS慢了約50％。

以下是我在L1緩存上進行「內核」方陣矩陣乘法的代碼，在Goto的工作中稱爲「DGEBB」（一般塊塊操作），它將兩個NB * NB方陣（NB限制爲是4的倍數）。我在GCC 4.8下檢查了它的彙編輸出，意識到編譯器在調度展開的最內層循環方面做得不好：kk-loop。 我希望編譯器能夠優化寄存器分配以實現寄存器重用，併爲流水線調度計算交織乘法，加法和內存操作;然而，編譯器沒有做到這一點。由於這個原因，我想用一些內聯彙編替換最裏面的循環。

我完全是x86組裝的新手。雖然已經閱讀了海灣合作委員會的擴展asm幾個小時，但我仍然不確定如何正確地做到這一點。我附加了一個我可以寫得最好的愚蠢版本，但知道它是錯誤的。該版本從編譯器的kk循環的原始程序集輸出中修改。因爲我知道如何使用「movl」，「movapd」等來分配寄存器，所以我按照我想象的順序重新安排了計算。但它還沒有工作。 1）在我看來，寄存器％eax，％ebx，％ecx在組件內部和外部都被使用，這是討厭的。 2）另外，我通過輸入和輸出操作數的方式不起作用。 3）最後，我真的想要一個可以內聯整個kk循環的版本。謝謝如果有人能幫助我！

爲DGEBB C代碼（稱爲DGEBB_SSE2_x86，如我的筆記本電腦是32位x86機器，具有SSE2 - SSE4.1支持）：

#include <stdint.h> /* type define of "uintptr_t" */ 
#include <emmintrin.h> /* double precision computation support since SSE2 */ 
#include <R.h> /* use R's error handling error() */ 

void DGEBB_SSE2_x86 (int *NB, double *ALPHA, double *A, double *B, double *C) { 
/* check "nb", must be a multiple of 4 */ 
int TWO=2, FOUR=4, nb=*NB; if (nb%FOUR) error("error in DGEBB_SSE2_x86: nb is not a multiple of 4!\n"); 
/* check memory alignment of A, B, C, 16 Byte alignment is mandatory (as XMM registers are 128-bit in length) */ 
uintptr_t sixteen_bytes=0xF; 
if ((uintptr_t)A & sixteen_bytes) error("error in DGEBB_SSE2_x86: A is not 16 Bytes aligned in memory!"); 
if ((uintptr_t)B & sixteen_bytes) error("error in DGEBB_SSE2_x86: B is not 16 Bytes aligned in memory!"); 
if ((uintptr_t)C & sixteen_bytes) error("error in DGEBB_SSE2_x86: C is not 16 Bytes aligned in memory!"); 
/* define vector variables */ 
__m128d C1_vec_reg=_mm_setzero_pd(), C2_vec_reg=C1_vec_reg, C3_vec_reg=C1_vec_reg, C4_vec_reg=C1_vec_reg,A1_vec_reg, A2_vec_reg, B_vec_reg, U_vec_reg; 
/* define scalar variables */ 
int jj, kk, ii, nb2=nb+nb, nb_half=nb/TWO; 
double *B1_copy, *B1, *C1, *a, *b, *c, *c0; 
/* start triple loop nest */ 
C1=C;B1=B; /* initial column tile of C and B */ 
jj=nb_half; 
while (jj--) { 
    c=C1;B1_copy=B1;C1+=nb2;B1+=nb2;b=B1_copy; 
    for (ii=0; ii<nb; ii+=FOUR) { 
    a=A+ii;b=B1_copy; 
    kk=nb_half; 
    while (kk--) { 
    /* [kernel] amortize pointer arithmetic! */ 
    A1_vec_reg=_mm_load_pd(a); /* [fetch] */ 
    B_vec_reg=_mm_load1_pd(b); /* [fetch] */ 
    U_vec_reg=_mm_mul_pd(A1_vec_reg,B_vec_reg);C1_vec_reg=_mm_add_pd(C1_vec_reg,U_vec_reg); /* [daxpy] */ 
    A2_vec_reg=_mm_load_pd(a+TWO);a+=nb; /* [fetch] */ 
    U_vec_reg=_mm_mul_pd(A2_vec_reg,B_vec_reg);C2_vec_reg=_mm_add_pd(C2_vec_reg,U_vec_reg); /* [daxpy] */ 
    B_vec_reg=_mm_load1_pd(b+nb);b++; /* [fetch] */ 
    U_vec_reg=_mm_mul_pd(A1_vec_reg,B_vec_reg);C3_vec_reg=_mm_add_pd(C3_vec_reg,U_vec_reg); /* [daxpy] */ 
    A1_vec_reg=_mm_load_pd(a); /* [fetch] */ 
    U_vec_reg=_mm_mul_pd(A2_vec_reg,B_vec_reg);C4_vec_reg=_mm_add_pd(C4_vec_reg,U_vec_reg); /* [daxpy]*/ 
    B_vec_reg=_mm_load1_pd(b); /* [fetch] */ 
    U_vec_reg=_mm_mul_pd(A1_vec_reg,B_vec_reg);C1_vec_reg=_mm_add_pd(C1_vec_reg,U_vec_reg); /* [daxpy] */ 
    A2_vec_reg=_mm_load_pd(a+TWO);a+=nb; /* [fetch] */ 
    U_vec_reg=_mm_mul_pd(A2_vec_reg,B_vec_reg);C2_vec_reg=_mm_add_pd(C2_vec_reg,U_vec_reg); /* [daxpy] */ 
    B_vec_reg=_mm_load1_pd(b+nb);b++; /* [fetch] */ 
    U_vec_reg=_mm_mul_pd(A1_vec_reg,B_vec_reg);C3_vec_reg=_mm_add_pd(C3_vec_reg,U_vec_reg); /* [daxpy] */ 
    U_vec_reg=_mm_mul_pd(A2_vec_reg,B_vec_reg);C4_vec_reg=_mm_add_pd(C4_vec_reg,U_vec_reg); /* [daxpy] */ 
    } /* [end of kk-loop] */ 
    /* [write-back] amortize pointer arithmetic! */ 
    A2_vec_reg=_mm_load1_pd(ALPHA); 
    U_vec_reg=_mm_load_pd(c);c0=c+nb;C1_vec_reg=_mm_mul_pd(C1_vec_reg,A2_vec_reg); /* [fetch] */ 
    A1_vec_reg=U_vec_reg;C1_vec_reg=_mm_add_pd(C1_vec_reg,A1_vec_reg);U_vec_reg=_mm_load_pd(c0); /* [fetch] */ 
    C3_vec_reg=_mm_mul_pd(C3_vec_reg,A2_vec_reg);_mm_store_pd(c,C1_vec_reg);c+=TWO; /* [store] */ 
    A1_vec_reg=U_vec_reg;C3_vec_reg=_mm_add_pd(C3_vec_reg,A1_vec_reg);U_vec_reg=_mm_load_pd(c); /* [fetch] */ 
    C2_vec_reg=_mm_mul_pd(C2_vec_reg,A2_vec_reg);_mm_store_pd(c0,C3_vec_reg);c0+=TWO; /* [store] */ 
    A1_vec_reg=U_vec_reg;C2_vec_reg=_mm_add_pd(C2_vec_reg,A1_vec_reg);U_vec_reg=_mm_load_pd(c0); /* [fetch] */ 
    C4_vec_reg=_mm_mul_pd(C4_vec_reg,A2_vec_reg);_mm_store_pd(c,C2_vec_reg);c+=TWO; /* [store] */ 
    C4_vec_reg=_mm_add_pd(C4_vec_reg,U_vec_reg);_mm_store_pd(c0,C4_vec_reg); /* [store] */ 
    C1_vec_reg=_mm_setzero_pd();C3_vec_reg=C1_vec_reg;C2_vec_reg=C1_vec_reg;C4_vec_reg=C1_vec_reg; 
    } /* [end of ii-loop] */ 
} /* [end of jj-loop] */ 
} 

我笨的KK-環內聯組件的版本是在這裏：

 while (kk--) { 
    asm("movapd %0, %%xmm3\n\t"  /* C1_vec_reg -> xmm3 */ 
     "movapd %1, %%xmm1\n\t"  /* C2_vec_reg -> xmm1 */ 
     "movapd %2, %%xmm2\n\t"  /* C3_vec_reg -> xmm2 */ 
     "movapd %3, %%xmm0\n\t"  /* C4_vec_reg -> xmm0 */ 
     "movl %4, %%eax\n\t" /* pointer a -> %eax */ 
     "movl %5, %%edx\n\t" /* pointer b -> %edx */ 
     "movl %6, %%ecx\n\t" /* block size nb -> %ecx */ 
     "movapd (%%eax), %%xmm5\n\t" /* A1_vec_reg -> xmm5 */ 
    "movsd (%%edx), %%xmm4\n\t"  /* B_vec_reg -> xmm4 */ 
    "unpcklpd %%xmm4, %%xmm4\n\t" 
     "movapd %%xmm5, %%xmm6\n\t"  /* xmm5 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
    "addpd %%xmm6, %%xmm3\n\t"  /* xmm3 += xmm6 */ 
     "movapd 16(%%eax), %%xmm7\n\t"  /* A2_vec_reg -> xmm7 */ 
     "movapd %%xmm7, %%xmm6\n\t"  /* xmm7 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
    "addpd %%xmm6, %%xmm1\n\t"  /* xmm1 += xmm6 */ 
    "movsd (%%edx,%%ecx), %%xmm4\n\t"  /* B_vec_reg -> xmm4 */ 
    "addl $8, %%edx\n\t"   /* b++ */ 
    "movsd (%%edx), %%xmm4\n\t"  /* B_vec_reg -> xmm4 */ 
    "unpcklpd %%xmm4, %%xmm4\n\t" 
     "movapd %%xmm5, %%xmm6\n\t"  /* xmm5 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
    "addpd %%xmm6, %%xmm2\n\t"  /* xmm2 += xmm6 */ 
    "addl %%ecx, %%eax\n\t"   /* a+=nb */ 
    "movapd (%%eax), %%xmm5\n\t"  /* A1_vec_reg -> xmm5 */ 
     "movapd %%xmm7, %%xmm6\n\t"  /* xmm7 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
    "addpd %%xmm6, %%xmm0\n\t"  /* xmm0 += xmm6 */ 
    "movsd (%%edx), %%xmm4\n\t"  /* B_vec_reg -> xmm4 */ 
    "unpcklpd %%xmm4, %%xmm4\n\t" 
     "movapd %%xmm5, %%xmm6\n\t"  /* xmm5 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
     "addpd %%xmm6, %%xmm3\n\t"  /* xmm3 += xmm6 */ 
     "movapd 16(%%eax), %%xmm7\n\t"  /* A2_vec_reg -> xmm7 */ 
     "movapd %%xmm7, %%xmm6\n\t"  /* xmm7 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
    "addpd %%xmm6, %%xmm1\n\t"  /* xmm1 += xmm6 */ 
    "movsd (%%edx,%%ecx), %%xmm4\n\t"  /* B_vec_reg -> xmm4 */ 
    "addl $8, %%edx\n\t"   /* b++ */ 
    "movsd (%%edx), %%xmm4\n\t"  /* B_vec_reg -> xmm4 */ 
    "unpcklpd %%xmm4, %%xmm4\n\t" 
     "movapd %%xmm5, %%xmm6\n\t"  /* xmm5 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
    "addpd %%xmm6, %%xmm2\n\t"  /* xmm2 += xmm6 */ 
     "movapd %%xmm7, %%xmm6\n\t"  /* xmm7 -> xmm6 */ 
     "mulpd %%xmm4, %%xmm6\n\t"  /* xmm6 *= xmm4 */ 
    "addpd %%xmm6, %%xmm0\n\t"  /* xmm0 += xmm6 */ 
    "addl %%ecx, %%eax" 
     : "+x"(C1_vec_reg), "+x"(C2_vec_reg), "+x"(C3_vec_reg), "+x"(C4_vec_reg), "+m"(a), "+m"(b) 
     : "x"(C1_vec_reg), "x"(C2_vec_reg), "x"(C3_vec_reg), "x"(C4_vec_reg), "4"(a), "5"(b), "rm"(nb)); 
} 

下面是代碼的一些解釋：

Unrolling out loops to expose a micro "dger" kernel for register resue: 
(c11 c12) += (a1) * (b1 b2) 
(c21 c22) (a2) 
(c31 c32) (a3) 
(c41 c42) (a4) 
This can be implemented as 4 vectorized "daxpy": 
(c11) += (a1) * (b1) , (c31) += (a3) * (b1) , (c12) += (a1) * (b2) , (c32) += (a3) * (b2) . 
(c21) (a2) (b1)  (c41) (a4) (b1)  (c22) (a2) (b2)  (c42) (a4) (b2) 
4 micor C-vectors are held constantly in XMM registers named C1_vec_reg, C2_vec_reg, C3_vec_reg, C4_vec_reg. 
2 micro A-vectors are loaded into XMM registers named A1_vec_reg, A2_vec_reg. 
2 micro B-vectors can reuse a single XMM register named B_vec_reg. 
1 additional XMM register, U_vec_reg, will store temporary values. 
The above scheduling exploits all 8 XMM registers on x84 architectures with SIMD unit, and each XMM is used twice after loaded. 

PS：我是從統計的R用戶組。頭文件允許使用R的錯誤處理功能error（）。這隻會終止C程序而不是整個R程序。如果您不使用R，請刪除代碼中的此行和相應的行。

Answer 1

這是一個老問題，回到HPC Cholesky分解例程開發的早期階段。 C代碼已過時，程序集天真不正確。稍後的帖子按照這個線程。

(inline assembly in C) Assembler messages: Error: unknown pseudo-op:給出了內聯彙編的正確實現。

How to ask GCC to completely unroll this loop (i.e., peel this loop)?給出更好的C代碼。

在編寫GCC內聯彙編時，需要注意狀態標誌的潛在變化。 (inline assembly in C) Funny memory segmentation fault對我來說是一個教訓。

矢量化是HPC的關鍵。 SSE instruction MOVSD (extended: floating point scalar & vector operations on x86, x86-64)包含關於英特爾SSE2/3的一些討論，而FMA instruction _mm256_fmadd_pd(): "132", "231" and "213"?包含有關英特爾AVX的FMA指令的一些信息。

當然，所有這些只與計算內核有關。還有很多其他的工作涉及到最終的高性能Cholesky分解例程的所有方面。 第一個我的例程發佈的性能是在Why can't my CPU maintain peak performance in HPC。

目前我正在升級內核例程以獲得更高的性能。可能在這個線程上會有更多的帖子。感謝堆棧溢出社區，特別是Z boson,Peter Cordes和nominal animal回答了我的各種問題。在這個過程中，我學到了很多，並感到非常高興。 [當然，同時，我學會了成爲更好的SO成員。]