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

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的操作，优化这个乘法例程的逻辑下一步是什么？

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

- Alexander Büse

---

@tobi303，你确定复杂度吗？这是向量-矩阵乘法而不是矩阵-矩阵。只有N * N次操作。 - Alexander Büse

抱歉，可能错过了这个信息 ;) - 463035818_is_not_a_number

这可能会鼓励优化器将vector_a和matrix_b声明为const double*。如果您正在使用C语言，我建议同时使用'restrict';虽然这不是C++，但一些编译器--例如gcc--实现了其形式作为扩展。 - dmuir

@dmuir 感谢您的建议。我刚试了一下，但时间基本相同。 - Alexander Büse

奇怪的是，只有在循环计数器为有符号时，gcc/clang 才会自动向量化。如果循环计数器为无符号，则 i + col 可能会溢出，使得加载不连续。icc13 可以自动向量化原始代码，但是 gcc 和 clang 只执行标量双精度操作，例如 mulsd。这是我进行自动向量化的版本。请注意，它不需要 -ffast-math 就可以进行向量化，因为矢量*矩阵被卡在了循环多次遍历结果向量或从矩阵进行跨步加载的尴尬位置。虽然可能不值得反转。 - Peter Cordes

@PeterCordes 感谢您的评论！是的，有符号计数器很奇怪。我刚试图并行化内部循环，现在微软编译器也要求我更改为有符号计数器和边界。现在它可以工作了，并且对于最大的N，性能稍微有所提高。 - Alexander Büse

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];
        }
    }
}

- Chiel

---

非常感谢您的回答！您是正确的，向量化确实帮了很大的忙（但是是的，我在问题中展示的数据已经使用了向量化循环）。关于别名提示很有用。我找到了一些有趣的东西可以阅读。但是，在我的特定情况下，使用限定符“restrict”并没有改善情况。 - Alexander Büse

@AlexanderBüse，请检查您的编译器中的限定符是什么。Intel使用restrict，gcc和clang使用__restrict__，如果我没记错的话，微软使用__restrict。 - Chiel

是的，你说得对，它是__restrict，但不幸的是，这并没有改善运行时间。 - Alexander Büse

@AlexanderBüse 如果你已经确认它是矢量化的，我想就没有更多可以获得的了。 - Chiel

是的，可能吧。谢谢你的帮助！ - Alexander Büse

Answer 2

对于标准化，我选择了

elapsed_time/(VEC_LENGTH*VEC_LENGTH*ITERATIONS)

这始于 N=64，持续时间为 6 纳秒，结束于 N=8192，持续时间为 7 纳秒。

ITERATIONS=20

对于所有情况，唯一缓存的是"vec_a"，因此对于大矩阵，只有矩阵元素从内存中读取。

内存带宽约为20 GB/s，这意味着每秒超过2 G个双精度浮点数。核心频率为3.7 GHz，因此最多可以进行3.7 G次乘法运算。

核心可以发出3.7 G个双精度浮点数每秒，但内存每秒提供2 G个双精度浮点数，这意味着。

当然，这仅适用于64位浮点操作。还有

 i + col

这必须在乘法之前完成，因此这是一种串行执行。3.7 GHz的2个指令意味着实际上近1.8 G /秒。接近2。 即使缓存发挥了作用，CPU核心也缺乏计算能力来执行此串行代码。

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

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

- huseyin tugrul buyukisik

---

这很有趣。非常感谢！奇怪的是，每次迭代的时间几乎是恒定的。我也尝试将循环展开4次，现在所有N的运行时间也几乎是恒定的，但是随着N的增加... - Alexander Büse

尝试使用所有核心来克服RAM带宽限制并释放L3缓存优势。 - huseyin tugrul buyukisik

是的，我刚刚尝试了内部循环并行化，对于最大的N，我获得了约20％的改进。在N = 8192以下，并行化的开销实际上会增加运行时间。 - Alexander Büse

并行化外层循环会更高效、更友好的缓存，因为核心现在共享所有向量。但这需要添加本地累加器来获取目标向量第i个元素的总和。 - huseyin tugrul buyukisik