使用内联汇编加速矩阵乘法
Using inline assembly to speed up Matrix multiplication
我一直在尝试通过寄存器阻塞、SSE2 向量化和 L1 缓存阻塞来加速矩阵-矩阵乘法 C <- C + alpha * A * B（请注意，我特别选择了转置设置 op(A) =A 和 op(B)=B)。经过一些努力，我编写的代码在单线程模式下 仍然比 GotoBLAS 慢大约 50%。
以下是我在 L1 缓存上进行 "kernel" 方阵-矩阵乘法的代码，在 Goto 的作品中称为 "DGEBB"（一般块-块操作），它将两个 NB*NB 平方相乘矩阵（NB 限制为 4 的倍数）。我检查了它在 GCC 4.8 下的汇编输出，意识到编译器在调度展开的最内层循环：kk-loop 方面做得不好。 希望编译器优化寄存器分配，实现寄存器重用，将计算交错的乘法、加法和内存操作调度为流水线；但是，编译器无法做到这一点。 因此，我想用一些内联汇编 替换最内层的循环。
我对 x86 汇编完全陌生。虽然已经阅读了 GCC 的 extended asm 几个小时，但我仍然不确定如何正确地完成它。我附上了一个愚蠢的版本，我可以尽我所能写，但我知道它是错误的。此版本是从编译器的 kk 循环的原始汇编输出修改而来的。因为我知道如何使用 "movl"、"movapd" 等分配寄存器，所以我按照我喜欢的顺序重新安排了计算。但它还行不通。 1) 在我看来，寄存器 %eax、%ebx、%ecx 在组件内部和外部都使用，这很讨厌。 2）另外，我传递输入和输出操作数的方式不起作用。 3）最后，我真的想要一个可以内联整个kk-loop的版本。如果有人能帮助我，谢谢！
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 , %%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 , %%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，请删除此行和代码中的相应行。
这是一个老问题，可以追溯到我的 HPC Cholesky 分解例程开发的早期阶段。 C 代码已过时，汇编不正确。以后的帖子将遵循此主题。
给出了内联汇编的正确实现。
提供更好的 C 代码。
编写GCC内联汇编时，需要注意状态标志的潜在变化。 给我上了一课
矢量化是 HPC 的关键。 contains some discussion on Intel SSE2/3, while 有一些关于 Intel AVX 的 FMA 指令的信息。
当然这些都只与计算内核有关。还有很多其他工作与如何为最终的高性能 Cholesky 分解例程包装所有内容有关。 第一个 我例程的表现在。
目前我正在升级内核例程以获得更高的性能。可能在这个线程上会有更多的帖子。感谢堆栈溢出社区，特别是 Z boson, Peter Cordes and nominal animal 回答了我的各种问题。在这个过程中我学到了很多，也觉得很开心。 [当然，与此同时，我学会了成为一名更好的 SO 成员。]
我一直在尝试通过寄存器阻塞、SSE2 向量化和 L1 缓存阻塞来加速矩阵-矩阵乘法 C <- C + alpha * A * B（请注意，我特别选择了转置设置 op(A) =A 和 op(B)=B)。经过一些努力，我编写的代码在单线程模式下 仍然比 GotoBLAS 慢大约 50%。
以下是我在 L1 缓存上进行 "kernel" 方阵-矩阵乘法的代码，在 Goto 的作品中称为 "DGEBB"（一般块-块操作），它将两个 NB*NB 平方相乘矩阵（NB 限制为 4 的倍数）。我检查了它在 GCC 4.8 下的汇编输出，意识到编译器在调度展开的最内层循环：kk-loop 方面做得不好。 希望编译器优化寄存器分配，实现寄存器重用，将计算交错的乘法、加法和内存操作调度为流水线；但是，编译器无法做到这一点。 因此，我想用一些内联汇编 替换最内层的循环。
我对 x86 汇编完全陌生。虽然已经阅读了 GCC 的 extended asm 几个小时，但我仍然不确定如何正确地完成它。我附上了一个愚蠢的版本，我可以尽我所能写，但我知道它是错误的。此版本是从编译器的 kk 循环的原始汇编输出修改而来的。因为我知道如何使用 "movl"、"movapd" 等分配寄存器，所以我按照我喜欢的顺序重新安排了计算。但它还行不通。 1) 在我看来，寄存器 %eax、%ebx、%ecx 在组件内部和外部都使用，这很讨厌。 2）另外，我传递输入和输出操作数的方式不起作用。 3）最后，我真的想要一个可以内联整个kk-loop的版本。如果有人能帮助我，谢谢！
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 , %%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 , %%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，请删除此行和代码中的相应行。
这是一个老问题，可以追溯到我的 HPC Cholesky 分解例程开发的早期阶段。 C 代码已过时，汇编不正确。以后的帖子将遵循此主题。
编写GCC内联汇编时，需要注意状态标志的潜在变化。
矢量化是 HPC 的关键。
当然这些都只与计算内核有关。还有很多其他工作与如何为最终的高性能 Cholesky 分解例程包装所有内容有关。 第一个 我例程的表现在
目前我正在升级内核例程以获得更高的性能。可能在这个线程上会有更多的帖子。感谢堆栈溢出社区，特别是 Z boson, Peter Cordes and nominal animal 回答了我的各种问题。在这个过程中我学到了很多，也觉得很开心。 [当然，与此同时，我学会了成为一名更好的 SO 成员。]