OpenMP 对数组的性能不佳

Question

我有以下问题：我正在尝试使用 openMP 在 c++ 中并行化一个非常简单的 PDE 求解器，但如果我增加线程数，性能并没有提高。该方程是一个简单的一维对流热方程。由于我在每个时间步都需要解决方案，因此我决定使用 2D 数组

double solution[iterationsTime][numPoints];

其中每一行都包含特定时间步的离散化函数。更新是通过 for 循环完成的

#pragma omp parallel default(shared) private(t, i, iBefore, iAfter)
{
for(t=0; t<iterationsTime; t++)

#pragma omp for schedule(auto) 
   for(i=0; i<numPoints; i++) {
       iBefore = (i==0)?numPoints-2:i-1;
       iAfter = (i==numPoints-1)?1:i;
       solution[t+1][i] = solution[t][iAfter] - solution[t][iBefore];
}

使用值 iBefore 和 iAfter 是因为我基本上将数组视为环形缓冲区，因此 PDE 具有周期性边界条件并且域被视为环形。无论如何，对 solution[t+1] 的每次更新都需要对 solution[t] 进行一些计算，就像上面代码中显示的那样。 我知道可扩展性差的原因很可能是错误共享，因此我将 2D 矩阵转换为 3D 矩阵

double solution[iterationsTime][numPoints][PAD];

这使我可以确保不会在共享缓存行上执行任何写操作，因为我可以改变 PAD 的大小。由于现在每个值都将存储在

中，因此代码有所更改

solution[t][i][0];

还有下面的

solution[t][i+1][0];

请注意，所需的内存是在并行区域之外使用 new 运算符在堆上分配的。该代码运行良好，但不可扩展。我尝试过不同的时间表，如静态、动态、自动……我用

编译它

g++ code.cpp -fopenmp -march=native -O3 -o out

我尝试删除或添加 -march 和 -O3 标志，但没有看到任何改进。我尝试了不同大小的 PAD 和环境变量，如 OMP_PROC_BIND 但没有改进。我不知道是什么导致了此时的性能损失。 这是代码

const int NX = 500; //DOMAIN DISCRETIZATION
const int PAD = 8; //PADDING TO AVOID FALSE SHARING
const double DX = 1.0/(NX-1.0); //STEP IN SPACE
const double DT = 0.01*DX; //STEP IN TIME
const int NT = 1000; //MAX TIME ITERATIONS
const double C = 10.0; //CONVECTION VELOCITY
const double K = 0.01; //DIFFUSION COEFFICIENT

int main(int argc, char **argv) {
  omp_set_num_threads(std::atoi(argv[1]));  //SET THE REQUIRED NUMBER OF THREADS

  //INTIALING MEMORY --> USING STD::VECTOR INSTEAD OF DOUBLE***
  std::vector<std::vector<std::vector<double>>> solution(NT, std::vector<std::vector<double>>(NX, std::vector<double>(PAD,0)));
  for (int i=0; i<NX; i++){
    solution[0][i][0] = std::sin(i*DX*2*M_PI); //INITIAL CONDITION
  }

  int numThreads, i, t, iBefore, iAfter;
  double energy[NT]{0.0}; //ENERGY of the solution --> e(t)= integral from 0 to 1 of ||u(x,t)||^2 dx

  //SOLVE THE PDE ON A RING
  double start = omp_get_wtime();
  #pragma omp parallel default(none) shared(solution, energy, numThreads, std::cout) private(i, t, iBefore, iAfter)
  {
    #pragma omp master
    numThreads = omp_get_num_threads();

    for(t=0; t<NT-1; t++){

      #pragma omp for schedule(static, 8) nowait
      for(i=0; i<NX; i++){
        iBefore = (i==0)?NX-2:i-1;
        iAfter = (i==NX-1)?1:i+1;
        solution[t+1][i][0]=solution[t][i][0] 
          + DT*(  -C*((solution[t][iAfter][0]-solution[t][iBefore][0])/(2*DX))
            + K*(solution[t][iAfter][0]-2*solution[t][i][0]+ solution[t][iBefore][0])/(DX*DX)  );
      }

      // COMPUTE THE ENERGY OF PREVOIUS TIME ITERATION
      #pragma omp for schedule(auto) reduction(+:energy[t])
      for(i=0; i<NX; i++) {
        energy[t] += DX*solution[t][i][0]*solution[t][i][0];
      }
    }
  }
  std::cout << "numThreads: " <<numThreads << ". Elapsed Time: "<<(omp_get_wtime()-start)*1000 << std::endl;
  return 0;
}

还有表演

numThreads: 1. Elapsed Time: 9.65456
numThreads: 2. Elapsed Time: 9.1855
numThreads: 3. Elapsed Time: 9.85965
numThreads: 4. Elapsed Time: 8.9077
numThreads: 5. Elapsed Time: 15.5986
numThreads: 6. Elapsed Time: 15.5627
numThreads: 7. Elapsed Time: 16.204
numThreads: 8. Elapsed Time: 17.5612

Answer 1

分析

首先，您正在处理的粒度太小，以至于多线程效率不高。 实际上，您的连续时间是 9.6 毫秒，并且有 999 个时间步长。 因此，每个时间步长大约需要 9.6 us，这相当小。

另外，内存访问效率不高：

• 一方面使用std::vector<std::vector<std::vector<double>>> 在内部生成一个数组，该数组包含指向数组的指针，该数组包含指向双精度数组的指针（全部动态分配）。数组在内存中可能不连续，也可能对齐不当。这会显着减慢执行速度，因为处理器更难以从内存中预取数据。考虑分配一个大数组而不是许多小数组（例如，一个大的扁平 std::vector）。
• 另一方面，使用填充的方式会导致非常低效的内存访问模式。实际上，您只使用 1 倍于 8，因此浪费了 7/8 的内存使用量（可能更多，因为 std::vector 可以分配额外的内存）。此外，由于添加了填充，读/写的一个不连续，处理器很难预取数据并有效地使用内存（因为读/写是按缓存行执行的，即多个双精度标量）。考虑应用填充行或矩阵块（不是标量）。

最后，使用块大小为 8 的时间表似乎太小了。在这里指定 schedule(static) 对于并行处理和减少可能会更好（如果您使用 nowait 并且我认为您想要正确的结果，那么两者的时间表应该是相同的和静态的）。

因此，您可能正在测量延迟和内存开销。

改进版

以下是修正后的代码，其中包含最重要的修复（错误共享效果被忽略）：

#include <iostream>
#include <vector>
#include <cmath>
#include <omp.h>

const int NX = 500; //DOMAIN DISCRETIZATION
const int PAD = 8; //PADDING TO AVOID FALSE SHARING
const double DX = 1.0/(NX-1.0); //STEP IN SPACE
const double DT = 0.01*DX; //STEP IN TIME
const int NT = 1000; //MAX TIME ITERATIONS
const double C = 10.0; //CONVECTION VELOCITY
const double K = 0.01; //DIFFUSION COEFFICIENT

int main(int argc, char **argv) {
  omp_set_num_threads(std::atoi(argv[1]));  //SET THE REQUIRED NUMBER OF THREADS

  //INTIALING MEMORY --> USING A FLATTEN DOUBLE ARRAY
  std::vector<double> solution(NT * NX);
  for (int i=0; i<NX; i++){
    solution[0*NX+i] = std::sin(i*DX*2*M_PI); //INITIAL CONDITION
  }

  int numThreads, i, t, iBefore, iAfter;
  double energy[NT]{0.0}; //ENERGY of the solution --> e(t)= integral from 0 to 1 of ||u(x,t)||^2 dx

  //SOLVE THE PDE ON A RING
  double start = omp_get_wtime();
  #pragma omp parallel default(none) shared(solution, energy, numThreads, std::cout) private(i, t, iBefore, iAfter)
  {
    #pragma omp master
    numThreads = omp_get_num_threads();

    for(t=0; t<NT-1; t++){

      #pragma omp for schedule(static) nowait
      for(i=0; i<NX; i++){
        iBefore = (i==0)?NX-2:i-1;
        iAfter = (i==NX-1)?1:i+1;
        solution[(t+1)*NX+i]=solution[t*NX+i]
          + DT*(  -C*((solution[t*NX+iAfter]-solution[t*NX+iBefore])/(2*DX))
            + K*(solution[t*NX+iAfter]-2*solution[t*NX+i]+ solution[t*NX+iBefore])/(DX*DX)  );
      }

      // COMPUTE THE ENERGY OF PREVOIUS TIME ITERATION
      #pragma omp for schedule(static) reduction(+:energy[t])
      for(i=0; i<NX; i++) {
        energy[t] += DX*solution[t*NX+i]*solution[t*NX+i];
      }
    }
  }
  std::cout << "numThreads: " <<numThreads << ". Elapsed Time: "<<(omp_get_wtime()-start)*1000 << std::endl;
  return 0;
}

结果

在我的 6 核 (Intel i5-9600KF) 机器上。我得到以下结果。

之前：

1 thread : 3.35 ms
2 threads: 2.90 ms
3 threads: 2.89 ms
4 threads: 2.83 ms
5 threads: 3.07 ms
6 threads: 2.90 ms

之后：

1 thread : 1.62 ms
2 threads: 1.03 ms
3 threads: 0.87 ms
4 threads: 0.95 ms
5 threads: 1.00 ms
6 threads: 1.16 ms

使用新版本，顺序时间更快，并且成功扩展到 3 个核心。然后同步开销变得很大，并减慢了整体执行速度（请注意，每个时间步长在这里持续不到 1 us，这是非常小的）。