多维数组和数组的性能表现

6

我一直认为,只需要通过乘法一次索引的多维数组比通过两个指针解除引用进行索引的数组更快,因为前者具有更好的局部性和空间节省。

不久前,我进行了一项小测试,结果相当令人惊讶。至少我的callgrind分析器报告说,使用数组的相同函数运行速度略快。

我想知道是否应该更改我的矩阵类的定义,以在内部使用数组。这个类在我的仿真引擎中被几乎所有地方使用,我确实希望找到最好的方法来节省几秒钟。

test_matrix的成本为350 200 020test_array_array的成本为325 200 016。代码是由clang++编译的,并且根据分析器,所有成员函数都已经被内联。

#include <iostream>
#include <memory>

template<class T>
class BasicArray : public std::unique_ptr<T[]> {
public:
    BasicArray() = default;
    BasicArray(std::size_t);
};

template<class T>
BasicArray<T>::BasicArray(std::size_t size)
: std::unique_ptr<T[]>(new T[size]) {}

template<class T>
class Matrix : public BasicArray<T> {
public:
    Matrix() = default;
    Matrix(std::size_t, std::size_t);
    T &operator()(std::size_t, std::size_t) const;
    std::size_t get_index(std::size_t, std::size_t) const;
    std::size_t get_size(std::size_t) const;

private:
    std::size_t sizes[2];
};

template<class T>
Matrix<T>::Matrix(std::size_t i, std::size_t j)
: BasicArray<T>(i * j)
, sizes {i, j} {}

template<class T>
T &Matrix<T>::operator()(std::size_t i, std::size_t j) const {
    return (*this)[get_index(i, j)];
}

template<class T>
std::size_t Matrix<T>::get_index(std::size_t i, std::size_t j) const {
    return i * get_size(2) + j;
}

template<class T>
std::size_t Matrix<T>::get_size(std::size_t d) const {
    return sizes[d - 1];
}

template<class T>
class Array : public BasicArray<T> {
public:
    Array() = default;
    Array(std::size_t);
    std::size_t get_size() const;

private:
    std::size_t size;
};

template<class T>
Array<T>::Array(std::size_t size)
: BasicArray<T>(size)
, size(size) {}

template<class T>
std::size_t Array<T>::get_size() const {
    return size;
}

static void __attribute__((noinline)) test_matrix(const Matrix<int> &m) {
    for (std::size_t i = 0; i < m.get_size(1); ++i) {
        for (std::size_t j = 0; j < m.get_size(2); ++j) {
            static_cast<volatile void>(m(i, j) = i + j);
        }
    }
}

static void __attribute__((noinline))
test_array_array(const Array<Array<int>> &aa) {
    for (std::size_t i = 0; i < aa.get_size(); ++i) {
        for (std::size_t j = 0; j < aa[0].get_size(); ++j) {
            static_cast<volatile void>(aa[i][j] = i + j);
        }
    }
}

int main() {
    constexpr int N = 1000;
    Matrix<int> m(N, N);
    Array<Array<int>> aa(N);
    for (std::size_t i = 0; i < aa.get_size(); ++i) {
        aa[i] = Array<int>(N);
    }
    test_matrix(m);
    test_array_array(aa);
}

3
你所做的和许多人一样:猜测你应该担心什么。相反,让程序告诉你应该关注什么——它可能完全不同于你想象的。这里有一个例子 - Mike Dunlavey
1
也许你想使用不同于按存储顺序递增索引的访问模式进行基准测试,这样可以让优化编译器删除乘法和大部分间接引用。(还有@MikeDunlavey所说的。) - rici
1
@xiver77:无论如何,猜测应该优化什么都无法与良好的诊断相比。如果您是正确的,它会告诉您。如果是其他原因,它也会告诉您。如果是其他问题,并且您修复了它,那么可能数组索引是下一个最大的问题,也可能是以后的事情。期望在新编写的程序中仅获得10%就是自我限制,在猜测问题所在和知道实际问题之间存在巨大差异。 - Mike Dunlavey
1
如果你交换main函数中的最后两行,你会得到不同的结果吗?这可能只是缓存和内存碎片化等问题。每个测试应该在一个程序中被多次调用,交替使用这两种类型。 - Aaron McDaid
1
@AaronMcDaid 是的!我在for循环中运行了相同的测试10次,一个非常有趣的事实是,在第一次迭代中,矩阵版本总是比较慢的,但之后矩阵版本确实像预期的那样运行得更快。我从未想过顺序会导致速度差异。知道这个真的很有趣。谢谢! - user3810155
显示剩余10条评论
3个回答

6
两种方法的性能几乎相同,因为在两种情况下,最内层循环可以采用相同的优化方式,并且计算很可能受到内存限制。这意味着间接寻址的开销在其余的计算中被稀释,而这些计算占据了大部分时间,并且受到实际上可能比开销更大的变化的影响。因此,基准测试对于两种方法之间的差异不是非常敏感。以下是最内层循环的汇编代码(左侧:矩阵,右侧:数组的数组):
.LBB0_17:                                            .LBB1_30:
        movdqa  xmm5, xmm1                                   movdqa  xmm5, xmm1
        paddq   xmm5, xmm4                                   paddq   xmm5, xmm4
        movdqa  xmm6, xmm0                                   movdqa  xmm6, xmm0
        paddq   xmm6, xmm4                                   paddq   xmm6, xmm4
        shufps  xmm5, xmm6, 136                              shufps  xmm5, xmm6, 136 
        movdqa  xmm6, xmm3                                   movdqa  xmm6, xmm3
        paddq   xmm6, xmm1                                   paddq   xmm6, xmm1
        movdqa  xmm7, xmm3                                   movdqa  xmm7, xmm3
        paddq   xmm7, xmm0                                   paddq   xmm7, xmm0
        shufps  xmm6, xmm7, 136                              shufps  xmm6, xmm7, 136 
        movups  xmmword ptr [rdi + 4*rbx - 48], xmm5         movups  xmmword ptr [rsi + 4*rcx], xmm5
        movups  xmmword ptr [rdi + 4*rbx - 32], xmm6         movups  xmmword ptr [rsi + 4*rcx + 16], xmm6
        movdqa  xmm5, xmm0                                   movdqa  xmm5, xmm0
        paddq   xmm5, xmm10                                  paddq   xmm5, xmm10
        movdqa  xmm6, xmm1                                   movdqa  xmm6, xmm1
        paddq   xmm6, xmm10                                  paddq   xmm6, xmm10
        movdqa  xmm7, xmm3                                   movdqa  xmm7, xmm3
        paddq   xmm7, xmm6                                   paddq   xmm7, xmm6
        paddq   xmm6, xmm4                                   paddq   xmm6, xmm4
        movdqa  xmm2, xmm3                                   movdqa  xmm2, xmm3
        paddq   xmm2, xmm5                                   paddq   xmm2, xmm5
        paddq   xmm5, xmm4                                   paddq   xmm5, xmm4
        shufps  xmm6, xmm5, 136                              shufps  xmm6, xmm5, 136 
        shufps  xmm7, xmm2, 136                              shufps  xmm7, xmm2, 136 
        movups  xmmword ptr [rdi + 4*rbx - 16], xmm6         movups  xmmword ptr [rsi + 4*rcx + 32], xmm6
        movups  xmmword ptr [rdi + 4*rbx], xmm7              movups  xmmword ptr [rsi + 4*rcx + 48], xmm7
        add     rbx, 16                                      add     rcx, 16
        paddq   xmm1, xmm11                                  paddq   xmm1, xmm11
        paddq   xmm0, xmm11                                  paddq   xmm0, xmm11
        add     rbp, 2                                       add     rax, 2
        jne     .LBB0_17                                     jne     .LBB1_30

正如我们所看到的,这两种方法的循环基本上包含相同的指令。存储的顺序(movups)不同,但这不应该影响执行时间(特别是如果数组在内存中对齐)。不同寄存器名称也适用于同样的情况。循环使用SIMD指令(SSE)进行向量化,并展开4次,因此它可以非常快速(每个SIMD单元可以计算4个项目,每次迭代可以计算16个项目)。需要大约62个迭代才能完成最内层循环。
话虽如此,在这两种情况下,循环都会写入4*1000*1000=3.81 MiB的数据。这通常适用于较新处理器上的L3缓存(或旧处理器上的RAM)。L3 / RAM的吞吐量受到核心的限制(远低于L1甚至L2缓存),因此一个核心可能会等待内存层次结构准备就绪而停滞不前。因此,循环不是很快,因为它们大部分时间都在等待内存层次结构。硬件预取程序在现代x86-64处理器上非常有效,因此它们可以在核心实际请求数据之前预取数据,特别是对于存储和如果编写的数据是连续的。
数组的数组方法通常效率较低,因为每个子数组不能保证连续分配。现代内存分配器通常使用基于桶的策略来查找适合请求大小的内存块。在像这个基准测试程序中,请求的内存可以是连续的(或非常接近),因为所有数组都是原始分配的,而且当程序启动时,桶内存通常不会碎片化。然而,当内存被分段时,数组往往位于非连续区域,导致一种称为内存扩散的效应。内存扩散使预取器难以有效地工作,从而导致更低效的加载/存储。这通常特别适用于加载,但由于写分配缓存策略,在大多数x86-64处理器(Intel处理器或最近的AMD处理器)上存储也会导致加载。简而言之,这是数组的数组方法通常在应用中效率较低的主要原因之一。但这不是唯一的原因:另一个原因来自间接引用。
在这个基准测试中,额外的间接引用开销非常小,主要是由于内存绑定的内部循环。子数组的指针被连续存储,因此它们可以适合L1缓存并被有效地预取。这意味着间接引用可以很快,因为它们不太可能导致缓存未命中。间接引用指令会导致额外的加载指令,但由于大多数时间都在等待L3缓存或RAM,因此这些指令的开销非常小,甚至可以忽略不计。事实上,现代处理器以并行方式乱序方式执行指令,因此L1访问可以与L3/RAM加载/存储重叠。例如,英特尔处理器专门为此提供了单元:线填充缓冲区(位于L1和L2缓存之间)、超级队列(位于L2和L3缓存之间)和集成内存控制器(位于L3和RAM之间)。大多数操作都是异步完成的。话虽如此,当核心停止等待传入数据或缓冲区/队列饱和时,事情开始变得同步。

如果内层循环较小或者2D数组的遍历是非连续的,这是可能的。确实,如果内层循环只计算少量项目,或者甚至用1个语句替换,那么间接寻址的开销就更加明显。处理器无法(轻松地)重叠开销和数组方法,使得数组方法比基于矩阵的方法更慢。 这里 是这个新基准测试的结果。两种方法之间的差距似乎很小,但应该记住,在基准测试期间,缓存是热的,而在实际应用中可能不是这样。拥有冷缓存有利于基于矩阵的方法,它需要从缓存中加载更少的数据(不需要加载指针数组)。

为了理解为什么差距并不是很大,我们需要再次分析汇编代码。完整的汇编代码可以在Godbolt上看到。Clang使用了三种不同的策略来加速循环(SIMD、标量+展开和标量),但展开的策略才应该在这种情况下使用。这里是基于矩阵的方法的热点循环部分:
.LBB0_27:
        mov     dword ptr [r12 + rdi], esi
        lea     ebx, [rsi + 1]
        mov     dword ptr [r12 + rdx], ebx
        lea     ebx, [rsi + 2]
        mov     dword ptr [r12 + rcx], ebx
        lea     ebx, [rsi + 3]
        mov     dword ptr [r12 + rax], ebx
        add     rsi, 4
        add     r12, r8
        cmp     rsi, r9
        jne     .LBB0_27

这是一个关于二维数组的例子:

.LBB1_28:
        mov     rbp, qword ptr [rdi - 48]
        mov     dword ptr [rbp], edx
        mov     rbp, qword ptr [rdi - 32]
        lea     ebx, [rdx + 1]
        mov     dword ptr [rbp], ebx
        mov     rbp, qword ptr [rdi - 16]
        lea     ebx, [rdx + 2]
        mov     dword ptr [rbp], ebx
        mov     rbp, qword ptr [rdi]
        lea     ebx, [rdx + 3]
        mov     dword ptr [rbp], ebx
        add     rdx, 4
        add     rdi, 64
        cmp     rsi, rdx
        jne     .LBB1_28

乍一看,第二个方法的效率显然较低,因为需要执行更多的指令。但正如先前所述,现代处理器可以并行执行指令。因此,指令依赖性,特别是关键路径,在结果性能中起着重要作用(例如依赖链),更不用说处理器单元的饱和度以及计算核心后端端口的饱和度。由于这个循环的性能强烈依赖于目标架构,我们应该考虑选择一个特定的处理器架构,以便分析每种方法在这种情况下的速度有多快。让我们选择一个相对较新的主流架构:Intel CoffeeLake。
第一个循环显然由存储指令 (mov dword ptr [...], ...) 限制,因为该架构只有1个存储端口,而 leaadd 指令可以在多个端口上执行 (而且 cmp+jne 很便宜,因为它可以 宏融合预测)。除非受内存层次结构的限制,否则每次迭代循环应该需要 4个周期
第二个循环更为复杂,但也由存储指令mov dword ptr [rbp], edx限制。实际上,CofeeLake有两个加载端口,因此每个周期可以执行2个mov rbp, qword ptr [...]指令;同样的情况也适用于可以在2个端口上执行的leaaddcmp+jne仍然很便宜。指令数量不足以使前端饱和,因此端口是瓶颈。最终,假设内存层次结构没有问题,该循环每次迭代还需要4个周期。问题在于,在实践中,指令的调度并不总是完美的,因此对于加载指令的依赖如果出现问题,可能会引入显著的延迟。由于内存层次结构的压力更大,缓存未命中将导致第二个循环停滞多个周期,而第一个循环(仅进行写操作)则不会如此。更何况,由于存在一个8KB的指针缓存区,以便在L1缓存中保持这个计算快速进行,所以在第二种情况下发生缓存未命中的可能性更大:从L2加载项目需要几个周期,并将数据加载到L3可能会导致某些缓存行被驱逐。这就是为什么第二个循环在这个新基准测试中略微较慢的原因。

如果我们使用另一个处理器会怎样?结果可能会有很大差异,特别是对于IceLake(Intel)和Zen2(AMD),因为它们有2个存储端口。在这种处理器上分析问题非常困难(因为不止一个端口可能是瓶颈,实际上后端可能根本不是瓶颈)。对于Zen2/Zen3来说更是如此,因为它们有2个共享的load/store端口和1个仅用于stores的专用端口(意味着每个周期可以调度2个loads+1个store,或1个load+2个stores,或没有load+3个stores)。因此,最好是在这些平台上运行实际的基准测试,并注意避免基准测试偏差。

注意,子数组的内存对齐也非常关键。当N = 1024时,基于矩阵的方法可能会明显变慢。这是因为在这种情况下,基于矩阵的方法的内存布局很可能会导致高速缓存清除,而基于数组的方法通常会添加一些填充来防止在这种非常特定的情况下出现此问题。事实上,对于主流基于桶的分配器,添加的填充通常为sizeof(size_t),因此该问题仅会在N的另一个值上发生,并没有真正预防。实际上,对于N = 1022,基于数组的方法明显变慢。这与上述解释完全相符,因为在目标机器(64位)上sizeof(size_t) = 2*sizeof(int) = 8。因此,两种方法都存在这个问题,但是可以通过向基于矩阵的方法添加一些填充来轻松控制,而使用基于数组的方法则无法轻松控制,因为分配器的实现默认依赖于平台。

这能解释为什么关闭优化时时间相同吗?(请参见我的非答案)。 - alfC
1
当关闭优化时,代码似乎在两种情况下大部分时间都在调用函数(这非常低效且相当大)。例如,get_sizeoperator[]operator()在每次迭代中都会被调用。尽管如此,矩阵类的开销并不像数组-数组那样大,因为后者必须调用更多的函数,这些函数效率较低(operator[]被调用两次,其内部代码不太高效,因为它应该被优化)。这在快速运行中得到了证实。 - Jérôme Richard
1
因此,简而言之:我无法重现“当关闭优化时时间相同”的事实,而我实际上证明了相反的情况。 - Jérôme Richard
你是对的,我读错了结果。矩阵比数组更适合(非)优化。 - alfC

1
我没有详细查看您的代码。相反,我测试了您的实现,针对一个非常简单的围绕std::vector的包装器进行了测试,然后添加了一些时间代码,以便我不必在分析器下运行即可获得有意义的结果。哦,而且我真的不喜欢代码采用对const的引用,然后使用强制转换为void来允许代码修改矩阵。我肯定无法想象在正常使用中期望人们这样做。
结果如下:
#include <chrono>
#include <iomanip>
#include <iostream>
#include <memory>
#include <vector>

template <class T>
class BasicArray : public std::unique_ptr<T[]> {
public:
    BasicArray() = default;
    BasicArray(std::size_t);
};

template <class T>
BasicArray<T>::BasicArray(std::size_t size)
    : std::unique_ptr<T[]>(new T[size])
{
}

template <class T>
class Matrix : public BasicArray<T> {
public:
    Matrix() = default;
    Matrix(std::size_t, std::size_t);
    T& operator()(std::size_t, std::size_t) const;
    std::size_t get_index(std::size_t, std::size_t) const;
    std::size_t get_size(std::size_t) const;

private:
    std::size_t sizes[2];
};

template <class T>
Matrix<T>::Matrix(std::size_t i, std::size_t j)
    : BasicArray<T>(i * j)
    , sizes { i, j }
{
}

template <class T>
T& Matrix<T>::operator()(std::size_t i, std::size_t j) const
{
    return (*this)[get_index(i, j)];
}

template <class T>
std::size_t Matrix<T>::get_index(std::size_t i, std::size_t j) const
{
    return i * get_size(2) + j;
}

template <class T>
std::size_t Matrix<T>::get_size(std::size_t d) const
{
    return sizes[d - 1];
}

template <class T>
class Array : public BasicArray<T> {
public:
    Array() = default;
    Array(std::size_t);
    std::size_t get_size() const;

private:
    std::size_t size;
};

template <class T>
Array<T>::Array(std::size_t size)
    : BasicArray<T>(size)
    , size(size)
{
}

template <class T>
std::size_t Array<T>::get_size() const
{
    return size;
}

static void test_matrix(Matrix<int>& m)
{
    for (std::size_t i = 0; i < m.get_size(1); ++i) {
        for (std::size_t j = 0; j < m.get_size(2); ++j) {
            m(i, j) = i + j;
        }
    }
}

static void
test_array_array(Array<Array<int>>& aa)
{
    for (std::size_t i = 0; i < aa.get_size(); ++i) {
        for (std::size_t j = 0; j < aa[0].get_size(); ++j) {
            aa[i][j] = i + j;
        }
    }
}

namespace JVC {
template <class T>
class matrix {
    std::vector<T> data;
    size_t cols;
    size_t rows;

public:
    matrix(size_t y, size_t x)
        : cols(x)
        , rows(y)
        , data(x * y)
    {
    }
    T& operator()(size_t y, size_t x)
    {
        return data[y * cols + x];
    }

    T operator()(size_t y, size_t x) const
    {
        return data[y * cols + x];
    }
    std::size_t get_rows() const { return rows; }
    std::size_t get_cols() const { return cols; }
};

static void test_matrix(matrix<int>& m)
{
    for (std::size_t i = 0; i < m.get_rows(); ++i) {
        for (std::size_t j = 0; j < m.get_cols(); ++j) {
            m(i, j) = i + j;
        }
    }
}
}

template <class F, class C>
void do_test(F f, C &c, std::string const &label) {
    using namespace std::chrono;

    auto start = high_resolution_clock::now();
    f(c);
    auto stop = high_resolution_clock::now();

    std::cout << std::setw(20) << label << " time: ";
    std::cout << duration_cast<milliseconds>(stop - start).count() << " ms\n";
}

int main()
{
    std::cout.imbue(std::locale(""));

    constexpr int N = 20000;
    Matrix<int> m(N, N);
    Array<Array<int>> aa(N);
    JVC::matrix<int> m2 { N, N };

    for (std::size_t i = 0; i < aa.get_size(); ++i) {
        aa[i] = Array<int>(N);
    }
    using namespace std::chrono;

    do_test(test_matrix, m, "Matrix");
    do_test(test_array_array, aa, "array of arrays");
    do_test(JVC::test_matrix, m2, "JVC Matrix");
}

结果看起来像这样:

              Matrix time: 1,893 ms
     array of arrays time: 1,812 ms
          JVC Matrix time: 620 ms

因此,一个对于std::vector的简单封装比你两种实现都快大约3倍。

我建议,在这么大的开销下,很难确定你所看到的时间差异是否源自存储布局。


基准测试是有偏差的。编译器可以内联函数并完全优化掉一个函数。Clang在我的机器上实际上就这么做了(显然在GodBolt上也是如此)。这导致了“矩阵时间:0毫秒”,因此从技术上讲它是最快的;此外,基准测试只计算一次,因此受分配器和操作顺序的影响而产生变化。如果JVC时间是在我的机器上最后一个基准测试启动的话,那么它才会很快。如果我改变do_test的顺序,它就不再快了。这肯定是由于页面错误引起的。 - Jérôme Richard

0

令我惊讶的是,你的测试基本上是正确的。 它们也违反了历史知识。(参见C语言中的动态数组-错误方法)

我用Quickbench验证了结果,两个时间几乎相同。 https://quick-bench.com/q/FhhJTV8IdIym0rUMkbUxvgnXPeA

我别无选择,只能说由于你的代码如此规则,编译器正在确定你要求连续的等大小分配,这可以被一个单一的块所替代,进而,硬件可以预测访问模式。

然而,我尝试将N设置为易失性,并在初始化时插入一堆随机交错的分配,仍然得到相同的结果。 我甚至将优化降低到-Og-Ofast,并增加了N,但仍然得到相同的结果。

只有当我使用了benchmark::ClobberMemory时,我才看到了一个非常小但明显的差异(在clang中,但在GCC中没有)。 所以这可能与内存访问模式有关。

quickbench https://quick-bench.com/q/FhhJTV8IdIym0rUMkbUxvgnXPeA

另一个在实际应用中起到(小)差异但重要的事情是将初始化步骤包含在计时内部,但令人惊讶的是,它只有5%到10%的差距(对单个块数组更有利)。
结论:编译器,或者更可能是硬件,一定在做一些非常神奇的事情。 指针间接版本从未真正比块数组更快,这让我觉得某种效果将一个情况转化为另一个情况。 这值得进一步研究。 如果有人感兴趣,这里是机器代码https://godbolt.org/z/ssGj7aq7j
事后想法:在放弃连续数组之前,我至少会怀疑这个结果是否只适用于2维,并且对于3或4维的结构无效。
免责声明:对我来说,这很有趣,因为我正在实现一个多维数组库,我关心性能。 该库是您类的通用版本Matrix,适用于任意维度https://gitlab.com/correaa/boost-multi

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接