更快地归一化下三角矩阵
Normalize lower triangular matrix more quickly
下面的代码似乎不是瓶颈。
我只是想知道是否有更快的方法在 cpu 上使用 SSE4.2 完成此操作。
该代码适用于在 ar_tri 中以以下形式存储为一维数组的矩阵的下三角项：
[ (1,0),
(2,0),(2,1),
(3,0),(3,1),(3,2),
...,
(n,0)...(n,n-1) ]
其中 (x,y) 是矩阵中第 x 行第 y 列的条目。
还有ar_rdia中如下形式的矩阵对角线的平方根倒数(rsqrt):
[ rsqrt(0,0), rsqrt(1,1), ... ,rsqrt(n,n) ]
gcc6.1 -O3 on the Godbolt compiler explorer 使用 SIMD 指令 (mulps) 自动向量化两个版本。三角形版本在每一行的末尾都有清理代码，所以也有一些标量指令。
使用在连续内存中存储为一维数组的矩形矩阵会提高性能吗？
// Triangular version
#include <iostream>
#include <stdlib.h>
#include <stdint.h>
using namespace std;
int main(void){
size_t n = 10000;
size_t n_tri = n*(n-1)/2;
size_t repeat = 10000;
// test 10000 cycles of the code
float* ar_rdia = (float*)aligned_alloc(16, n*sizeof(float));
//reciprocal square root of diagonal
float* ar_triangular = (float*)aligned_alloc(16, n_tri*sizeof(float));
//lower triangular matrix
size_t i,j,k;
float a,b;
k = 0;
for(i = 0; i < n; ++i){
for(j = 0; j < i; ++j){
ar_triangular[k] *= ar_rdia[i]*ar_rdia[j];
++k;
}
}
cout << k;
free((void*)ar_rdia);
free((void*)ar_triangular);
}
// Square version
#include <iostream>
#include <stdlib.h>
#include <stdint.h>
using namespace std;
int main(void){
size_t n = 10000;
size_t n_sq = n*n;
size_t repeat = 10000;
// test 10000 cycles of the code
float* ar_rdia = (float*)aligned_alloc(16, n*sizeof(float));
//reciprocal square root of diagonal
float* ar_square = (float*)aligned_alloc(16, n_sq*sizeof(float));
//lower triangular matrix
size_t i,j,k;
float a,b;
k = 0;
for(i = 0; i < n; ++i){
for(j = 0; j < n; ++j){
ar_square[k] *= ar_rdia[i]*ar_rdia[j];
++k;
}
}
cout << k;
free((void*)ar_rdia);
free((void*)ar_square);
}
汇编输出：
## Triangular version
main:
...
call aligned_alloc
movl , %edi
movq %rax, %rbp
xorl %esi, %esi
xorl %eax, %eax
.L2:
testq %rax, %rax
je .L3
leaq -4(%rax), %rcx
leaq -1(%rax), %r8
movss (%rbx,%rax,4), %xmm0
shrq , %rcx
addq , %rcx
cmpq , %r8
leaq 0(,%rcx,4), %rdx
jbe .L9
movaps %xmm0, %xmm2
leaq 0(%rbp,%rsi,4), %r10
xorl %r8d, %r8d
xorl %r9d, %r9d
shufps [=14=], %xmm2, %xmm2 # broadcast ar_rdia[i]
.L6: # vectorized loop
movaps (%rbx,%r8), %xmm1
addq , %r9
mulps %xmm2, %xmm1
movups (%r10,%r8), %xmm3
mulps %xmm3, %xmm1
movups %xmm1, (%r10,%r8)
addq , %r8
cmpq %rcx, %r9
jb .L6
cmpq %rax, %rdx
leaq (%rsi,%rdx), %rcx
je .L7
.L4: # scalar cleanup
movss (%rbx,%rdx,4), %xmm1
leaq 0(%rbp,%rcx,4), %r8
leaq 1(%rdx), %r9
mulss %xmm0, %xmm1
cmpq %rax, %r9
mulss (%r8), %xmm1
movss %xmm1, (%r8)
leaq 1(%rcx), %r8
jnb .L7
movss (%rbx,%r9,4), %xmm1
leaq 0(%rbp,%r8,4), %r8
mulss %xmm0, %xmm1
addq , %rdx
addq , %rcx
cmpq %rax, %rdx
mulss (%r8), %xmm1
movss %xmm1, (%r8)
jnb .L7
mulss (%rbx,%rdx,4), %xmm0
leaq 0(%rbp,%rcx,4), %rcx
mulss (%rcx), %xmm0
movss %xmm0, (%rcx)
.L7:
addq %rax, %rsi
cmpq 000, %rdi
je .L16
.L3:
addq , %rax
addq , %rdi
jmp .L2
.L9:
movq %rsi, %rcx
xorl %edx, %edx
jmp .L4
.L16:
... print and free
ret
方盒组装的有趣部分：
main:
... allocate both arrays
call aligned_alloc
leaq 40000(%rbx), %rsi
movq %rax, %rbp
movq %rbx, %rcx
movq %rax, %rdx
.L3: # loop over i
movss (%rcx), %xmm2
xorl %eax, %eax
shufps [=15=], %xmm2, %xmm2 # broadcast ar_rdia[i]
.L2: # vectorized loop over j
movaps (%rbx,%rax), %xmm0
mulps %xmm2, %xmm0
movups (%rdx,%rax), %xmm1
mulps %xmm1, %xmm0
movups %xmm0, (%rdx,%rax)
addq , %rax
cmpq 000, %rax
jne .L2
addq , %rcx # no scalar cleanup: gcc noticed that the row length is a multiple of 4 elements
addq 000, %rdx
cmpq %rsi, %rcx
jne .L3
... print and free
ret
存储到三角数组的循环应该可以向量化，但每行末尾的效率低下。根据您发布的 asm，gcc 实际上对两者都进行了自动矢量化。我希望我先看看它，而不是相信你的话，它需要手动矢量化。 :(
.L6: # from the first asm dump.
movaps (%rbx,%r8), %xmm1
addq , %r9
mulps %xmm2, %xmm1
movups (%r10,%r8), %xmm3
mulps %xmm3, %xmm1
movups %xmm1, (%r10,%r8)
addq , %r8
cmpq %rcx, %r9
jb .L6
这看起来与我的手动矢量化版本将编译成的内循环一模一样。 .L4 是对一行中最后最多 3 个元素的完全展开标量清理。 （所以它可能不如我的代码好）。尽管如此，它还是相当不错的，自动矢量化可以让您利用 AVX 和 AVX512 而无需更改源代码。
我编辑了你的问题，将 link 添加到 Godbolt 的代码中，两个版本都作为单独的函数。我没有花时间将它们转换为将数组作为函数参数，因为那样我就必须花时间来正确设置所有 __restrict__ 关键字，并告诉 gcc 数组在 a 上对齐4B * 16 = 64 字节边界，因此它可以根据需要使用对齐加载。
在一行中，您每次都使用相同的 ar_rdia[i]，因此您在行的开头将其广播到向量中一次。然后你只需在源 ar_rdia[j + 0..3] 和目标 ar_triangular[k + 0..3].
之间进行垂直操作
要处理行末尾的最后几个不是矢量大小倍数的元素，我们有两个选择：
- 向量化循环后的标量（或更窄的向量）回退/清理，处理每行的最后最多 3 个元素。
- 将
i 上的循环展开 4，并使用最佳序列处理最后剩下的奇数 0、1、2 和 3 元素一排。因此 j 上的循环将重复 4 次，每次之后都有固定的清理。 这可能是最佳方法。
让最终向量迭代超过一行的末尾，而不是在最后一个完整向量之后停止。所以我们重叠了下一行的开头。由于您的操作不是幂等的，因此此选项效果不佳。此外，确保 k 正确更新为下一行的开头需要一些额外的代码。
不过，这将是可能的，方法是让行的最终向量混合乘数，使当前行末尾以外的元素乘以 1.0（乘法恒等式）。这应该是可行的 blendvpswith a vector of 1.0 to replace some elements of ar_rdia[i] * ar_rdia[j + 0..3]. We'd also have to create a selector mask (maybe by indexing into an array of int32_t row_overshoot_blend_window {0, 0, 0, 0, -1, -1, -1} using j-i as the index, to take a window of 4 elements). Another option is branching to select either no blend or one of three immediate blends (blendps 更快，并且不需要矢量控制掩码，并且分支将有一个容易预测的 table 模式）。
当来自 ar_triangular 的负载与最后一行末尾的存储重叠时，这会导致每 4 行中的 3 行开始时存储转发失败。 IDK 哪个表现最好。
另一个可能更好的选择是执行超过行尾的加载，并使用打包的 SIMD 进行数学计算，然后有条件地存储 1 到 4 个元素。
不读取您分配的内存之外的内容可能需要在缓冲区的末尾留下填充，例如如果最后一行不是 4 个元素的倍数。
/****** Normalize a triangular matrix using SIMD multiplies,
handling the ends of rows with narrower cleanup code *******/
// size_t i,j,k; // don't do this in C++ or C99. Put declarations in the narrowest scope possible. For types without constructors/destructors, it's still a style / human-readability issue
size_t k = 0;
for(size_t i = 0; i < n; ++i){
// maybe put this inside the for() loop and let the compiler hoist it out, to avoid doing it for small rows where the vector loop doesn't even run once.
__m128 vrdia_i = _mm_set1_ps(ar_rdia[i]); // broadcast-load: very efficient with AVX, load+shuffle without. Only done once per row anyway.
size_t j = 0;
for(j = 0; j < (i-3); j+=4){ // vectorize over this loop
__m128 vrdia_j = _mm_loadu_ps(ar_rdia + j);
__m128 scalefac = _mm_mul_ps(vrdia_j, v_rdia_i);
__m128 vtri = _mm_loadu_ps(ar_triangular + k);
__m128 normalized = _mm_mul_ps(scalefac , vtri);
_mm_storeu_ps(ar_triangular + k, normalized);
k += 4;
}
// scalar fallback / cleanup for the ends of rows. Alternative: blend scalefac with 1.0 so it's ok to overlap into the next row.
/* Fine in theory, but gcc likes to make super-bloated code by auto-vectorizing cleanup loops. Besides, we can do better than scalar
for ( ; j < i; ++j ){
ar_triangular[k] *= ar_rdia[i]*ar_rdia[j]; ++k; }
*/
if ((i-j) >= 2) { // load 2 floats (using movsd to zero the upper 64 bits, so mulps doesn't slow down or raise exceptions on denormals or NaNs
__m128 vrdia_j = _mm_castpd_ps( _mm_load_sd(static_cast<const double*>(ar_rdia+j)) );
__m128 scalefac = _mm_mul_ps(vrdia_j, v_rdia_i);
__m128 vtri = _mm_castpd_ps( _mm_load_sd(static_cast<const double*>(ar_triangular + k) ));
__m128 normalized = _mm_mul_ps(scalefac , vtri);
_mm_storel_pi(static_cast<__m64*>(ar_triangular + k), normalized); // movlps. Agner Fog's table indicates that Nehalem decodes this to 2 uops, instead of 1 for movsd. Bizarre!
j+=2;
k+=2;
}
if (j<i) { // last single element
ar_triangular[k] *= ar_rdia[i]*ar_rdia[j];
++k;
//++j; // end of the row anyway. A smart compiler would still optimize it away...
}
// another possibility: load 4 elements and do the math, then movss, movsd, movsd + extractps (_mm_extractmem_ps), or movups to store the last 1, 2, 3, or 4 elements of the row.
// don't use maskmovdqu; it bypasses cache
}
movsd 和 movlps 等同于存储，但不等同于负载。参见 . Update: Agner Fog's insn tables indicate that Nehalem decodes MOVH/LPS/D to 2 fused-domain uops。他们还说 SnB 将其解码为 1，但 IvB 将其解码为 2 uops。那一定是错的。对于 Haswell，他的 table 将内容拆分为 movlps/d（1 个微融合 uop）和 movhps/d（也是 1 个微融合 uop）的单独条目。 movlps 的存储形式是 2 uops 并且任何东西都需要 shuffle port 是没有意义的；它与 movsd 商店的功能完全相同。
如果您的矩阵非常大，不必太担心行尾处理。如果它们很小，那么总时间的更多部分将花费在行的末尾，因此值得尝试多种方法，并仔细查看 asm.
如果源数据是连续的，您可以在此处轻松地即时计算 rsqrt。否则，是的，只将对角线复制到数组中（并在执行该复制时计算 rsqrt，而不是 。在从矩阵的对角线复制到数组时使用标量 rsqrtss 并且没有 NR 步骤，或手动将元素收集到 SIMD 向量中（使用 _mm_set_ps(a[i][i], a[i+1][i+1], a[i+2][i+2], a[i+3][i+3]) 让编译器选择混洗）并执行 rsqrtps + NR 步骤，然后将 4 个结果的向量存储到数组中。
小问题：避免在行尾不做全向量造成的浪费
矩阵的开头是一个特例，因为三个 "ends" 在前 6 个元素中是连续的。 （第 4 行有 4 个元素）。可能值得对此进行特殊处理，并使用两个 SSE 向量来处理前 3 行。或者也许只是前两行在一起，然后第三行作为一个单独的 3 组。实际上，一组 4 和一组 2 是最佳的，因为 SSE 可以做那些 8B 和 16B loads/stores ，但不是 12B。
前 6 个比例因子是 ar_rdia 的前三个元素的乘积，因此我们可以进行单个矢量加载并以多种方式对其进行打乱。
ar_rdia[0]*ar_rdia[0]
ar_rdia[1]*ar_rdia[0], ar_rdia[1]*ar_rdia[1],
ar_rdia[2]*ar_rdia[0], ar_rdia[2]*ar_rdia[1], ar_rdia[2]*ar_rdia[2]
^
end of first vector of 4 elems, start of 2nd.
事实证明，编译器并不擅长发现和利用这里的模式，因此要在这里获得前 10 个元素的最佳代码，我们需要剥离这些迭代并手动优化洗牌和乘法。我决定做前 4 行，因为第 4 行仍然重用 ar_rdia[0..3] 的 SIMD 向量。该向量甚至仍然被第 4 行（第五行）的第一个向量宽度使用。
还值得考虑：做 2、4、4 而不是这个 4、2、4。
void triangular_first_4_rows_manual_shuffle(float *tri, const float *ar_rdia)
{
__m128 vr0 = _mm_load_ps(ar_rdia); // we know ar_rdia is aligned
// elements 0-3 // row 0, row 1, and the first element of row 2
__m128 vi0 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(2, 1, 1, 0));
__m128 vj0 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(0, 1, 0, 0));
__m128 sf0 = vi0 * vj0; // equivalent to _mm_mul_ps(vi0, vj0); // gcc defines __m128 in terms of GNU C vector extensions
__m128 vtri = _mm_load_ps(tri);
vtri *= sf0;
_mm_store_ps(tri, vtri);
tri += 4;
// elements 4 and 5, last two of third row
__m128 vi4 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(3, 3, 2, 2)); // can compile into unpckhps, saving a byte. Well spotted by clang
__m128 vj4 = _mm_movehl_ps(vi0, vi0); // save a mov by reusing a previous shuffle output, instead of a fresh _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(2, 1, 2, 1)); // also saves a code byte (no immediate)
// actually, a movsd from ar_ria+1 would get these two elements with no shuffle. We aren't bottlenecked on load-port uops, so that would be good.
__m128 sf4 = vi4 * vj4;
//sf4 = _mm_movehl_ps(sf4, sf4); // doesn't save anything compared to shuffling before multiplying
// could use movhps to load and store *tri to/from the high half of an xmm reg, but each of those takes a shuffle uop
// so we shuffle the scale-factor down to the low half of a vector instead.
__m128 vtri4 = _mm_castpd_ps(_mm_load_sd((const double*)tri)); // elements 4 and 5
vtri4 *= sf4;
_mm_storel_pi((__m64*)tri, vtri4); // 64bit store. Possibly slower than movsd if Agner's tables are right about movlps, but I doubt it
tri += 2;
// elements 6-9 = row 4, still only needing elements 0-3 of ar_rdia
__m128 vi6 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(3, 3, 3, 3)); // broadcast. clang puts this ahead of earlier shuffles. Maybe we should put this whole block early and load/store this part of tri, too.
//__m128 vi6 = _mm_movehl_ps(vi4, vi4);
__m128 vj6 = vr0; // 3, 2, 1, 0 already in the order we want
__m128 vtri6 = _mm_loadu_ps(tri+6);
vtri6 *= vi6 * vj6;
_mm_storeu_ps(tri+6, vtri6);
tri += 4;
// ... first 4 rows done
}
gcc 和 clang 编译它与 -O3 -march=nehalem 非常相似（启用 SSE4.2 但不启用 AVX）。 See the code on Godbolt, with some other versions that don't compile as nicely:
# gcc 5.3
movaps xmm0, XMMWORD PTR [rsi] # D.26921, MEM[(__v4sf *)ar_rdia_2(D)]
movaps xmm1, xmm0 # tmp108, D.26921
movaps xmm2, xmm0 # tmp111, D.26921
shufps xmm1, xmm0, 148 # tmp108, D.26921,
shufps xmm2, xmm0, 16 # tmp111, D.26921,
mulps xmm2, xmm1 # sf0, tmp108
movhlps xmm1, xmm1 # tmp119, tmp108
mulps xmm2, XMMWORD PTR [rdi] # vtri, MEM[(__v4sf *)tri_5(D)]
movaps XMMWORD PTR [rdi], xmm2 # MEM[(__v4sf *)tri_5(D)], vtri
movaps xmm2, xmm0 # tmp116, D.26921
shufps xmm2, xmm0, 250 # tmp116, D.26921,
mulps xmm1, xmm2 # sf4, tmp116
movsd xmm2, QWORD PTR [rdi+16] # D.26922, MEM[(const double *)tri_5(D) + 16B]
mulps xmm1, xmm2 # vtri4, D.26922
movaps xmm2, xmm0 # tmp126, D.26921
shufps xmm2, xmm0, 255 # tmp126, D.26921,
mulps xmm0, xmm2 # D.26925, tmp126
movlps QWORD PTR [rdi+16], xmm1 #, vtri4
movups xmm1, XMMWORD PTR [rdi+48] # tmp129,
mulps xmm0, xmm1 # vtri6, tmp129
movups XMMWORD PTR [rdi+48], xmm0 #, vtri6
ret
前 4 行总共只有 22 条指令，其中 4 条是 movaps reg-reg moves。 （clang 只用 3 条指令管理，总共有 21 条指令）。我们可能会通过从 ar_rdia+1 将 [ x x 2 1 ] 放入带有 movsd 的向量中来保存一个，而不是另一个 movaps + shuffle。并降低洗牌端口（和一般的 ALU 微处理器）上的压力。
对于 AVX，clang 使用 vpermilps 进行大多数随机播放，但这只会浪费一个字节的代码大小。除非它省电（因为它只有 1 个输入），否则没有理由比 shufps 更喜欢它的直接形式，除非你可以将负载折叠到它里面。
我考虑过使用 palignr 总是一次 4 次通过三角矩阵，但这几乎肯定更糟。你会一直需要那些 palignr，而不仅仅是在最后。
我认为行尾的额外复杂性/更窄 loads/stores 只会让无序执行有事可做。对于大型问题，您将花费大部分时间在内循环中一次执行 16B。这可能会成为内存瓶颈，因此只要无序执行尽可能快地从内存中拉取缓存行，行末的内存密集型工作基本上是免费的。
所以三角矩阵对于这个用例来说仍然很好；保持你的工作集密集和连续内存似乎很好。根据您下一步要做什么，这可能是理想的，也可能不是。
下面的代码似乎不是瓶颈。
我只是想知道是否有更快的方法在 cpu 上使用 SSE4.2 完成此操作。
该代码适用于在 ar_tri 中以以下形式存储为一维数组的矩阵的下三角项：
[ (1,0),
(2,0),(2,1),
(3,0),(3,1),(3,2),
...,
(n,0)...(n,n-1) ]
其中 (x,y) 是矩阵中第 x 行第 y 列的条目。
还有ar_rdia中如下形式的矩阵对角线的平方根倒数(rsqrt):
[ rsqrt(0,0), rsqrt(1,1), ... ,rsqrt(n,n) ]
gcc6.1 -O3 on the Godbolt compiler explorer 使用 SIMD 指令 (mulps) 自动向量化两个版本。三角形版本在每一行的末尾都有清理代码，所以也有一些标量指令。
使用在连续内存中存储为一维数组的矩形矩阵会提高性能吗？
// Triangular version
#include <iostream>
#include <stdlib.h>
#include <stdint.h>
using namespace std;
int main(void){
size_t n = 10000;
size_t n_tri = n*(n-1)/2;
size_t repeat = 10000;
// test 10000 cycles of the code
float* ar_rdia = (float*)aligned_alloc(16, n*sizeof(float));
//reciprocal square root of diagonal
float* ar_triangular = (float*)aligned_alloc(16, n_tri*sizeof(float));
//lower triangular matrix
size_t i,j,k;
float a,b;
k = 0;
for(i = 0; i < n; ++i){
for(j = 0; j < i; ++j){
ar_triangular[k] *= ar_rdia[i]*ar_rdia[j];
++k;
}
}
cout << k;
free((void*)ar_rdia);
free((void*)ar_triangular);
}
// Square version
#include <iostream>
#include <stdlib.h>
#include <stdint.h>
using namespace std;
int main(void){
size_t n = 10000;
size_t n_sq = n*n;
size_t repeat = 10000;
// test 10000 cycles of the code
float* ar_rdia = (float*)aligned_alloc(16, n*sizeof(float));
//reciprocal square root of diagonal
float* ar_square = (float*)aligned_alloc(16, n_sq*sizeof(float));
//lower triangular matrix
size_t i,j,k;
float a,b;
k = 0;
for(i = 0; i < n; ++i){
for(j = 0; j < n; ++j){
ar_square[k] *= ar_rdia[i]*ar_rdia[j];
++k;
}
}
cout << k;
free((void*)ar_rdia);
free((void*)ar_square);
}
汇编输出：
## Triangular version
main:
...
call aligned_alloc
movl , %edi
movq %rax, %rbp
xorl %esi, %esi
xorl %eax, %eax
.L2:
testq %rax, %rax
je .L3
leaq -4(%rax), %rcx
leaq -1(%rax), %r8
movss (%rbx,%rax,4), %xmm0
shrq , %rcx
addq , %rcx
cmpq , %r8
leaq 0(,%rcx,4), %rdx
jbe .L9
movaps %xmm0, %xmm2
leaq 0(%rbp,%rsi,4), %r10
xorl %r8d, %r8d
xorl %r9d, %r9d
shufps [=14=], %xmm2, %xmm2 # broadcast ar_rdia[i]
.L6: # vectorized loop
movaps (%rbx,%r8), %xmm1
addq , %r9
mulps %xmm2, %xmm1
movups (%r10,%r8), %xmm3
mulps %xmm3, %xmm1
movups %xmm1, (%r10,%r8)
addq , %r8
cmpq %rcx, %r9
jb .L6
cmpq %rax, %rdx
leaq (%rsi,%rdx), %rcx
je .L7
.L4: # scalar cleanup
movss (%rbx,%rdx,4), %xmm1
leaq 0(%rbp,%rcx,4), %r8
leaq 1(%rdx), %r9
mulss %xmm0, %xmm1
cmpq %rax, %r9
mulss (%r8), %xmm1
movss %xmm1, (%r8)
leaq 1(%rcx), %r8
jnb .L7
movss (%rbx,%r9,4), %xmm1
leaq 0(%rbp,%r8,4), %r8
mulss %xmm0, %xmm1
addq , %rdx
addq , %rcx
cmpq %rax, %rdx
mulss (%r8), %xmm1
movss %xmm1, (%r8)
jnb .L7
mulss (%rbx,%rdx,4), %xmm0
leaq 0(%rbp,%rcx,4), %rcx
mulss (%rcx), %xmm0
movss %xmm0, (%rcx)
.L7:
addq %rax, %rsi
cmpq 000, %rdi
je .L16
.L3:
addq , %rax
addq , %rdi
jmp .L2
.L9:
movq %rsi, %rcx
xorl %edx, %edx
jmp .L4
.L16:
... print and free
ret
方盒组装的有趣部分：
main:
... allocate both arrays
call aligned_alloc
leaq 40000(%rbx), %rsi
movq %rax, %rbp
movq %rbx, %rcx
movq %rax, %rdx
.L3: # loop over i
movss (%rcx), %xmm2
xorl %eax, %eax
shufps [=15=], %xmm2, %xmm2 # broadcast ar_rdia[i]
.L2: # vectorized loop over j
movaps (%rbx,%rax), %xmm0
mulps %xmm2, %xmm0
movups (%rdx,%rax), %xmm1
mulps %xmm1, %xmm0
movups %xmm0, (%rdx,%rax)
addq , %rax
cmpq 000, %rax
jne .L2
addq , %rcx # no scalar cleanup: gcc noticed that the row length is a multiple of 4 elements
addq 000, %rdx
cmpq %rsi, %rcx
jne .L3
... print and free
ret
存储到三角数组的循环应该可以向量化，但每行末尾的效率低下。根据您发布的 asm，gcc 实际上对两者都进行了自动矢量化。我希望我先看看它，而不是相信你的话，它需要手动矢量化。 :(
.L6: # from the first asm dump.
movaps (%rbx,%r8), %xmm1
addq , %r9
mulps %xmm2, %xmm1
movups (%r10,%r8), %xmm3
mulps %xmm3, %xmm1
movups %xmm1, (%r10,%r8)
addq , %r8
cmpq %rcx, %r9
jb .L6
这看起来与我的手动矢量化版本将编译成的内循环一模一样。 .L4 是对一行中最后最多 3 个元素的完全展开标量清理。 （所以它可能不如我的代码好）。尽管如此，它还是相当不错的，自动矢量化可以让您利用 AVX 和 AVX512 而无需更改源代码。
我编辑了你的问题，将 link 添加到 Godbolt 的代码中，两个版本都作为单独的函数。我没有花时间将它们转换为将数组作为函数参数，因为那样我就必须花时间来正确设置所有 __restrict__ 关键字，并告诉 gcc 数组在 a 上对齐4B * 16 = 64 字节边界，因此它可以根据需要使用对齐加载。
在一行中，您每次都使用相同的 ar_rdia[i]，因此您在行的开头将其广播到向量中一次。然后你只需在源 ar_rdia[j + 0..3] 和目标 ar_triangular[k + 0..3].
要处理行末尾的最后几个不是矢量大小倍数的元素，我们有两个选择：
- 向量化循环后的标量（或更窄的向量）回退/清理，处理每行的最后最多 3 个元素。
- 将
i上的循环展开 4，并使用最佳序列处理最后剩下的奇数 0、1、2 和 3 元素一排。因此j上的循环将重复 4 次，每次之后都有固定的清理。 这可能是最佳方法。 让最终向量迭代超过一行的末尾，而不是在最后一个完整向量之后停止。所以我们重叠了下一行的开头。由于您的操作不是幂等的，因此此选项效果不佳。此外，确保
k正确更新为下一行的开头需要一些额外的代码。不过，这将是可能的，方法是让行的最终向量混合乘数，使当前行末尾以外的元素乘以 1.0（乘法恒等式）。这应该是可行的
blendvpswith a vector of1.0to replace some elements ofar_rdia[i] * ar_rdia[j + 0..3]. We'd also have to create a selector mask (maybe by indexing into an array ofint32_t row_overshoot_blend_window {0, 0, 0, 0, -1, -1, -1}usingj-ias the index, to take a window of 4 elements). Another option is branching to select either no blend or one of three immediate blends (blendps更快，并且不需要矢量控制掩码，并且分支将有一个容易预测的 table 模式）。当来自
ar_triangular的负载与最后一行末尾的存储重叠时，这会导致每 4 行中的 3 行开始时存储转发失败。 IDK 哪个表现最好。另一个可能更好的选择是执行超过行尾的加载，并使用打包的 SIMD 进行数学计算，然后有条件地存储 1 到 4 个元素。
不读取您分配的内存之外的内容可能需要在缓冲区的末尾留下填充，例如如果最后一行不是 4 个元素的倍数。
/****** Normalize a triangular matrix using SIMD multiplies,
handling the ends of rows with narrower cleanup code *******/
// size_t i,j,k; // don't do this in C++ or C99. Put declarations in the narrowest scope possible. For types without constructors/destructors, it's still a style / human-readability issue
size_t k = 0;
for(size_t i = 0; i < n; ++i){
// maybe put this inside the for() loop and let the compiler hoist it out, to avoid doing it for small rows where the vector loop doesn't even run once.
__m128 vrdia_i = _mm_set1_ps(ar_rdia[i]); // broadcast-load: very efficient with AVX, load+shuffle without. Only done once per row anyway.
size_t j = 0;
for(j = 0; j < (i-3); j+=4){ // vectorize over this loop
__m128 vrdia_j = _mm_loadu_ps(ar_rdia + j);
__m128 scalefac = _mm_mul_ps(vrdia_j, v_rdia_i);
__m128 vtri = _mm_loadu_ps(ar_triangular + k);
__m128 normalized = _mm_mul_ps(scalefac , vtri);
_mm_storeu_ps(ar_triangular + k, normalized);
k += 4;
}
// scalar fallback / cleanup for the ends of rows. Alternative: blend scalefac with 1.0 so it's ok to overlap into the next row.
/* Fine in theory, but gcc likes to make super-bloated code by auto-vectorizing cleanup loops. Besides, we can do better than scalar
for ( ; j < i; ++j ){
ar_triangular[k] *= ar_rdia[i]*ar_rdia[j]; ++k; }
*/
if ((i-j) >= 2) { // load 2 floats (using movsd to zero the upper 64 bits, so mulps doesn't slow down or raise exceptions on denormals or NaNs
__m128 vrdia_j = _mm_castpd_ps( _mm_load_sd(static_cast<const double*>(ar_rdia+j)) );
__m128 scalefac = _mm_mul_ps(vrdia_j, v_rdia_i);
__m128 vtri = _mm_castpd_ps( _mm_load_sd(static_cast<const double*>(ar_triangular + k) ));
__m128 normalized = _mm_mul_ps(scalefac , vtri);
_mm_storel_pi(static_cast<__m64*>(ar_triangular + k), normalized); // movlps. Agner Fog's table indicates that Nehalem decodes this to 2 uops, instead of 1 for movsd. Bizarre!
j+=2;
k+=2;
}
if (j<i) { // last single element
ar_triangular[k] *= ar_rdia[i]*ar_rdia[j];
++k;
//++j; // end of the row anyway. A smart compiler would still optimize it away...
}
// another possibility: load 4 elements and do the math, then movss, movsd, movsd + extractps (_mm_extractmem_ps), or movups to store the last 1, 2, 3, or 4 elements of the row.
// don't use maskmovdqu; it bypasses cache
}
movsd 和 movlps 等同于存储，但不等同于负载。参见 MOVH/LPS/D to 2 fused-domain uops。他们还说 SnB 将其解码为 1，但 IvB 将其解码为 2 uops。那一定是错的。对于 Haswell，他的 table 将内容拆分为 movlps/d（1 个微融合 uop）和 movhps/d（也是 1 个微融合 uop）的单独条目。 movlps 的存储形式是 2 uops 并且任何东西都需要 shuffle port 是没有意义的；它与 movsd 商店的功能完全相同。
如果您的矩阵非常大，不必太担心行尾处理。如果它们很小，那么总时间的更多部分将花费在行的末尾，因此值得尝试多种方法，并仔细查看 asm.
如果源数据是连续的，您可以在此处轻松地即时计算 rsqrt。否则，是的，只将对角线复制到数组中（并在执行该复制时计算 rsqrt，而不是 rsqrtss 并且没有 NR 步骤，或手动将元素收集到 SIMD 向量中（使用 _mm_set_ps(a[i][i], a[i+1][i+1], a[i+2][i+2], a[i+3][i+3]) 让编译器选择混洗）并执行 rsqrtps + NR 步骤，然后将 4 个结果的向量存储到数组中。
小问题：避免在行尾不做全向量造成的浪费
矩阵的开头是一个特例，因为三个 "ends" 在前 6 个元素中是连续的。 （第 4 行有 4 个元素）。可能值得对此进行特殊处理，并使用两个 SSE 向量来处理前 3 行。或者也许只是前两行在一起，然后第三行作为一个单独的 3 组。实际上，一组 4 和一组 2 是最佳的，因为 SSE 可以做那些 8B 和 16B loads/stores ，但不是 12B。
前 6 个比例因子是 ar_rdia 的前三个元素的乘积，因此我们可以进行单个矢量加载并以多种方式对其进行打乱。
ar_rdia[0]*ar_rdia[0]
ar_rdia[1]*ar_rdia[0], ar_rdia[1]*ar_rdia[1],
ar_rdia[2]*ar_rdia[0], ar_rdia[2]*ar_rdia[1], ar_rdia[2]*ar_rdia[2]
^
end of first vector of 4 elems, start of 2nd.
事实证明，编译器并不擅长发现和利用这里的模式，因此要在这里获得前 10 个元素的最佳代码，我们需要剥离这些迭代并手动优化洗牌和乘法。我决定做前 4 行，因为第 4 行仍然重用 ar_rdia[0..3] 的 SIMD 向量。该向量甚至仍然被第 4 行（第五行）的第一个向量宽度使用。
还值得考虑：做 2、4、4 而不是这个 4、2、4。
void triangular_first_4_rows_manual_shuffle(float *tri, const float *ar_rdia)
{
__m128 vr0 = _mm_load_ps(ar_rdia); // we know ar_rdia is aligned
// elements 0-3 // row 0, row 1, and the first element of row 2
__m128 vi0 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(2, 1, 1, 0));
__m128 vj0 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(0, 1, 0, 0));
__m128 sf0 = vi0 * vj0; // equivalent to _mm_mul_ps(vi0, vj0); // gcc defines __m128 in terms of GNU C vector extensions
__m128 vtri = _mm_load_ps(tri);
vtri *= sf0;
_mm_store_ps(tri, vtri);
tri += 4;
// elements 4 and 5, last two of third row
__m128 vi4 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(3, 3, 2, 2)); // can compile into unpckhps, saving a byte. Well spotted by clang
__m128 vj4 = _mm_movehl_ps(vi0, vi0); // save a mov by reusing a previous shuffle output, instead of a fresh _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(2, 1, 2, 1)); // also saves a code byte (no immediate)
// actually, a movsd from ar_ria+1 would get these two elements with no shuffle. We aren't bottlenecked on load-port uops, so that would be good.
__m128 sf4 = vi4 * vj4;
//sf4 = _mm_movehl_ps(sf4, sf4); // doesn't save anything compared to shuffling before multiplying
// could use movhps to load and store *tri to/from the high half of an xmm reg, but each of those takes a shuffle uop
// so we shuffle the scale-factor down to the low half of a vector instead.
__m128 vtri4 = _mm_castpd_ps(_mm_load_sd((const double*)tri)); // elements 4 and 5
vtri4 *= sf4;
_mm_storel_pi((__m64*)tri, vtri4); // 64bit store. Possibly slower than movsd if Agner's tables are right about movlps, but I doubt it
tri += 2;
// elements 6-9 = row 4, still only needing elements 0-3 of ar_rdia
__m128 vi6 = _mm_shuffle_ps(vr0, vr0, _MM_SHUFFLE(3, 3, 3, 3)); // broadcast. clang puts this ahead of earlier shuffles. Maybe we should put this whole block early and load/store this part of tri, too.
//__m128 vi6 = _mm_movehl_ps(vi4, vi4);
__m128 vj6 = vr0; // 3, 2, 1, 0 already in the order we want
__m128 vtri6 = _mm_loadu_ps(tri+6);
vtri6 *= vi6 * vj6;
_mm_storeu_ps(tri+6, vtri6);
tri += 4;
// ... first 4 rows done
}
gcc 和 clang 编译它与 -O3 -march=nehalem 非常相似（启用 SSE4.2 但不启用 AVX）。 See the code on Godbolt, with some other versions that don't compile as nicely:
# gcc 5.3
movaps xmm0, XMMWORD PTR [rsi] # D.26921, MEM[(__v4sf *)ar_rdia_2(D)]
movaps xmm1, xmm0 # tmp108, D.26921
movaps xmm2, xmm0 # tmp111, D.26921
shufps xmm1, xmm0, 148 # tmp108, D.26921,
shufps xmm2, xmm0, 16 # tmp111, D.26921,
mulps xmm2, xmm1 # sf0, tmp108
movhlps xmm1, xmm1 # tmp119, tmp108
mulps xmm2, XMMWORD PTR [rdi] # vtri, MEM[(__v4sf *)tri_5(D)]
movaps XMMWORD PTR [rdi], xmm2 # MEM[(__v4sf *)tri_5(D)], vtri
movaps xmm2, xmm0 # tmp116, D.26921
shufps xmm2, xmm0, 250 # tmp116, D.26921,
mulps xmm1, xmm2 # sf4, tmp116
movsd xmm2, QWORD PTR [rdi+16] # D.26922, MEM[(const double *)tri_5(D) + 16B]
mulps xmm1, xmm2 # vtri4, D.26922
movaps xmm2, xmm0 # tmp126, D.26921
shufps xmm2, xmm0, 255 # tmp126, D.26921,
mulps xmm0, xmm2 # D.26925, tmp126
movlps QWORD PTR [rdi+16], xmm1 #, vtri4
movups xmm1, XMMWORD PTR [rdi+48] # tmp129,
mulps xmm0, xmm1 # vtri6, tmp129
movups XMMWORD PTR [rdi+48], xmm0 #, vtri6
ret
前 4 行总共只有 22 条指令，其中 4 条是 movaps reg-reg moves。 （clang 只用 3 条指令管理，总共有 21 条指令）。我们可能会通过从 ar_rdia+1 将 [ x x 2 1 ] 放入带有 movsd 的向量中来保存一个，而不是另一个 movaps + shuffle。并降低洗牌端口（和一般的 ALU 微处理器）上的压力。
对于 AVX，clang 使用 vpermilps 进行大多数随机播放，但这只会浪费一个字节的代码大小。除非它省电（因为它只有 1 个输入），否则没有理由比 shufps 更喜欢它的直接形式，除非你可以将负载折叠到它里面。
我考虑过使用 palignr 总是一次 4 次通过三角矩阵，但这几乎肯定更糟。你会一直需要那些 palignr，而不仅仅是在最后。
我认为行尾的额外复杂性/更窄 loads/stores 只会让无序执行有事可做。对于大型问题，您将花费大部分时间在内循环中一次执行 16B。这可能会成为内存瓶颈，因此只要无序执行尽可能快地从内存中拉取缓存行，行末的内存密集型工作基本上是免费的。
所以三角矩阵对于这个用例来说仍然很好；保持你的工作集密集和连续内存似乎很好。根据您下一步要做什么，这可能是理想的，也可能不是。