向量矩阵乘法、浮点向量、二进制矩阵
Vector matrix multiplication, float vector, binary matrix
我想将大小为 N 的浮点向量与大小为 NxM 的矩阵相乘。
该矩阵是一个二元矩阵（只包含0和1）并且比较稀疏：非零值的密度在1-5%之间。
目前，我正在将其形成为密集向量和稀疏浮点矩阵乘法。
但这只是矫枉过正，不是吗？
如果我将矩阵的列存储为 bitset，然后乘法只是简单地使用 bitset 来索引向量，然后将其求和，会怎么样。
我假设我可以将其形成为 SSE/AVX 中的向量化运算，类似于 load + and + sum 或 load + mask + sum
如果您能指出正确的内在函数来执行此操作，我将不胜感激，主要问题是处理解包位集的最佳方法是什么？
所以你的结果向量的每个元素都是输入向量的掩码和？这些掩码来自矩阵的列，因此它们不是连续的位。
使用连续位图的屏蔽总和对于 AVX512 来说很简单（只需使用合并屏蔽添加或零屏蔽加载）。对于 SSE/AVX2，您将使用 is there an inverse instruction to the movemask instruction in intel avx2? + _mm256_and_ps。或者对跨掩码向量进行优化的一些变体，例如使用 32 位广播负载，然后将其转移到下一步。而不是为每个字节进行另一个未对齐的双字广播。
但是你的掩码位不连续你有一个选择：
- 分别做每个输出向量元素，最后水平求和。需要收集位并制作矢量掩码。可能很难，除了 M=32 的情况，其中位跨度已经将它们与连续的 32 位浮点数对齐。
- 使用 4 或 8 个掩码位的连续组累积 4 或 8 个输出元素的向量。因此，您在外循环上进行矢量化，在内循环中对输入向量进行广播加载。 使用这个。您实际上应该展开多个向量和以隐藏 FP 添加延迟。
像__m256 v = _mm256_set1_ps(invec[i])这样的广播加载基本上是免费的（vbroadcastss是一个纯加载uop，没有ALU shuffle uop）。你不需要任何其他的浮点数洗牌，只需要纯垂直 SIMD，即使在循环结束时也是如此：你只需 _mm256_storeu_ps 进入输出向量。
并且您使用的是连续的掩码位组，因此通常的反向移动掩码问答很有用。
// untested, rough outline of what it might look like
uint8_t matrix[rows * cols]; // bit matrix in chunks of 8 bits
float invec[N], outvec[N]; // A normal function will just take pointer inputs.
constexpr int unroll = 4;
for(int outpos = 0 ; outpos < M-8*unroll+1 ; outpos += 8 * unroll) {
__m256 sum0, sum1, sum2, sum3; //optionally use an array of accumulators, sums[unroll];
sum0 = sum1 = sum2 = sum3 = _mm256_setzero_ps();
// optionally peel the first inner iteration to just load+mask without adding to 0.0
for (int inpos = 0 ; in < N ; in++ ){
__m256 inv = _mm256_set1_ps(invec[inpos]);
__m256 mask0 = inverse_movemask(matrix[outpos*stride + inpos + 0]); // 8 bits -> 8 vector elements
__m256 mask1 = inverse_movemask(matrix[outpos*stride + inpos + 1]);
...
sum0 = _mm256_add_ps(sum0, _mm256_and_ps(inv, mask0) ); // add in[i] or 0.0 according to mask
sum1 = _mm256_add_ps(sum1, _mm256_and_ps(inv, mask1) );
...
}
__m256_storeu_ps(&outvec[outpos + 0*8], sum0);
__m256_storeu_ps(&outvec[outpos + 1*8], sum1);
__m256_storeu_ps(&outvec[outpos + 2*8], sum2);
...
}
not-unrolled __m256 and/or __m128 cleanup for M % (8*unroll) != 0
cleanup for M % 4 != 0 using __m128 broadcast loads
for the last 1..3 rows of masks
maybe use a masked store (AVX2 vmaskmov) or pad your output vector
每个内循环迭代掩码一个浮点数8 * unroll不同的方式，并累积到相应的8 * unroll不同的运行总数。（跨unroll 个向量，每个包含 8 个浮点数。）
这对内存带宽也很有帮助
在 vec*mat 产品中，您只能读取每个位图位一次，但输入向量被有效使用 M 次。遍历连续的位图行提供了良好的局部性，不需要多次加载任何这些缓存行。
即使使用 AVX512 和每个时钟 2x _mm512_mask_add_ps，每个 FP 元素添加 1 位对于位图加载来说带宽也不多。
但是，您循环输入向量 M/(8*unroll) 次。每个和向量的掩码加法使用不同的掩码位，但使用相同的广播输入 float。由于矩阵元素比向量元素小 32 倍，这还不错。
每 4x 或 8x vaddps 指令加载一个浮点数是非常好的计算强度。尤其是在没有 AVX512 的情况下，位图 -> 矢量掩码将花费周期。
为了进一步提高缓存/内存带宽，缓存阻塞/L2 缓存大小 (256kiB) 的循环平铺可能有助于重用输入向量元素.但我不确定您是否可以有效地阻止输入和输出。与 mat*mat 产品不同，只有 O(n^2) 的工作要做。重新读取输入并只编写一个输出流可能没问题，但您可以找到一个中间地带，将部分结果添加到输出向量的部分块中。但是之后不再读取一个连续流中的位矩阵。只要您停在高速缓存行边界，它可能就没问题。
如果您的 NxM 矩阵恰好有 M = 32 那么它与 float 的大小完全匹配，并且 _mm256_loadu_si256 将得到一个具有掩码位的向量for outvec[0] 在每个元素的低位。以及高位中 outvec[31] 的掩码位。您可以使用 _mm256_blendv_ps 将它们应用于求和的输入，然后左移 1 以将下一位移动到顶部位置。 （vblendvps 的替代方法是 psrad 乘以 31 + andps：算术右移以将最高位广播到所有位置）。
但这可能并不比其他方式好多少，即使对于这种特殊情况也是如此。您可以展开不同向量中的多个输出元素，以便多次重复使用浮点向量。
使用 AVX512F，您可以将矩阵行用作 a masked add like _mm512_mask_add_ps 的 __mmask16 值。
sum = _mm512_mask_add_ps(sum, matrix[col*rowstride + row], sum, invec[i]); 如果 matrix 是 uint16_t 的数组。
或使用 AVX512BW，kmovq 64 位掩码进入 k 寄存器，然后 kshift 向下，以匹配 4 个向量累加器的展开。不幸的是，kmov k, [mem] 在 Skylake-X 上是 2 微指令：加载 + 端口 5，而不仅仅是可以写入掩码 regs 的加载微指令。因此，使用 kshift 加载 3x 解包与 4x kmovw k1, [mem] / kmovw k2, [mem+2] 等相比是纯粹的胜利。无法在 [=39= 的底部获取每 16 位掩码数据] 为每个注册时不带 port5 uop。因此它与具有 2 个 FMA 单元的 SKX 内核上的 512 位 FMA/add/mul 吞吐量竞争，否则只是前端吞吐量成本。
我想将大小为 N 的浮点向量与大小为 NxM 的矩阵相乘。
该矩阵是一个二元矩阵（只包含0和1）并且比较稀疏：非零值的密度在1-5%之间。
目前，我正在将其形成为密集向量和稀疏浮点矩阵乘法。
但这只是矫枉过正，不是吗？
如果我将矩阵的列存储为 bitset，然后乘法只是简单地使用 bitset 来索引向量，然后将其求和，会怎么样。
我假设我可以将其形成为 SSE/AVX 中的向量化运算，类似于 load + and + sum 或 load + mask + sum
如果您能指出正确的内在函数来执行此操作，我将不胜感激，主要问题是处理解包位集的最佳方法是什么？
所以你的结果向量的每个元素都是输入向量的掩码和？这些掩码来自矩阵的列，因此它们不是连续的位。
使用连续位图的屏蔽总和对于 AVX512 来说很简单（只需使用合并屏蔽添加或零屏蔽加载）。对于 SSE/AVX2，您将使用 is there an inverse instruction to the movemask instruction in intel avx2? + _mm256_and_ps。或者对跨掩码向量进行优化的一些变体，例如使用 32 位广播负载，然后将其转移到下一步。而不是为每个字节进行另一个未对齐的双字广播。
但是你的掩码位不连续你有一个选择：
- 分别做每个输出向量元素，最后水平求和。需要收集位并制作矢量掩码。可能很难，除了 M=32 的情况，其中位跨度已经将它们与连续的 32 位浮点数对齐。
- 使用 4 或 8 个掩码位的连续组累积 4 或 8 个输出元素的向量。因此，您在外循环上进行矢量化，在内循环中对输入向量进行广播加载。 使用这个。您实际上应该展开多个向量和以隐藏 FP 添加延迟。
像__m256 v = _mm256_set1_ps(invec[i])这样的广播加载基本上是免费的（vbroadcastss是一个纯加载uop，没有ALU shuffle uop）。你不需要任何其他的浮点数洗牌，只需要纯垂直 SIMD，即使在循环结束时也是如此：你只需 _mm256_storeu_ps 进入输出向量。
并且您使用的是连续的掩码位组，因此通常的反向移动掩码问答很有用。
// untested, rough outline of what it might look like
uint8_t matrix[rows * cols]; // bit matrix in chunks of 8 bits
float invec[N], outvec[N]; // A normal function will just take pointer inputs.
constexpr int unroll = 4;
for(int outpos = 0 ; outpos < M-8*unroll+1 ; outpos += 8 * unroll) {
__m256 sum0, sum1, sum2, sum3; //optionally use an array of accumulators, sums[unroll];
sum0 = sum1 = sum2 = sum3 = _mm256_setzero_ps();
// optionally peel the first inner iteration to just load+mask without adding to 0.0
for (int inpos = 0 ; in < N ; in++ ){
__m256 inv = _mm256_set1_ps(invec[inpos]);
__m256 mask0 = inverse_movemask(matrix[outpos*stride + inpos + 0]); // 8 bits -> 8 vector elements
__m256 mask1 = inverse_movemask(matrix[outpos*stride + inpos + 1]);
...
sum0 = _mm256_add_ps(sum0, _mm256_and_ps(inv, mask0) ); // add in[i] or 0.0 according to mask
sum1 = _mm256_add_ps(sum1, _mm256_and_ps(inv, mask1) );
...
}
__m256_storeu_ps(&outvec[outpos + 0*8], sum0);
__m256_storeu_ps(&outvec[outpos + 1*8], sum1);
__m256_storeu_ps(&outvec[outpos + 2*8], sum2);
...
}
not-unrolled __m256 and/or __m128 cleanup for M % (8*unroll) != 0
cleanup for M % 4 != 0 using __m128 broadcast loads
for the last 1..3 rows of masks
maybe use a masked store (AVX2 vmaskmov) or pad your output vector
每个内循环迭代掩码一个浮点数8 * unroll不同的方式，并累积到相应的8 * unroll不同的运行总数。（跨unroll 个向量，每个包含 8 个浮点数。）
这对内存带宽也很有帮助
在 vec*mat 产品中，您只能读取每个位图位一次，但输入向量被有效使用 M 次。遍历连续的位图行提供了良好的局部性，不需要多次加载任何这些缓存行。
即使使用 AVX512 和每个时钟 2x _mm512_mask_add_ps，每个 FP 元素添加 1 位对于位图加载来说带宽也不多。
但是，您循环输入向量 M/(8*unroll) 次。每个和向量的掩码加法使用不同的掩码位，但使用相同的广播输入 float。由于矩阵元素比向量元素小 32 倍，这还不错。
每 4x 或 8x vaddps 指令加载一个浮点数是非常好的计算强度。尤其是在没有 AVX512 的情况下，位图 -> 矢量掩码将花费周期。
为了进一步提高缓存/内存带宽，缓存阻塞/L2 缓存大小 (256kiB) 的循环平铺可能有助于重用输入向量元素.但我不确定您是否可以有效地阻止输入和输出。与 mat*mat 产品不同，只有 O(n^2) 的工作要做。重新读取输入并只编写一个输出流可能没问题，但您可以找到一个中间地带，将部分结果添加到输出向量的部分块中。但是之后不再读取一个连续流中的位矩阵。只要您停在高速缓存行边界，它可能就没问题。
如果您的 NxM 矩阵恰好有 M = 32 那么它与 float 的大小完全匹配，并且 _mm256_loadu_si256 将得到一个具有掩码位的向量for outvec[0] 在每个元素的低位。以及高位中 outvec[31] 的掩码位。您可以使用 _mm256_blendv_ps 将它们应用于求和的输入，然后左移 1 以将下一位移动到顶部位置。 （vblendvps 的替代方法是 psrad 乘以 31 + andps：算术右移以将最高位广播到所有位置）。
但这可能并不比其他方式好多少，即使对于这种特殊情况也是如此。您可以展开不同向量中的多个输出元素，以便多次重复使用浮点向量。
使用 AVX512F，您可以将矩阵行用作 a masked add like _mm512_mask_add_ps 的 __mmask16 值。
sum = _mm512_mask_add_ps(sum, matrix[col*rowstride + row], sum, invec[i]); 如果 matrix 是 uint16_t 的数组。
或使用 AVX512BW，kmovq 64 位掩码进入 k 寄存器，然后 kshift 向下，以匹配 4 个向量累加器的展开。不幸的是，kmov k, [mem] 在 Skylake-X 上是 2 微指令：加载 + 端口 5，而不仅仅是可以写入掩码 regs 的加载微指令。因此，使用 kshift 加载 3x 解包与 4x kmovw k1, [mem] / kmovw k2, [mem+2] 等相比是纯粹的胜利。无法在 [=39= 的底部获取每 16 位掩码数据] 为每个注册时不带 port5 uop。因此它与具有 2 个 FMA 单元的 SKX 内核上的 512 位 FMA/add/mul 吞吐量竞争，否则只是前端吞吐量成本。