矩阵向量乘法优化 - 缓存大小

Question

这个问题是关于 C++ 优化技术的。我有一个大尺寸的矩阵向量乘法，并希望减少运行时间。我知道有专门的线性代数库，但我实际上想了解一些关于底层处理器特性的知识。到目前为止，我正在使用 \O2 (Microsoft) 进行编译，并让编译器确认乘法的内部循环是矢量化的。

示例代码为：

#include <stdio.h>
#include <ctime>
#include <iostream>

#define VEC_LENGTH 64
#define ITERATIONS 4000000

void gen_vector_matrix_multiplication(double *vec_result, double *vec_a, double *matrix_B, unsigned int cols_B, unsigned int rows_B)
{
    // initialise result vector
    for (unsigned int i = 0; i < rows_B; i++)
    {
        vec_result[i] = 0;
    }
    // perform multiplication
    for (unsigned int j = 0; j < cols_B; j++)
    {
        const double entry = vec_a[j];
        const int col = j*rows_B;

        for (unsigned int i = 0; i < rows_B; i++)
        {
            vec_result[i] += entry * matrix_B[i + col];
        }
    }
}

int main()
{
    double *vec_a = new double[VEC_LENGTH];
    double *vec_result = new double[VEC_LENGTH];
    double *matrix_B = new double[VEC_LENGTH*VEC_LENGTH];

    // start clock
    clock_t begin = clock();

    // this outer loop is just for test purposes so that the timing becomes meaningful
    for (unsigned int i = 0; i < ITERATIONS; i++)
    {
        gen_vector_matrix_multiplication(vec_result, vec_a, matrix_B, VEC_LENGTH, VEC_LENGTH);
    }

    // stop clock
    double elapsed_time = static_cast<double>(clock() - begin) / CLOCKS_PER_SEC;
    std::cout << elapsed_time/(VEC_LENGTH*VEC_LENGTH) << std::endl;

    delete[] vec_a;
    delete[] vec_result;
    delete[] matrix_B;

    return 1;
}

多次执行乘法以获得对运行时间的可靠估计。我已经测量了许多不同向量长度的运行时间（在这个例子中只有一个元素N，这是向量的长度，同时定义了矩阵的大小NxN）和将测量的运行时间标准化为元素的数量。

您可以看到，对于足够小的N，每个操作的运行时间是恒定的。但是，在N=512 之上，运行时会跳起来。蓝色和红色数据点之间的区别在于处理器的负载。如果示例程序几乎是单独运行，则运行时间由蓝点给出，而当其他核心忙时，时间由红点表示。

我现在有几个与此相关的问题。

• 我是否正确假设N=512 和N=1024 之间的跳转与我的处理器（Ivy Bridge i5-3570）的 L3 缓存大小（应该为 6MB）有关？ 512*512*8byte 大约等于 2MB，1024*1024*8byte 大约是 8MB。所以矩阵不再适合缓存，因此从 RAM 中获取数据是执行时间较长的原因？
• 运行时间稳定增长超过此阈值的原因是什么？
• 对于繁忙和空闲处理器的曲线在阈值以上如此不同的原因有什么想法吗？
• 在优化此乘法例程以使用 N>1024 操作时，接下来的逻辑步骤是什么？

我很想听听你的想法。谢谢！

Answer 1

优化此类代码的一个重要方面是处理别名和矢量化，您的帖子表明您已经处理了后者。编译器经常需要一些帮助。在 GCC 5.3.0 上，使用下面的循环大大减少了运行时间。 __restrict__ 限定符告诉编译器不可能有别名，#pragma GCC ivdep 告诉 GCC 编译器可以向量化代码。此外，编译器标志也非常重要。我使用g++ -O3 -march=native -mtune=native matrix_example.cxx编译了代码。

void gen_vector_matrix_multiplication(double* const __restrict__ vec_result,
                                      const double* const __restrict__ vec_a,
                                      const double* const __restrict__ matrix_B,
                                      const int cols_B,
                                      const int rows_B)
{
    // initialise result vector
#pragma GCC ivdep
    for (int i = 0; i < rows_B; i++)
        vec_result[i] = 0;

    // perform multiplication
    for (int j = 0; j < cols_B; j++)
    {
        const double entry = vec_a[j];
        const int col = j*rows_B;

#pragma GCC ivdep
        for (int i = 0; i < rows_B; i++)
        {
            vec_result[i] += entry * matrix_B[i + col];
        }
    }
}

Answer 2

为了标准化，我选择了

elapsed_time/(VEC_LENGTH*VEC_LENGTH*ITERATIONS)

它以 6 纳秒开始，以 7 纳秒结束，从 N=64 到 N=8192

ITERATIONS=20

所有情况下唯一缓存的东西是“vec_a”，因此只有矩阵元素从内存中读取大矩阵。

内存带宽约为 20 GB/s，这意味着每秒超过 2 G 双倍。 核心频率为 3.7 GHz，因此这将是最大 3.7 G 倍频。

核心每秒可以发出 3.7 G 双倍，但内存每秒可提供 2 G 的含义。

当然这仅适用于 64 位 fp 操作。还有

 i + col

必须在乘法之前完成，所以这是一个串行执行。 3.7 GHz 的 2 条指令实际上意味着接近 1.8 G / 秒。接近2。 即使缓存完成了它的工作，cpu 核心也缺乏该串行代码的计算能力。

将循环展开 4 时发生了同样的事情。这将时间减少了一半！。现在每次操作是 3.4 纳秒，但对于所有 N 值，因为在 CPU 必须执行 1 个内存带宽单位之后仍然有 2 个指令（1 个整数和 1 个浮点）。

编辑：使用所有内核将超过内存带宽，并使 L3 缓存的效果更加明显。