为什么C++动态数组的乘法比std::vector版本更好？

Question

我正在实现使用不同数据结构和技术（向量、数组和OpenMP）进行矩阵乘法的C++程序，并发现了一个奇怪的情况...我的动态数组版本效果更好：

时间如下:

openmp mult_1: time: 5.882000 s

array mult_2: time: 1.478000 s

我的编译标志是：

/usr/bin/g++ -fopenmp -pthread -std=c++1y -O3

C++向量版本

typedef std::vector<std::vector<float>> matrix_f;
void mult_1 (const matrix_f &  matrixOne, const matrix_f & matrixTwo, matrix_f & result) {
    const int matrixSize = (int)result.size();
    #pragma omp parallel for simd
    for (int rowResult = 0; rowResult < matrixSize; ++rowResult) {
        for (int colResult = 0; colResult < matrixSize; ++colResult) {
            for (int k = 0; k < matrixSize; ++k) {
                result[rowResult][colResult] += matrixOne[rowResult][k] * matrixTwo[k][colResult];  
            }
        }
    }
}

动态数组版本

void mult_2 ( float *  matrixOne, float * matrixTwo,  float * result, int size)  {
    for (int row = 0; row < size; ++row) {
        for (int col = 0; col < size; ++col) {
            for (int k = 0; k < size; ++k) {
                (*(result+(size*row)+col)) += (*(matrixOne+(size*row)+k)) * (*(matrixTwo+(size*k)+col));
            }
        }
    }
}

测试：

C++向量版本

utils::ChronoTimer timer;
/* set Up simple matrix */
utils::matrix::matrix_f matr1 = std::vector<std::vector<float>>(size,std::vector<float>(size));
fillRandomMatrix(matr1);

utils::matrix::matrix_f matr2 = std::vector<std::vector<float>>(size,std::vector<float>(size));
fillRandomMatrix(matr2);

utils::matrix::matrix_f result = std::vector<std::vector<float>>(size,std::vector<float>(size));    
timer.init();
utils::matrix::mult_1(matr1,matr2,result);
std::printf("openmp mult_1: time: %f ms\n",timer.now() / 1000);

动态数组版本

utils::ChronoTimer timer;

float *p_matr1 = new float[size*size];
float *p_matr2 = new float[size*size];
float *p_result = new float[size*size];

fillRandomMatrixArray(p_matr1,size);
fillRandomMatrixArray(p_matr2,size);

timer.init();
utils::matrix::mult_2(p_matr1,p_matr2,p_result,size);
std::printf("array mult_2: time: %f ms\n",timer.now() / 1000);

delete [] p_matr1;
delete [] p_matr2;
delete [] p_result;

我检查了一些以前的帖子，但我没有找到任何与我的问题相关的内容链接，链接2，链接3： 更新： 我根据答案重构了测试，vector 的效果略好于数组：

vector mult: 时间：1.194000 秒

array mult_2: 时间：1.202000 秒

C++ vector 版本

void mult (const std::vector<float> &  matrixOne, const std::vector<float> & matrixTwo, std::vector<float> & result, int size) {
    for (int row = 0; row < size; ++row) {
        for (int col = 0; col < size; ++col) {
            for (int k = 0; k <size; ++k) {
                result[(size*row)+col] += matrixOne[(size*row)+k] * matrixTwo[(size*k)+col];
            }
        }
    }
}

动态数组版本

void mult_2 ( float *  matrixOne, float * matrixTwo,  float * result, int size)  {
    for (int row = 0; row < size; ++row) {
        for (int col = 0; col < size; ++col) {
            for (int k = 0; k < size; ++k) {
                (*(result+(size*row)+col)) += (*(matrixOne+(size*row)+k)) * (*(matrixTwo+(size*k)+col));
            }
        }
    }
}

此外，我的向量化版本效果更好（0.803秒）;

- carlos.baez

---

数据在内存中的排列方式不同。您的矩阵在内存中是连续的，而使用vector<vector>时，每个向量单独分配。如果大小在编译时固定，则可以尝试vector<array<float,N>>或执行其他操作以确保完整的矩阵在内存中是连续的。 - PeterT

请参考https://dev59.com/IGQm5IYBdhLWcg3w4CTA，了解为什么通常应该避免使用“真正的”二维结构(例如T **，vector<vector<T>>...)来存储密集矩阵。 - Pixelchemist

我猜内存布局不是你唯一的问题。请向我们展示你的计时器代码以及你在OpenMP版本上运行了多少个线程。 - jepio

@jepio 我正在逐步应用每个改进步骤...我更改分配错误，检查线程并发布线程配置。 - carlos.baez

Answer 1

一个vector vector类似于一个嵌套数组--数组中每个条目都是指针，每个指针指向另一个浮点数数组。

相比之下，原始数组版本是一个内存块，在其中进行数学计算以找到元素。

使用单个vector而不是vector vector，并手动进行数学计算。或者，使用固定大小的std :: array的vector vector。或编写一个帮助器类型，该类型接受（一维）float向量，并为您提供二维视图。

在连续缓冲区中的数据比在其他断开的缓冲区中的数据更加缓存和优化友好。

template<std::size_t Dim, class T>
struct multi_dim_array_view_helper {
  std::size_t const* dims;
  T* t;
  std::size_t stride() const {
    return
      multi_dim_array_view_helper<Dim-1, T>{dims+1, nullptr}.stride()
      * *dims;
  }
  multi_dim_array_view_helper<Dim-1, T> operator[](std::size_t i)const{
    return {dims+1, t+i*stride()};
  }
};
template<class T>
struct multi_dim_array_view_helper<1, T> {
  std::size_t stride() const{ return 1; }
  T* t;
  T& operator[](std::size_t i)const{
    return t[i];
  }
  multi_dim_array_view_helper( std::size_t const*, T* p ):t(p) {}
};
template<std::size_t Dims>
using dims_t = std::array<std::size_t, Dims-1>;
template<std::size_t Dims, class T>
struct multi_dim_array_view_storage
{
  dims_t<Dims> storage;
};
template<std::size_t Dims, class T>
struct multi_dim_array_view:
  multi_dim_array_view_storage<Dims, T>,
  multi_dim_array_view_helper<Dims, T>
{
  multi_dim_array_view( dims_t<Dims> d, T* t ):
    multi_dim_array_view_storage<Dims, T>{std::move(d)},
    multi_dim_array_view_helper<Dims, T>{
      this->storage.data(), t
    }
  {}
};

现在您可以这样做：

std::vector<float> blah = {
   01.f, 02.f, 03.f,
   11.f, 12.f, 13.f,
   21.f, 22.f, 23.f,
};

multi_dim_array_view<2, float> view = { {3}, blah.data() };
for (std::size_t i = 0; i < 3; ++i )
{
  std::cout << "[";
  for (std::size_t j = 0; j < 3; ++j )
    std::cout << view[i][j] << ",";
  std::cout << "]\n";
}

实时示例

视图类中不复制任何数据。它只提供了一个多维数组的扁平数组的视图。

- Yakk - Adam Nevraumont

Answer 2

你的两种方法有所不同：

• 在“动态数组”版本中，你为每个矩阵分配一个单一的内存块，并将矩阵的行映射到这个一维内存范围。

• 在“向量”版本中，你使用“真实”和“动态”的二维向量，这意味着矩阵的每一行的存储位置与其他行无关。

你可能想要做的是：

• 使用vector<float>(size*size)，并手动执行与“动态数组”示例中相同的映射，或者

• 编写一个在内部处理映射并提供2维访问接口（T& operator()(size_t, size_t)或某种row_proxy operator[](size_t)，其中row_proxy反过来又有T& operator[](size_t)）的类。

- Pixelchemist

Answer 3

这只是为了实践有关连续内存的理论。

在对使用g++(-O2)生成的代码进行分析后，可以在以下链接中找到源代码：https://gist.github.com/42be237af8e3e2b1ca03

数组版本的相关生成代码如下：

.L3:
    lea r9, [r13+0+rbx]                ; <-------- KEEPS THE ADDRESS
    lea r11, [r12+rbx]
    xor edx, edx
.L7:
    lea r8, [rsi+rdx]
    movss   xmm1, DWORD PTR [r9]
    xor eax, eax
.L6:
    movss   xmm0, DWORD PTR [r11+rax*4]
    add rax, 1
    mulss   xmm0, DWORD PTR [r8]
    add r8, r10
    cmp ecx, eax
    addss   xmm1, xmm0
    movss   DWORD PTR [r9], xmm1     ; <------------ ADDRESS IS USED
    jg  .L6
    add rdx, 4
    add r9, 4                        ; <--- ADDRESS INCREMENTED WITH SIZE OF FLOAT
    cmp rdx, rdi
    jne .L7
    add ebp, 1
    add rbx, r10
    cmp ebp, ecx
    jne .L3

看看使用r9的值如何反映目标数组的连续内存，r8用于其中一个输入数组。

另一方面，向量的向量生成的代码如下：

.L12:
    mov r9, QWORD PTR [r12+r11]
    mov rdi, QWORD PTR [rbx+r11]
    xor ecx, ecx
.L16:
    movss   xmm1, DWORD PTR [rdi+rcx]
    mov rdx, r10
    xor eax, eax
    jmp .L15
.L13:
    movaps  xmm1, xmm0
.L15:
    mov rsi, QWORD PTR [rdx]
    movss   xmm0, DWORD PTR [r9+rax]
    add rax, 4
    add rdx, 24
    cmp r8, rax
    mulss   xmm0, DWORD PTR [rsi+rcx]
    addss   xmm0, xmm1
    movss   DWORD PTR [rdi+rcx], xmm0   ; <------------ HERE
    jne .L13
    add rcx, 4
    cmp rcx, r8
    jne .L16
    add r11, 24
    cmp r11, rbp
    jne .L12

不难想到，编译器足够聪明，不会为所有operator []调用生成代码，并且在内联它们方面做得很好，但是请看它在将值存回结果向量时如何通过rdi + rcx跟踪不同的地址，以及各个子向量的额外内存访问(mov rsi, QWORD PTR [rdx])，这些都会产生一些开销。

- Ferenc Deak