logo标签列表

使用位运算和popcount而不是实际的整数或浮点数乘法进行大型（0,1）矩阵乘法？

c++ssematrix-multiplicationavxbitset
10
对于大型二进制矩阵（10Kx20K）的乘法，我通常会将矩阵转换为浮点数并执行浮点矩阵乘法，因为整数矩阵乘法速度相对较慢（可以看看这里）。
然而这一次，我需要执行数十万次这样的乘法运算，平均性能提升每毫秒都很重要。
我希望结果是一个 int 或 float 矩阵，因为乘积可能具有不为 0 或 1 的元素。输入矩阵元素都为 0 或 1，因此它们可以存储为单个位。
在行向量和列向量之间进行内积（以产生输出矩阵的一个元素）时，乘法简化为按位 AND。加法仍然是加法，但我们可以使用人口统计函数来添加位，而不必逐个循环。
一些其他的布尔值 / 二进制矩阵函数是通过 OR 每个位来产生位矩阵结果，但这不是我所需要的。
以下是一个示例代码，展示了将问题表示为 std::bitset、AND 和 count 操作比矩阵乘法更快。
#include <iostream>
using std::cout; using std::endl;
#include <vector>
    using std::vector;
#include <chrono>
#include <Eigen/Dense>
    using Eigen::Map; using Eigen::Matrix; using Eigen::MatrixXf;
#include <random>
    using std::random_device; using std::mt19937; using std::uniform_int_distribution;
#include <bitset>
    using std::bitset;

using std::floor;

const int NROW = 1000;
const int NCOL = 20000;

const float DENSITY = 0.4;
const float DENOMINATOR = 10.0 - (10*DENSITY);

void fill_random(vector<float>& vec) {
    random_device rd;
    mt19937 eng(rd());
    uniform_int_distribution<> distr(0, 10);
    int nnz = 0;
    for (int i = 0; i < NROW*NCOL; ++i)
        vec.push_back(floor(distr(eng)/DENOMINATOR));
}

void matmul(vector<float>& vec){
    float *p = vec.data();
    MatrixXf A = Eigen::Map<Eigen::Matrix<float, NROW, NCOL, Eigen::RowMajor>>(p);
    cout << "Eigen matrix has " << A.rows() << " rows and " << A.cols() << " columns." << endl;
    cout << "Total non-zero values : " << A.sum() << endl;
    cout << "The density of non-zero values is " <<  A.sum() * 1.0 / (A.cols()*A.rows()) << endl;

    auto start = std::chrono::steady_clock::now();
    MatrixXf B = A.transpose() * A;
    auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
    cout << "Mat mul took " << end << " ms"<< endl;

    // Just to make sure the operation is not skipped by compiler
    cout << "Eigen coo ";
    for (int i=0; i<10; ++i)
        cout << B(0,i) << " ";
    cout << endl;
}


void bitset_op(vector<float>& vec) {
    // yeah it's not a great idea to set size at compile time but have to
    vector<bitset<NROW>> col_major(NCOL);

    // right, multiple par for isn't a good idea, maybe in a parallel block
    // Doing this for simplicity to profile second loop timing 
    // converting row major float vec to col major bool vec
    #pragma omp parallel for
    for (int j=0; j < NCOL; ++j) {
        for (int i=0; i < NROW; ++i) {
            col_major[j].set(i, vec[i*NCOL + j] && 1);
        }
    }

    auto start = std::chrono::steady_clock::now();
    vector<int> coo;
    coo.assign(NCOL*NCOL, 0);
    #pragma omp parallel for
    for (int j=0; j < NCOL; ++j) {
        for (int k=0; k<NCOL; ++k) {
            coo[j*NCOL + k] = (col_major[j]&col_major[k]).count();
        }
    }
    auto end = (std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start)).count();
    cout << "bitset intersection took " << end << " ms"<< endl;

    // Just to make sure the operation is not skipped by compiler
    cout << "biset coo ";
    for (int i=0; i<10; ++i)
        cout << coo[i] << " ";
    cout << endl;
}


int main() {
    // Saving to float instead of int to speed up matmul
    vector<float> vec;
    fill_random(vec);
    matmul(vec);
    bitset_op(vec);
}

使用以下方式运行：

g++ -O3 -fopenmp -march=native -I. -std=c++11 code.cpp -o code

我理解为：

Eigen matrix has 1000 rows and 20000 columns.
Total non-zero values : 9.08978e+06
The density of non-zero values is 0.454489
Mat mul took 1849 ms
Eigen coo 458 206 208 201 224 205 204 199 217 210
bitset intersection took 602 ms
biset coo 458 206 208 201 224 205 204 199 217 210

如您所见，将matmul作为一组位集操作比Eigen的浮点matmul快约3倍，这是有道理的。

我想强调的是，我需要在100K上执行此操作（在HPC或云中），平均毫秒性能提高会有所不同。

我没有绑定到任何特定的库、C++标准等。因此，请随意提供您认为比使用GPU更快的任何解决方案，因为我由于多种原因无法使用它。

1
我认为如果编译器没有使用这些功能的话，你可以使用SSE和POPCNT来制作一个明显更快的版本。 - geza
1
你是否有可用的AVX2（Intel Haswell或更高版本）？我假设你使用的是Intel，因为这在当今HPC /云计算领域几乎是标准，但如果你使用AMD，请告诉我们。在Intel上，使用AVX2 vpshufb方法（4位半字节的LUT）进行大型数组的弹出计数比使用64位popcnt更快。 - Peter Cordes
1
希望你的编译器在使用-march=native选项编译std::bitset.count()时能够选择最佳策略。@geze：-march=native会启用支持它的CPU上的-mpopcnt。此外，gcc的std::bitset<64>确实使用了popcnt。 - Peter Cordes
@PeterCordes 是的，我确实有可用的AVX2。我主要使用Google云，获取更新的架构也很容易。 - NULL
@geza -mpopcnt 确实已启用 - NULL
1个回答
17
我不确定编译器在没有手动向量化（使用内置函数或C++向量类包装器，如Agner Fog's VCL，如果您的项目许可证与GPL兼容）的情况下能为您做多少工作。也有一些非GPL的包装器。
缓存阻塞矩阵乘法是一门艺术（并且在这里非常重要），如果您可以使用Eigen的现有模板，但使用一个不同的类，该类在整数上使用按位and，而不是浮点数上的乘法，则非常好。我不确定是否可能实现这一点。
我查找了一些资料，大部分关于二进制矩阵的文献都是关于产生布尔结果的（包括像这样的SO问题）。向量内积是使用AND作为乘法，但使用XOR或OR作为加法，而不是popcount。也许有一个我错过的搜索术语，描述了“正常”的矩阵，它们只是恰好是（0,1）矩阵，但乘积不会是。

由于每一毫秒都很重要，您可能需要手动进行矢量化。

这并不是说向量整数操作在一般情况下很慢，只是向量整数乘法很慢，特别是与最近的x86硬件（尤其是Intel）相比，向量浮点FMA每个时钟周期可处理2个256位的向量。

既然你不需要实际乘以布尔元素，只需要AND（每个时钟周期3个向量的吞吐量），那对你来说就不是问题了。做更多元素的向量应该会带来更高的效率收益，从而弥补每个向量的额外成本。

当然，这假定一个整数矩阵乘法实现使用与等效FP矩阵乘法相同的缓存块和其他优化，并且如果您不想（或不知道如何）自己编写它，并且找不到可以为您执行此操作的库，则存在问题。

我只是回答如何通过最佳实现使其有效率的问题。 标题问题的答案非常肯定：使用实际乘法，特别是32位元素，是一种极大的时间浪费。

存储格式选项：

每个字节一个元素（0/1）：

  • float 的密度是 4 倍（缓存占用 / 内存带宽 / 向量中的元素数）
  • 易于通过字节洗牌进行转置
  • 如果需要，垂直加法很容易（例如，对多行或多列进行矢量化。如果您的数据以使此操作无需额外洗牌的方式交错，则可以避免在最后进行水平求和。）

每个字节 4 个元素，打包到低四位：

  • 比单独字节密度高4倍
  • 使用AVX2 vpshufb 进行popcount非常高效。在L1D缓存中输入热点的情况下，理论上可以以每个时钟周期（每个核心）128个AND结果元素的吞吐量进行加载/AND/累加-popcount。每个时钟周期的4个融合域uop饱和了SKL/HSW前端问题带宽，每个时钟周期发出4个，不会在3个向量ALU端口上形成瓶颈，因为其中一个uop是纯负载。（另一个负载微粒与vpand融合）。当在L2带宽上形成瓶颈（每个周期约为32B负载），以64个元素每个时钟周期运行。请参见下文。
  • 从整数或打包位图创建较慢（但如果您按交错顺序将位放入向量以有效地打包/解包到有序字节中，而不是强制位按顺序排列，则不会很差）。
  • 难以转置（可能比完全打包更差）

打包位：

  • 每个AVX2向量具有256个元素，比单独字节密度高8倍。
  • 可以使用pmovmskb从非交错存储顺序的向量创建。(对于即时创建来说并不是很有用，因为这会将结果放置在整数寄存器而不是向量中。交错位顺序可能是最佳选择，特别是在转置期间进行解包)。
  • 使用AVX2进行popcount相当有效：掩码/移位+掩码/2xvpshufb。(9个融合域uops(8个向量ALU uops)用于AND+累加256个元素的popcount(来自2个行/列向量)，而4个字节策略(来自4个行/列向量)需要8个uops(6个向量ALU uops))。ALU端口瓶颈将其限制为每个时钟周期从L1D或L2读取96个元素。因此，在理论上，当它在L2带宽上出现瓶颈时，其内积吞吐量约为pack4策略的1.5倍，或者数据在L1D中热时的吞吐量的3/4，仅计算内部循环。这只是内积部分，没有考虑不同的pack/unpack成本。
  • 难以转置(但使用pmovmskb从每个字节中提取1位并使它们连续可能不是很糟糕)。
每个字节6个元素，0xxx0xxx（在HSW/SKL上可能没有优势，但考虑这一点很有趣）：
  • 比单独字节高6倍的密度。
  • 通过移位/OR，以交错的方式从0/1字节中相对容易地创建，与每个字节4位的格式相同。

  • 针对AVX2的vpshufb进行了高效的popcount优化。在2xvpshufb之前不需要进行掩码处理，只需进行1次右移操作。（vpshufb会将高位设置为0的字节清零，否则就将低半字作为索引。这就是它需要掩码的原因。）将该格式右移4位（vpsrld ymm0,4）仍会在每个字节的高位保留一个零。加载+AND->累积popcount每个向量需要7个融合域uop(vmovdqa/vpand ymm,[mem]/vpsrld ymm,4/2xvpshufb/2xvpaddb)，其中只有6个需要ALU端口。因此，HSW/SKL理论吞吐量为每2个时钟一个向量（192个元素），或每个时钟96个元素。这需要平均每个时钟周期传输一个256位向量，因此正好达到L2带宽瓶颈。

    理论上与完全打包相同，但实际上可能会稍微快一些或慢一些，这取决于哪个调度更好（比如较少的AND/ADD uop从随机移位器窃取端口5）。完全打包可能更有可能接近理论速度，因为它的更多uop可以在多个端口上运行。无序调度缺陷不太可能出现。

  • pmovmskb转置技巧不干净利落。
  • 如果我们只需要popcount（A[]）而不是popcount（A[]＆B[]），或针对不同的微体系结构，其中ALU vs.负载吞吐量不同，则可能非常有用。
另一种变化是，每个字节的7个元素可以通过单个AVX512VBMI（Cannonlake？）进行popcount计算{{link1：vpermi2b（_mm512_permutex2var_epi8）}}，其中每个索引字节从另外两个寄存器的连接中选择128个字节之一。这样宽的洗牌可能会很慢，但希望其吞吐量比AVX512 vpshufb分离的nibble要好。
使用AVX512VBMI（但不使用AVX512VPOPCNTDQ）来计算打包的8，您可以使用vpermi2b计算低7位，然后移位+掩码顶部位并将其相加。（单个位的popcount =该位）。
元素更容易有效地洗牌（因为有类似于的字节洗牌），因此如果您必须即时转置，则可能值得考虑。或者只在转置时即时压缩到位？

32位整数也是一种选择，但不是一个好的选择。每个向量中较少的元素意味着在转置中有较少的洗牌指令，但不是4的倍数。转置中的洗牌数量可能随着每个向量中元素的对数而扩展。

这对于缓存占用/内存带宽也很重要。大小差异的因素为8可以意味着完成整行或整列仅占据L1的一部分，而不是溢出L1。因此，它可以使缓存块化更容易/不太重要。

每个矩阵占用 23.84MiB，使用压缩位元素。这比 L2 缓存（Haswell 上为 256kiB，Skylake-AVX512 上为 1MiB {{link1}}）大得多，但适合许多核 Xeon CPU 的 L3。但 L3 被所有核心（包括云环境中的其他虚拟机）竞争性地共享，并且比 L2 慢得多。（像你在 HPC/云系统中运行的许多核心 Xeons，由于与 L3 缓存的更高延迟以及并发性的不增加，其每个核心的内存带宽比四核桌面系统低，详见{{link2}}中的“延迟受限平台”部分回答。在 Xeon 上驱动相同数量的内存带宽需要更多的核心，尽管总吞吐量更高。但是，如果每个核心在其私有 L2 中工作，您将获得很多好处。）

将AND结果相加：您已经安排了循环，因此需要将布尔值的单个运行减少为非零计数。这是一件好事。

使用8位整数0/1元素，您可以执行多达255个vpaddb，然后一个元素可能会溢出。它具有良好的吞吐量：在Haswell上每个时钟2次，在Skylake上每个时钟3次。使用多个累加器，这覆盖了许多AND结果向量。对于一个全零向量，使用vpsadbw以将向量中的字节水平添加到64位整数中。然后使用vpaddq组合您的累加器，然后进行水平求和。

使用紧缩位时，您只需要对 AND 结果的向量进行 popcount。 使用 AVX2 并且您的数据已经是向量化的，则绝对要使用基于 VPSHUFB 的位切片 popcount。 （例如，请参见 http://wm.ite.pl/articles/sse-popcount.html。如果您必须手动矢量化它，则应使用内置函数而不是汇编语言来编写它。）

您可以考虑将数据以每字节4位的方式打包在低四位中。这意味着一个vpshufb可以计算AND结果中每个字节的位数，而不需要任何移位/掩码操作。在内部循环中，您将有2个加载操作、vpand、vpshufb和vpaddb。通过适当展开，这应该可以跟上每时钟周期2x 32B的L1D负载带宽，并饱和Haswell或Skylake上的所有三个向量执行端口。每128或255个向量左右中断一次，使用vpsadbw/vpaddq累加累加器的字节。(但是通过缓存块，您可能想要经常中断并进行不同的工作)。因此，如果您可以安排它读取L1D缓存中热数据，则最内层循环应以每字节4个元素* 32B每向量=每时钟周期128个元素的速度运行。从Haswell/Skylake的L2缓存预计会获得大约一半的带宽，或者从L3缓存中获得更差的性能。

使用值为0或1的uint8_t元素，您可以尝试使用一些整数乘加指令。它们的设计有点奇怪，旨在用于不同于FP FMA的用例。它们添加成对的乘积结果，产生更宽的元素。VPMADDUBSW将8位扩展为16位元素，并且在值为0或1时效果很好。由于每个元素只能在范围0..2内，因此仍然可以使用vpsadbw进行水平求和。但是如果您要使用vpsadbw，则与vpand相比，这并没有任何好处。它只在您想要使用vpaddw在向量累加器中使用16位元素而不是退出循环以避免字节溢出时才有用。vpmaddubsw在这里似乎没有用，因为vpsadbw是按字节水平相加的更好方法。

将0/1整数转换为位图可以使用SSE/AVX高效完成：对于32位整数元素，使用vpslld ymm0, 31将相关位左移到每个元素的顶部，然后使用vmovmskps eax, ymm0获取每个32位元素的高字节的8位掩码。对于8位整数元素，使用vpslld ymm0, 7/vpmovmskb eax, ymm0做同样的事情，但是对于每个字节，产生一个32位整数位图结果。（只有每个字节的符号位很重要，所以没有仅具有8位粒度的移位指令也没关系。您不需要对进入下一个元素的位做任何处理。）
这不是立即与向量一起使用的非常好的方法，因为你最终会得到整数寄存器中的结果。这不是一个很好的格式来生成和实时使用，但它是最紧凑的，所以如果您可以长期保留矩阵在这种格式中，这可能是有意义的。（如果在加载它们时受到内存带宽的限制。）
将32位整数转换为8位: 一种方法是使用2x vpackssdw + vpacksswb。因为它们在128b通道内操作，所以您的元素将被重新排序。但只要每行/列的排序相同，这就没有问题。仅当您想要获取不以32个元素的倍数开头的行/列的块时，才会出现问题。另一个选择是左移(8、16和24)，然后将向量OR在一起。实际上，通过使用偏移1、2或3个字节的非对齐加载，您可以免费进行移位。 
static inline
__m256i load_interleave4x32(const int32_t *input) {
  const char *p = (const char*)input;
  __m256i t0 = _mm256_load_si256((const __m256i*)(p));
  __m256i t1 = _mm256_load_si256((const __m256i*)(p+32*1-1));  // the 1/0 bits will be in the 2nd byte of each 32-bit element
  __m256i t2 = _mm256_load_si256((const __m256i*)(p+32*2-2));
  __m256i t3 = _mm256_load_si256((const __m256i*)(p+32*3-3));
  return t0 | t1 | t2 | t3;
  // or write this out with _mm256_or_si256, if you don't have overloaded operators like GNU C does.
  // this should compile to 1 load and 3 vpor ymm0, [rdi+31] ... instructions.
}

将每个字节压缩为半包 4 位: 我们可以使用与上面相同的思路。从 load_interleave4x32 中获取 4 个向量（或从一个 uint8_t 数组中获取，如果您始于 8 位元素）。将它们左移 0、1、2 和 3 位，并将它们全部 OR 在一起。这种交错的位顺序在我们只需要对行/列进行 AND 并弹出整个结果时是没问题的，因为顺序不重要。这种位顺序在解压回有序字节时相当高效，例如，AND with set1_epi8(1) 将会给你一个字节向量。
如果您将整个矩阵存储在此格式中，则可以将其用作转置的一部分，或者您可以使用此格式来存储临时副本以进行高速缓存块转置。由于 matmul 多次访问每行/列，因此进行额外的工作以首次创建紧凑格式可能是值得的，这样可以在后续传递中使每个向量执行 4 倍的工作。

使用AVX512BW (Skylake-AVX512)

我们希望使用向量而不是标量整数来执行AND和popcnt，因为向量的宽度是AVX2的两倍，所以它们比标量popcnt更快。尽管Skylake-AVX512在运行512b指令时关闭了端口1上的向量ALU（但未关闭标量ALU）。

@Harold 指出了一个有趣的等式，使我们可以进行2/3的向量popcounts，但需要进行额外的整数操作。

   popcnt(a) + popcnt(b) + popcnt(c)
 = popcnt(a ^ b ^ c) + 2 * popcnt((a ^ b) & c | (a & b))
a ^ b ^ c 和 (a ^ b) & c | (a & b) 每个都可以使用一个vpternlogd（因为每个都有3个布尔输入）。 如果我们使用单独的预先移位的vpshufb LUT向量，则2*是免费的。 另请参见此实现，它使用30个vpternlogd+ 1个矢量popcnt来处理16个512b向量，最后进行一些清理（仅在循环内执行16 * popcnt计数；其他所有内容都被链接）。
这对于计算完全填充的每字节8位元素非常值得，并使该格式比为popcounting优化而没有太多移位/掩码的低密度格式更具吸引力。

vpternlogd也可以作为位混合指令用于转置，如果字节粒度的VPBLENDMB zmm{k1}, zmm, zmm不够细粒度。

对于某些CPU上的AVX2而言，这可能是值得的，也许可以避免每4或5个向量popcounts中的1个，而不是3个之一？或者，如果它只增加了总执行端口压力，并且没有任何特定的瓶颈，那么它可能根本没有帮助。对于标量popcnt指令（可能在没有AVX2的CPU上），它将非常有用，因为这些指令在Intel CPU上会出现单个端口的瓶颈。

我们可以比AVX2更高效地将uint8_t布尔元素转换为非交错位图（甚至不需要移位），并且可以更高效地执行相反操作。测试到掩码或与set1_epi8（1）向量进行比较都可以完成任务，从64个字节的输入产生64位掩码。或者从32位整数开始，每次产生16位掩码。您可以使用kunpck指令有效地连接这些位。 _mm512_test_epi8_mask (vptestmb) 很有趣：将两个向量进行AND运算，并生成一个掩码寄存器结果，其中包含了真/假的字节元素。但这并不是我们想要的：如果我们要打包位，我们希望在输入矩阵上作为预处理步骤进行，而不是在执行内积时进行。

将位图转换为 0 / -1 的向量是快速的：__m512i _mm512_movm_epi8 (__mmask64 k) (vpmovm2b) 可以在一条指令中完成。你可以使用减法 -1 替代加法 1，但是在将一个字节内的多个位进行 OR 运算之前，必须先对它进行掩码处理。

没有AVX512BW或AVX512DQ（骑士着陆Xeon Phi），你没有512b的vpshufb，因此无法高效地进行向量popcnt。有一个AVX512 popcnt扩展可以直接进行向量popcnt，但是尚未宣布任何拥有它的硬件。（在KNL上，AVX2 vpshufb ymm非常慢：每12个周期一次，而psadbw ymm每9个周期一次，因此即使使用256b向量也不太吸引人）。您可能会使用基于32位整数元素的位操作popcnt，因为那只是AND/shift/ADD。32位元素将比64位元素更快地进行popcnt，但仍足够大，以便于合理的问题大小不会溢出（因此您可以将向量的水平总和推迟到循环外部）

在选择存储格式时，将多个比特打包到一个字节中可能不是 KNL 的好选择，但单字节整数元素是好的。vpandd zmm和vpaddd zmm都非常快且属于 AVX512F，我们可以使用它们，因为我们不想让单字节溢出。(当我们实际有 8 位元素不会相互进位时，使用打包的 32 位加法是一种SWAR技术)。我认为相对于 Skylake-AVX512，KNL 具有良好的内存带宽和较差的指令吞吐量。

位翻转：

BMI2 _pdep_u64 在这里可能很有用。它是一个标量指令/内嵌函数。如果它使位翻转比解包为字节更加高效，那么你可能想要在使用 AND + count 进行向量加载之前存储一块翻转结果。 （在标量存储后立即重新加载向量将导致存储转发停顿。）

另一个有用的选项是vpmovmskb可以从32字节向量中切割出32位，每个字节一个。这为转置提供了一个构建块，也许可以与字节洗牌结合使用以获得正确顺序的字节。更多内容，请参见这篇博客文章，还有如何转置二进制矩阵？。

在matmul中使用

你的选择取决于输入数据的格式以及你将重复使用相同矩阵的频率。如果一个矩阵将被多次使用，则事先将其压缩为每字节4或8位是有意义的。(或者在第一次使用时动态压缩)。保留其转置副本也可能是有意义的，特别是如果它总是需要转置的一侧进行乘法。(如果你有时需要一种方式，有时需要另一种方式，那么动态重做可能更适合L3缓存占用。但这些足够大，你可能不会得到很多L3命中，所以只保留一个转置副本可能是好的。)

或者在从输入格式转换时甚至编写一个转置和非转置版本。

你肯定想要缓存块乘法，这样相同的数据在L1中热重复使用。我脑海中没有任何有用的话可说。与缓存块正常FP matmul相同的原则适用，因此请阅读相关内容。

对你的C++实现的评论：

使用&位集(bitset)处理整个列会将值放回内存，然后在结果上再次循环遍历.count()。我怀疑编译器是否会将其优化为一次使用基于VPSHUFB的比特切片popcnt的循环，用于每个VPAND结果的向量。但是这将更好。(例如参见http://wm.ite.pl/articles/sse-popcount.html。如果必须手动向量化，则应使用intrinsic而不是asm进行编写。)

对于你的矩阵大小，至少那个内部循环可能会命中L1D缓存，但从两次循环中的额外load/store指令产生的开销更大，并且它还干扰了有价值数据的预取。

动态大小位图的高效popcnt编译器（无需手动矢量化）并不容易。唯一不那么糟糕的是clang++ -stdlib=libc++与vector<bool>，它将std::count(v.begin(), v.end(), true);编译为矢量化的vpshufb+vpsadbw+vpaddq循环，这相当不错。如果仅在展开的循环内部使用vpaddb，并且每次迭代仅使用vpsadbw + vpaddq，则速度会更快，但对于自动矢量化的代码来说已经很好了。
g++的vector<bool>也是位图，但std::count(v.begin(), v.end(), true);非常糟糕：它使用完全天真的循环，逐位测试。而且甚至没有高效地执行。使用默认的libstdc++而不是新的libc++的clang++也是如此。

boost::dynamic_bitset有一个.count()成员函数，但它没有利用popcnt指令或AVX2。它执行逐字节的LUT查找。这比没有libc++的std::count(vector<bool>)要好得多，但对于HPC来说还远远不够好。

这是测试代码在Godbolt编译器浏览器上，使用了gcc和clang的汇编输出。所有的代码都使用了-march=haswell。
但不幸的是，似乎没有一种有效的方法可以对两个std::vector<bool>进行按位与操作。这个答案展示了如何获取g++的libstdc++vector<bool>的底层实现，但该代码不支持自动向量化。对于libc++，做同样的事情并调整它以使其自动向量化可能可以让您获得手动向量化可能具有的良好性能的很大一部分（除了转置），但您可能需要将整个矩阵保存在一个vector<bool>中，因为向量的向量是一个糟糕的额外间接级别。如果问题的转置部分也是性能关键，则使用标准容器来访问高效的popcount不能解决整个问题。
对于 std::bitset<1024*1024>.count()，Clang 使用相同高效的 AVX2 popcount 无论是否使用 libc++。g++ 使用标量循环，使用 64 位的 popcnt 指令，这在 Haswell 和 Skylake 上对于小位集比好的 AVX2 popcnt 稍快，但对于大位集稍慢（根据 this）。另请参见：关于 vector<bool> — Howard Hinnant，了解有关 C++ 标准库的评论以及为什么位数组是一种有用的数据结构，但 vector<bool> 是一个糟糕的名称。此外，还有一些基准测试，用于优化计数/查找等操作的位向量与每字节 1 个 bool 的 bool[] 数组相比，与朴素的 vector<bool>（如 gcc 和 clang 中没有 libc++ 的那样）。
1
@NULL: “一次完成”这个短语有点无稽之谈。我忘记了缓存阻塞需要你多次回到同一行/列。但是，如果你真的想节省每一毫秒，我认为你将不得不手动进行向量化。我更新了我的答案，并提供了更多关于您希望编译器生成什么样的汇编代码（可能使用内部函数）的详细信息。 - Peter Cordes
1
@NULL: 你可以尝试使用clang和vector<bool>来获得相当不错的结果。它可以执行popcnt，但我没有查看两个向量之间的 & 运算符。但是，当然，“相当不错”只适用于时间循环内部的部分，不包括转置或打包成位图的部分。这就是我不指望你能在没有内部函数的情况下让编译器很好地完成的部分（尤其是因为你可以通过位交错使其更有效率，就像我在我的回答中提到的那样）。 - Peter Cordes
1
@NULL：如果您有兴趣雇用某人手动矢量化此代码（并针对Haswell / Skylake进行优化），以适应您特定的用例，我很乐意为您提供服务。无论如何，如果您还有更多问题，我很乐意回答。 - Peter Cordes
1
我不知道这个是否有多大/有没有帮助，但您可以通过 popcnt(x) + popcnt(b) + popcnt(c) = popcnt(a ^ b ^ c) + 2 * popcnt((a ^ b) & c | (a & b)) 将每3个popcnt变为2个。 这将交换一个popcnt来执行6个简单的指令，如果使用 vpternlogd 将更好。这可以扩展，但是评论篇幅有限，无法展开。 - harold
1
我发现了一个非常好的实现，甚至更进一步：https://github.com/WojciechMula/sse-popcount/blob/master/popcnt-avx512-harley-seal.cpp - harold
显示剩余9条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接