二进制矩阵向量乘法的内在函数

Posted 2023-02-16

【中文标题】二进制矩阵向量乘法的内在函数

【英文标题】：Intrinsics for binary matrix vector multiplication

【发布时间】：2019-01-03 19:12:14

【问题描述】：

我正在尝试在二进制字段上实现矩阵向量乘法。向量 x 的维度为 1xa，矩阵 M 的维度为 axb，结果 y = a * M 的大小为 1xb。现在，我实现了它，使得 x 和 M 的类型为 uint8_t*，即，我连接 M 的列，因为它们也被连续访问。函数如下：

void mul(uint8_t M, size_t a, size_t b, uint8_t* x, uint8_t* y) 
    uint8_t val;
    uint8_t *ptr;
    for(size_t i = 0; i < b; i++) 
        val = 0;
        ptr = M + i * a;
        for(size_t j = 0; j < a; j++) 
            val ^= (x[j] & *ptr++);
        
        y[i] = bit;
    

M 和 x 已被调用者分配为

M = malloc(sizeof(uint8_t) * a * b);
x = malloc(sizeof(uint8_t) * a);
y = malloc(sizeof(uint8_t) * b);

由于这个例程被调用了数十亿次，我需要优化它；）为此，我正在考虑

不是将每个 0/1 表示为单独的 uint8_t（即 8 位），我可以将“x”和“M”中的所有位打包到更小尺寸的 uint64_t 数组中，例如 ap 和 Mp，其中

ap = (size_t) ceil ((double) a / 64); mp = (size_t) ceil ((double) (a*b) / 64);

使用向量内在函数。

到目前为止，我完成了 M 的（左对齐）打包（正确对齐）和乘法

typedef uint64_t word;
#define WORD_BITS      (CHAR_BIT * sizeof (word))

void mul_fast(word *M, size_t Mlen, word *x, size_t xlen, size_t b, word *y) 

    for(size_t i = 0; i < Mlen; i++) 
        y[i/xlen] ^= (M[i] & x[i % xlen]);
    
    for(size_t i = 0; i < b; i++) 
        y[i] = __builtin_popcountll(y[i]) & 1;
    

然而，事实证明，上面的方法比直接实现 mul() 慢得多。

您有什么想法或参考吗？我不是汇编专家，所以比较 gcc -S 的输出并不能告诉我太多：/

谢谢你，最好的问候，汤姆。

【问题讨论】：

我不是超级专家，但每当我在应该过度优化的代码中看到% 或/ 时，我都会怀疑。除法和模数非常慢，据我所知，在比语言更基础的层面上。

@kabanus：如果ulen 是 2 的编译时常数幂，那很好。 （size_t 是无符号类型，所以它实际上只是移位和与。有符号除法/余数与移位和掩码具有不同的舍入语义，因此需要一些额外的指令）。但如果它是一个运行时变量，那么是的，你真的搞砸了。

@PeterCordes 我明白了，谢谢你的启发，我什至没有想到这一点，很好。

相关：Large (0,1) matrix multiplication using bitwise AND and popcount instead of actual int or float multiplies? 可能是重复的，我认为它在做完全相同的问题。使用 SSE2 或 AVX2 将字节打包为位图应使用 _mm_movemask_epi8 作为构建块。

感谢您的参考！我想我了解如何使用 _mm_movemask_epi8 从相应的 uint8_t 数组创建位图。但是，我仍然坚持实际的矩阵向量乘法：与其他解决方案相反，我更愿意将我的矩阵简单地作为位图的一维向量。这不应该更有效吗？此外，二进制矩阵不是那么大，例如通常为 200x1000 位，即 ceil(200*1000/64 ) = 3125 uin64_t 整数。你能看出上面的代码有什么问题吗？

【参考方案1】：

汇编输出的相关区别是：

.L26: - movq %r10, %rax - xorl %edx, %edx - divq %rcx - movq (%r11,%rdx,8), %rdx - andq (%rdi,%r10,8), %rdx - addq $1, %r10 - xorq %rdx, (%r9,%rax,8) - cmpq %r10, %rsi + movq %rax, %rcx + movq %rax, %r10 + andl $1, %ecx + shrq %r10 + movq (%rdx,%rcx,8), %rcx + andq (%rdi,%rax,8), %rcx + addq $1, %rax + xorq %rcx, (%r9,%r10,8) + cmpq %rax, %rsi 你能看出罪魁祸首是什么吗？

【讨论】：

将 div 替换为编译时常量 xlen 中的 shift/AND 使其速度提高约 30 到 90 倍。整数除法是一个大量瓶颈，尤其是对于 64 位整数（而不是 32 位）。 C++ code for testing the Collatz conjecture faster than hand-written assembly - why?/.