将向量的向量转换为具有相反存储顺序的单个连续向量的更快方法

Question

我std::vector<std::vector<double>>正在尝试尽快将其转换为单个连续向量。我的向量的形状大致为4000 x 50。

问题是，有时我的输出向量需要以列为主的连续顺序（仅将2d输入向量的内部向量连接起来），有时我的输出向量需要以行为主的连续顺序，实际上需要转置。

我发现，朴素的for循环对于转换为以列为主的向量非常快：

auto to_dense_column_major_naive(std::vector<std::vector<double>> const & vec)
    -> std::vector<double>
{
    auto n_col = vec.size();
    auto n_row = vec[0].size();
    std::vector<double> out_vec(n_col * n_row);
    for (size_t i = 0; i < n_col; ++i)
        for (size_t j = 0; j < n_row; ++j)
            out_vec[i * n_row + j] = vec[i][j];
    return out_vec;
}

但是显然，由于所有的高速缓存未命中，因此类似的方法对于逐行转换非常慢。因此，对于逐行转换，我认为提高缓存局部性的阻止策略可能是我最好的选择：

auto to_dense_row_major_blocking(std::vector<std::vector<double>> const & vec)
    -> std::vector<double>
{
    auto n_col = vec.size();
    auto n_row = vec[0].size();
    std::vector<double> out_vec(n_col * n_row);
    size_t block_side = 8;

    for (size_t l = 0; l < n_col; l += block_side) {
        for (size_t k = 0; k < n_row; k += block_side) {
            for (size_t j = l; j < l + block_side && j < n_col; ++j) {
                auto const &column = vec[j];
                for (size_t i = k; i < k + block_side && i < n_row; ++i)
                    out_vec[i * n_col + j] = column[i];
            }
        }
    }
    return out_vec;
}

这比行主要转换的朴素循环快得多，但在输入大小上仍比朴素列主要循环慢了一个数量级。

我的问题是，是否有更快的方法将双精度向量的（列主）向量转换为单个连续的行主向量？我正在努力思考这段代码的速度限制应该是多少，从而质疑我是否缺少明显的东西。我的假设是，阻断会给我多少，然后更大的加速似乎实际上给。

---

该图表是使用QuickBench生成的（并通过我的计算机上的GBench进行了某种程度的验证），其代码如下：（Clang 7，C ++ 20，-O3）

auto to_dense_column_major_naive(std::vector<std::vector<double>> const & vec)
    -> std::vector<double>
{
    auto n_col = vec.size();
    auto n_row = vec[0].size();
    std::vector<double> out_vec(n_col * n_row);
    for (size_t i = 0; i < n_col; ++i)
        for (size_t j = 0; j < n_row; ++j)
            out_vec[i * n_row + j] = vec[i][j];
    return out_vec;
}

auto to_dense_row_major_naive(std::vector<std::vector<double>> const & vec)
    -> std::vector<double>
{
    auto n_col = vec.size();
    auto n_row = vec[0].size();
    std::vector<double> out_vec(n_col * n_row);
    for (size_t i = 0; i < n_col; ++i)
        for (size_t j = 0; j < n_row; ++j)
            out_vec[j * n_col + i] = vec[i][j];
    return out_vec;
}

auto to_dense_row_major_blocking(std::vector<std::vector<double>> const & vec)
    -> std::vector<double>
{
    auto n_col = vec.size();
    auto n_row = vec[0].size();
    std::vector<double> out_vec(n_col * n_row);
    size_t block_side = 8;

    for (size_t l = 0; l < n_col; l += block_side) {
        for (size_t k = 0; k < n_row; k += block_side) {
            for (size_t j = l; j < l + block_side && j < n_col; ++j) {
                auto const &column = vec[j];
                for (size_t i = k; i < k + block_side && i < n_row; ++i)
                    out_vec[i * n_col + j] = column[i];
            }
        }
    }
    return out_vec;
}

auto to_dense_column_major_blocking(std::vector<std::vector<double>> const & vec)
    -> std::vector<double>
{
    auto n_col = vec.size();
    auto n_row = vec[0].size();
    std::vector<double> out_vec(n_col * n_row);
    size_t block_side = 8;

    for (size_t l = 0; l < n_col; l += block_side) {
        for (size_t k = 0; k < n_row; k += block_side) {
            for (size_t j = l; j < l + block_side && j < n_col; ++j) {
                auto const &column = vec[j];
                for (size_t i = k; i < k + block_side && i < n_row; ++i)
                    out_vec[j * n_row + i] = column[i];
            }
        }
    }
    return out_vec;
}

auto make_vecvec() -> std::vector<std::vector<double>>
{
  std::vector<std::vector<double>> vecvec(50, std::vector<double>(4000));
  std::mt19937 mersenne {2019};
  std::uniform_real_distribution<double> dist(-1000, 1000);
  for (auto &vec: vecvec)
   for (auto &val: vec)
       val = dist(mersenne);
  return vecvec;
}

static void NaiveColumnMajor(benchmark::State& state) {
  // Code before the loop is not measured

  auto vecvec = make_vecvec();
  for (auto _ : state) {
    benchmark::DoNotOptimize(to_dense_column_major_naive(vecvec));
  }
}
BENCHMARK(NaiveColumnMajor);

static void NaiveRowMajor(benchmark::State& state) {
  // Code before the loop is not measured

  auto vecvec = make_vecvec();
  for (auto _ : state) {
    benchmark::DoNotOptimize(to_dense_row_major_naive(vecvec));
  }
}
BENCHMARK(NaiveRowMajor);

static void BlockingRowMajor(benchmark::State& state) {
  // Code before the loop is not measured

  auto vecvec = make_vecvec();
  for (auto _ : state) {
    benchmark::DoNotOptimize(to_dense_row_major_blocking(vecvec));
  }
}
BENCHMARK(BlockingRowMajor);

static void BlockingColumnMajor(benchmark::State& state) {
  // Code before the loop is not measured

  auto vecvec = make_vecvec();
  for (auto _ : state) {
    benchmark::DoNotOptimize(to_dense_column_major_blocking(vecvec));
  }
}
BENCHMARK(BlockingColumnMajor);

Answer 1

首先，只要某物被认定为“显而易见”，我就会畏缩。这个词经常被用来掩饰自己推论的一个缺点。

但是显然，由于所有的高速缓存未命中，因此类似的方法对于逐行转换非常慢。

我不确定哪个应该是显而易见的：将显示按行转换，或者由于高速缓存未命中而导致速度变慢。无论哪种情况，我都觉得不太明显。毕竟，这里有两个缓存注意事项，不是吗？一种阅读，一种写作？让我们从阅读的角度来看代码：

row_major_naive

for (size_t i = 0; i < n_col; ++i)
    for (size_t j = 0; j < n_row; ++j)
        out_vec[j * n_col + i] = vec[i][j];

连续读取vec是连续内存的读取：vec[i][0]后跟vec[i][1]，等等。非常适合缓存。所以...缓存未命中？慢？:)也许不是那么明显。

尽管如此，还是有一些需要收集的东西。仅通过“明显”主张才是错误的主张。存在非本地性问题，但它们发生在写作端。（成功写入被50个double值的空格所抵消。）经验测试证实了这种慢速性。因此，也许一个解决方案是翻转所谓的“显而易见”的东西？

行专业翻转

for (size_t j = 0; j < n_row; ++j)
    for (size_t i = 0; i < n_col; ++i)
        out_vec[j * n_col + i] = vec[i][j];

我在这里所做的只是扭转循环。从字面上交换那两行代码的顺序，然后调整缩进。现在，连续的读取可能遍及整个地方，因为它们是从不同的向量读取的。但是，连续写入现在是对连续的内存块的。从某种意义上说，我们处于与以前相同的情况。但是就像以前一样，应该在假设“快速”或“慢速”之前测量性能。

NaiveColumnMajor：3.4秒
NaiveRowMajor：7.7秒
FlippedRowMajor：4.2秒
BlockingRowMajor：4.4秒
BlockingColumnMajor：3.9秒

仍然比朴素的专栏主要转换慢。但是，这种方法不仅比朴素的行专业更快，而且比阻塞的行专业更快。至少在我的计算机上（使用gcc -O3并且显然：P反复执行数千次）。里程可能会有所不同。我不知道花哨的剖析工具会说些什么。问题是有时越简单越好。

对于有趣的事情，我进行了尺寸互换的测试（从4000个元素的50个向量更改为50个元素的4000个向量）。所有方法都受到这种伤害，但是“ NaiveRowMajor”受到的打击最大。值得一提的是，“翻转行专业”落后于阻止版本。因此，正如人们可能期望的那样，完成工作的最佳工具取决于工作的确切程度。

NaiveColumnMajor：3.7秒
NaiveRowMajor：16秒
FlippedRowMajor：5.6秒
BlockingRowMajor：4.9秒
BlockingColumnMajor：4.5秒

（顺便说一句，我还尝试了阻塞版本的翻转技巧。更改很小-约为0.2，与翻转朴素版本相反。也就是说，对于问题的答案，“翻转阻止”比“阻止”慢4000个向量中的50个，但对于我的4000个向量中的5000个向量，速度更快。微调可能会改善结果。）

更新：我对阻塞版本的翻转技巧做了更多测试。这个版本有四个循环，因此“翻转”不像只有两个循环时那样简单。似乎交换两个外部循环的顺序对性能不利，而交换两个内部循环的循环则有利。（最初，我同时完成了这两项工作，并且得到了混合的结果。）当我只交换内部循环时，我测量了3.8秒（在4000-of-50情况下为4.1秒），这使其成为测试中最佳的行优先选项。

行主要混合

for (size_t l = 0; l < n_col; l += block_side)
    for (size_t i = 0; i < n_row; ++i)
        for (size_t j = l; j < l + block_side && j < n_col; ++j)
            out_vec[i * n_col + j] = vec[j][i];

（交换内部循环后，我合并了中间循环。）

至于背后的理论，我想这相当于试图一次写入一个缓存块。写入块后，请尝试重新使用向量（vec[j]），然后再将其从缓存中弹出。耗尽所有源向量后，请移至新的一组源向量，再次一次写入完整的块。