使用按位 AND 和 popcount 而不是实际的 int 或浮点数乘法进行大 (0,1） 矩阵乘法

Large (0,1) matrix multiplication using bitwise AND and popcount instead of actual int or float multiplies?

本文关键字：浮点数 int AND popcount      更新时间：2023-10-16

对于乘以大型二进制矩阵(10Kx20K)，我通常要做的是将矩阵转换为浮点矩阵并执行浮点矩阵乘法，因为整数矩阵乘法非常慢(看看这里)。

不过，这一次，我需要执行超过数十万次这样的乘法，甚至平均而言，性能提高一毫秒对我来说很重要。

因此，我想要一个int或float矩阵，因为产品可能具有不是 0 或 1 的元素。 输入矩阵元素均为 0 或 1，因此它们可以存储为单个位。

在行向量和列向量之间的内积中(以产生输出矩阵的一个元素)，乘法简化为按位 AND。 加法仍然是加法，但我们可以使用总体计数函数添加位，而不是单独循环它们。

其他一些布尔/二进制矩阵函数或位而不是计数它们，产生位矩阵结果，但这不是我需要的。

下面是一个示例代码，表明将问题形成为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 的解决方案更快的任何解决方案，因为我由于多种原因无法使用它。

我不确定你可以让编译器为你做多少(如果有的话)，而无需使用内部函数或C++向量类包装器手动矢量化(如 Agner Fog 的 VCL，如果你的项目许可证与 GPL 兼容)。 也有一些非GPL包装器。

缓存阻塞矩阵乘法是一门艺术(在这里很重要)，如果你可以使用 Eigen 的现有模板，但使用不同的类，对整数使用按位and，而不是在浮点数上乘法，那就太好了。 我不确定这是否可能。

我做了一些搜索，大多数关于二进制矩阵的文献都是关于产生布尔结果的(包括像这样的 SO 问题)。 向量内积是用 AND 作为乘法完成的，但用 XOR 或 OR 作为加法，而不是 popcount。 也许我缺少一个搜索词来描述恰好是(0,1)矩阵的"正常"矩阵，但乘积不会。

由于每一毫秒都很重要，因此您可能必须手动对此进行矢量化。

这并不是说矢量整数的东西一般很慢，只是矢量整数乘法很慢，特别是与最近的 x86 硬件上的矢量floatFMA 相比(尤其是英特尔，它在 Haswell 及更高版本上的 FP FMA 吞吐量为每时钟 2x 256b 矢量)。

由于您不需要使用布尔元素进行实际乘法，只需要一个 AND(每个时钟吞吐量 3 个向量)，这对您来说不是问题。 每个向量做更多元素所带来的效率增益应该足以弥补每个向量的任何额外成本。

当然，这假设一个整数 matmul 实现使用与等效的 FP matmul 相同的所有缓存阻塞和其他优化，如果您不想(或不知道如何)自己编写它，并且找不到可以为您完成的库，这就是问题所在。

我只是在回答一个问题，即通过最佳实现，它可以有多高效。标题问题的答案是非常肯定的，使用实际乘法是浪费大量时间，尤其是对于 32 位元素。

存储格式选项：

每个字节一个元素 (0/1)：

  • 4 倍的密度float(缓存占用空间/内存带宽/每个矢量的元素数)
  • 使用字节随机播放易于转置
  • 垂直 ADD 很容易，以防万一(例如，用于在外部循环上进行矢量化，以及一次处理多行或多列。 如果数据交错的方式无需额外洗牌即可工作，则可能很好(避免最后的水平和)。

每字节 4 个元素，打包到低半字节中：

  • 4 倍于独立字节密度
  • 使用 AVX2vpshufb弹出计数非常有效。 在 L1D 缓存中输入热时，理论上，您可以加载/AND/累积-a-popcount，每个时钟周期(每个内核)的吞吐量为 128 个 AND-result 元素。 每个时钟 4 个融合域 uops 使每个时钟 4 个的 SKL/HSW 前端问题带宽饱和，并且不会在 3 个矢量 ALU 端口上造成瓶颈，因为其中一个 uops 是纯负载。(另一个负载微保险丝与vpand)。 当L2带宽出现瓶颈(每个周期~一个32B负载)时，以每时钟64个元件运行。见下文。
  • 从整数或打包位图创建速度较慢(但如果将位按交错顺序放入向量中以便有效地打包/解包到按顺序排列字节，而不是强制按顺序打包位)，则还不错)。
  • 难以转置(可能比完全包装更糟糕)

包装位：

  • 8 倍于单独字节的密度，每个 AVX2 向量 256 个元素。
  • 可以从具有非交错存储订单pmovmskb的向量创建。 (不过，对于动态创建不是很有用，因为这会将结果放在整数 reg 中，而不是向量中。 交错位序可能是最好的，尤其是在转置期间解包)。
  • 使用 AVX2 进行弹出计数相当有效：掩码/移位+掩码/2xvpshufb. (9 个融合域 uops(8 个向量 ALU uops)到 AND + 累积 256 个元素(来自 2 行/列向量)的弹出计数，而每字节 4 个策略(来自 4 行/列向量)为 8 个 uops(6 个向量 ALU uops)。 ALU 端口瓶颈将其限制为 L1D 或 L2 的每个时钟 96 个元件。 因此，当pack4 策略在 L2 带宽上遇到瓶颈时，其内部产品吞吐量约为 1.5 倍，或者理论上，仅计算内部循环时，L1D 中热数据吞吐量约为 3/4。 这只是产品内部部分，不考虑不同的包装/拆包成本。
  • 很难转置(但也许不可怕，pmovmskb从每个字节中提取 1 位并使它们连续)。

每字节 6 个元素，0xxx0xxx(在 HSW/SKL 上这个问题可能没有优势，但考虑起来很有趣)：

  • 6 倍于单独字节密度
  • 通过移位/ORing，以交错方式从0/1字节创建相当容易，与每字节4位格式相同。

  • 针对 AVX2vpshufb的高效爆音计数进行了优化。 无需在 2 倍vpshufb之前遮罩，只需 1 次右移。 (如果设置了高位，则vpshufb将字节归零，否则它将使用低半字节作为索引。 这就是为什么它需要屏蔽。 将这种格式右移4(vpsrld ymm0,4)仍然会在每个字节的高位上留下一个零。 加载+AND ->累积弹出计数是每个矢量 7 个融合域 uops (vmovdqa/vpand ymm,[mem]/vpsrld ymm,4/2xvpshufb/2xvpaddb)，其中只有 6 个需要 ALU 端口。 因此，HSW/SKL 吞吐量理论上是每 2 个时钟 1 个向量(192 个元素)，或每个时钟 96 个元素。 这需要每个时钟一个 256b 矢量的平均负载吞吐量，因此它正好解决了 L2 带宽瓶颈。

    从理论上讲，它与完全打包相同，但在实践中，它可能更快或更慢，具体取决于哪个计划更好(例如，更少的 AND/ADD uops 从随机播放中窃取端口 5)。 完全打包可能更有可能接近理论速度，因为它的更多 uops 可以在多个端口上运行。 无序调度缺陷的可能性较小。

  • pmovmskb转置技巧不能干净地工作。
  • 如果我们只需要popcount(A[])而不是popcount(A[] & B[])，可能会很有用。 或者对于不同的微架构，其中 ALU 与负载吞吐量不同。

另一个变体是，每个字节 7 个元素可以用一个AVX512VBMI弹出(Cannonlake？vpermi2b(_mm512_permutex2var_epi8)，其中每个索引字节从另外两个寄存器的串联中选择 128 个字节中的一个。 这么宽的洗牌可能会很慢，但它有望比 AVX512vpshufb单独半字节的东西具有更好的吞吐量。

要用 AVX512VBMI(但没有 AVX512VPOPCNTDQ)计算 packed-8，您可以使用vpermi2b来计算低 7，然后 shift+掩盖顶部位并添加它。 (单个位的弹出计数 = 该位)。

uint8_t元素更容易有效地洗牌(因为有像vpshufb这样的字节洗牌)，所以如果你必须即时转置，可能值得考虑。 还是在转调时只打包成位？

32 位整数也是一种选择，但不是一个好的选择。 每个向量的元素越少意味着转调中的随机播放指令越少，但不是 4 倍。 转置中的随机播放次数可能会像 log2(每个向量的元素)这样缩放。

这对于缓存占用空间/内存带宽来说也是一个大问题。 8 大小差异的因子可能意味着执行整行或整列仅占用 L1 的一部分，而不是溢出 L1。 因此，它可以使缓存阻塞更容易/不那么重要。

每个矩阵 10k * 20k/8 = 23.84MiB，使用打包位元素。 这比L2缓存(Haswell为256kiB，Skylake-AVX512为1MiB)大得多，但适合众核至强CPU上的L3。 但 L3 由所有核心(包括云环境中的其他 VM)竞争共享，并且比 L2 慢得多。 (像您这样的众核至强将在HPC/云系统中运行，其每核心内存带宽低于四核台式机，因为L3缓存的延迟更高，而并发性没有增加(请参阅本答案的"延迟受限平台"部分。 在至强上驱动相同数量的内存带宽需要更多的内核，即使总吞吐量更高。 但是，如果你可以让每个内核主要在其私有 L2 之外工作，你会获得很多。

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

使用 8 位整数 0/1 元素，您最多可以在元素溢出之前执行 255vpaddb。 它具有良好的吞吐量：Haswell上每个时钟2个，Skylake上每个时钟3个。 使用多个累加器，这涵盖了许多 AND 结果的向量。 对全零向量使用vpsadbw将向量中的字节水平相加为 64 位整数。 然后将累加器与vpaddq相结合，然后水平相加。

对于打包位，您只想弹出 AND 结果的向量。 对于 AVX2 并且您的数据已经在矢量中，您肯定希望使用VPSHUFB位切片弹出计数。 (例如，请参阅 http://wm.ite.pl/articles/sse-popcount.html。 如果你必须手动矢量化它，你会想用内联函数而不是 asm 来编写它。

您可以考虑将数据打包为每字节 4 位，以低半字节为单位。这意味着一个vpshufb可以计算 AND 结果的每个字节中的位，而无需任何移位/掩码。 在内部循环中，您将有 2 个负载，vpand、vpshufb、vpaddb。 通过适当的展开，这应该跟上每个时钟 2x 32B 的 L1D 负载带宽，并使所有三个矢量执行端口饱和(在 Haswell 或 Skylake 上)。 每 128 或 255 个向量或其他东西突破一次，以累积累加器的字节数vpsadbw/vpaddq. (但是对于缓存阻塞，您可能希望经常打破并执行不同的工作)。因此，最里面的循环应该以每字节 4 个元素 * 每个向量 32B = 每个时钟周期 128 个元素运行，如果你可以安排它读取 L1D 缓存中的热数据。 预计Haswell/Skylake上的L2缓存的带宽约为一半，或者L3缓存的带宽要差得多。

对于 0 或 1 的uint8_t元素，您可以使用一些整数乘加指令。 它们的设计有点奇怪，用于与FP FMA不同的用例。 它们添加水平乘法结果对，产生更广泛的元素。VPMADDUBSW从 8 位元素扩展到 16 位元素，并且可以很好地处理 0 和 1。 由于每个元素只能在 0..2 范围内，因此您仍然可以使用vpsadbw进行水平和。 但是如果你要vpsadbw，这对你没有任何好处vpand。 仅当您想使用vpaddw在矢量累加器中使用 16 位元素而不是中断循环以避免字节溢出时，它才有用。vpmaddubsw doesn't seem useful here, becausevpsadbw' 是水平添加字节的更好方法。

使用 SSE/AVX 可以有效地将 0/1 整数转换为位图：对于 32 位整数元素，vpslld ymm0, 31将相关位左移到每个元素的顶部，然后vmovmskps eax, ymm0获取每个 32 位元素的高字节的 8 位掩码。 对于 8 位整数元素，vpslld ymm0, 7/vpmovmskb eax, ymm0执行相同的操作，但对于每个字节，生成 32 位整数位图结果。 (只有每个字节的符号位很重要，所以没有只有 8 位粒度的移位指令是可以的。 您无需对进入下一个元素的位执行任何操作。

对于立即使用向量来说，这不是一个很好的方法，因为最终会得到整数寄存器中的结果。 这不是一种可以动态生成和使用的好格式，但它是最紧凑的，因此如果您可以长期保持这种格式的矩阵，它是有意义的。 (如果您在加载它们时会受到内存带宽的限制。

将 32 位整数转换为 8 位：一种方法是使用 2xvpackssdw+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(如果从 8 位元素开始，则从uint8_t数组中获取 4 个向量)。 将它们左移 0、1、2 和 3 位，以及 OR 一起移动。 当我们只需要 AND 一行/列并弹出整个结果时，这种交错位顺序很好，因为顺序无关紧要。 这种位顺序对于解包回按顺序的字节是相当有效的，例如，AND withset1_epi8(1)会得到一个字节向量。

如果以这种格式存储整个矩阵，则可以将其用作转置的一部分，也可以使用此格式存储缓存阻止转置的临时副本。 矩阵多次接触每一行/每一列，因此第一次做一个紧凑的格式可能是值得的，因为这可以让你在后续的传递中为每个向量做 4 倍的工作。

带AVX512BW (Skylake-AVX512)

我们真的想用向量而不是标量整数来做 AND 和 popcnt，因为向量的宽度是 AVX2 的两倍，所以它们比标量popcnt更靠前。 (即使 Skylake-AVX512 在运行 512b 指令时关闭端口 1 上的矢量 ALU(但不是标量)。

@Harold指出了一个有趣的恒等式，它允许我们执行 2/3 的向量弹出计数，但代价是额外的整数运算。

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 个布尔输入)。 如果我们使用单独的预移位vpshufbLUT 矢量，则2*是免费的。 另请参阅此实现，它使用 30xvpternlogd+ 1 向量 popcnt 来处理 16 个 512b 向量，最后进行一些清理(仅在循环内执行16*popcnt计数;其他所有内容都是链接的)。

这对于计算每字节 8 位元素的完全打包来说很可能是值得的，并且与针对弹出计数优化的密度较低的格式相比，该格式对 AVX512 更具吸引力，而无需太多的移位/掩蔽。

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

对于某些 CPU 上的 AVX2 来说，这可能是值得的，也许每 4 或 5 个矢量弹出计数中就有 1 个，而不是 3 个中的 1 个？ 或者，如果它只是增加总执行端口压力，并且任何特定端口都没有瓶颈，则可能根本无济于事。 这对于标量popcnt指令(可能在没有 AVX2 的 CPU 上)很有用，因为这些指令在英特尔 CPU 上的单个端口上存在瓶颈。

我们可以uint8_t布尔元素转换为非交错位图，比 AVX2 稍微更有效(甚至不需要移位)，并且更有效地执行相反的操作。 测试掩码或比较掩码与 set1_epi8(1) 向量都可以完成这项工作，从 64 字节的输入中生成 64 位掩码。 或者从 32 位整数开始，一次生成 16 位掩码。 您可以有效地将这些位与kunpck指令连接起来。

_mm512_test_epi8_mask(vptestmb)很有趣：将两个向量放在一起，并产生真/假字节元素的掩码寄存器结果。 但这并不是我们真正想要的：如果我们要打包我们的位，我们希望将其作为输入矩阵的预处理步骤，而不是在做内积时动态进行。

位图->0/-1的向量很快：__m512i _mm512_movm_epi8 (__mmask64 k)(vpmovm2b)在一个指令中做到这一点。 你可以减去-1而不是加1，但你必须屏蔽它，然后才能在一个字节内或将多个位放在一起。

如果没有AVX512BW或AVX512DQ(Knight's Landing Xeon Phi)，你就没有512bvpshufb所以你不能有效地矢量popcnt。 有一个AVX512 popcnt扩展直接用于矢量popcnt，但还没有宣布任何硬件。 (不过，AVX2vpshufb ymm在 KNL 上非常慢：每 12 个周期一个，psadbw ymm是每 9 个周期 1 个，因此即使使用 256b 向量也没有吸引力)。 您可以使用基于 32 位整数元素的 bithack popcnt，因为它只是 AND/shift/ADD。 与 64 位元素相比，32 位元素对 popcnt 执行的步骤更少，但仍然足够大，不会溢出合理的问题大小(因此您可以将向量的水平和推迟到循环之外)

考虑到存储格式的选择，对于 KNL 来说，每字节打包多个位可能不是一个好主意，但单字节整数元素是好的。vpandd zmm和vpaddd zmm都很快，并且是AVX512F的一部分，我们可以使用它们，因为我们无论如何都不想让我们的单字节溢出。 (当我们实际上有 8 位元素不会相互携带时，使用打包的 32 位添加是一种 SWAR 技术。 我认为，相对于Skylake-AVX512，KNL具有良好的内存带宽和较差的指令吞吐量。

转置位：

BMI2_pdep_u64在这里可能很有用。 这是一个标量指令/内在的。 如果它使位转置比解压缩为字节更有效，您可能希望在用 AND + count 的矢量加载重新加载它之前存储一个转置结果块。 (在标量存储后立即重新加载矢量将导致存储转发停止。

另一个有用的选择是vpmovmskb可以从 32 字节向量中切出 32 位，每字节一个。 这为您提供了一个转置的构建块，可能与字节随机播放结合使用，以使其按正确的顺序获取字节。 有关更多信息，请参阅此博客文章，以及您将如何转置二进制矩阵？。

在矩阵中使用它

您的某些选择取决于输入数据的格式，以及重复使用相同矩阵的频率。 如果矩阵将被多次使用，则提前将其打包为每字节 4 或 8 位是有意义的。 (或第一次使用时即时)。 保留它的转置副本也可能是有意义的，特别是如果它始终是需要转置的乘法的一侧。 (如果您有时需要一种方式，有时需要另一种方式，则动态重做可能更适合 L3 缓存占用空间。 但是这些足够大，你可能不会得到很多L3点击，所以只保留一个转置的副本可能会很好。

或者甚至可以在从输入格式转换时写出转置和非转置版本。

您肯定希望缓存阻止乘法，因此相同的数据在 L1 中热时会多次重用。 我对此没有任何有用的话要说。与缓存阻塞普通 FP 矩阵时相同的原则适用，因此请阅读。

对C++实施的评论：

对整个列使用位集&会将值放回内存中，然后再次循环访问它们以.count()结果。 我怀疑编译器是否会将其优化为一个单通道循环，该循环在每个VPAND结果向量上使用基于VPSHUFB位切片 popcnt，但这会好得多。 (例如，请参阅 http://wm.ite.pl/articles/sse-popcount.html。 如果你必须手动矢量化它，你会想用内联函数而不是 asm 来编写它。

对于矩阵大小，至少内部循环可能会在 L1D 缓存中命中，但循环两次的额外加载/存储指令会占用更多开销，并且还会干扰有价值数据的预取。

让编译器有效地弹出动态大小的位图(无需手动矢量化)并不容易。 唯一不糟糕的是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);非常糟糕：它使用了一个完全朴素的循环，一次测试1位。 它甚至不能有效地做到这一点。 对于使用默认libstdc++而不是新libc++的clang++也是如此。

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

这是 Godbolt 编译器资源管理器上的测试代码，带有 gcc 和 clang asm 输出。 他们都用-march=haswell.

但不幸的是，似乎没有一种有效的方法来按位 AND 两个std::vector<bool>。 这个答案显示了如何获取g++的libstdc++vector<bool>的底层实现，但该代码不会自动矢量化。 对libc++做同样的事情并对其进行调整以使其自动矢量化可能会让您获得手动矢量化可能的性能的很大一部分(转置除外)，但您可能必须将整个矩阵保留在一个vector<bool>中，因为矢量向量是一个糟糕的额外间接级别。 如果问题的转置部分对性能也很关键，那么使用标准容器来访问高效的popcount并不能解决整个问题。

对于std::bitset<1024*1024>.count()，clang 在有或没有libc++的情况下都会产生同样高效的 AVX2 弹出计数。 g++使用64位popcnt指令制作标量循环，根据这个指令，对于小位集，它比一个好的AVX2 popcnt快一些，但对于大位集来说，在Haswell和Skylake上要慢一些。

另请参阅：Onvector<bool>— Howard Hinnant，关于C++标准库的一些评论，以及为什么位数组是一种有用的数据结构，但vector<bool>是一个不好的名字。 此外，一些在位向量上适当优化计数/find_first/等的基准测试与每字节 1boolbool[]数组，与幼稚的vector<bool>(就像你从 gcc 和 clang 获得的没有 libc++ 一样)。

相关文章：