c++ - 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>>>在内部生成一个数组，该数组包含指向数组的指针，该指针包含指向 double 数组的指针(均动态分配)。数组可能在内存中不连续，也可能对齐不良。由于处理器从内存中预取数据更加困难，因此这可能会大大降低执行速度。考虑分配一个大数组而不是许多小数组(例如，一个大的扁平std::vector)。

另一方面，使用填充方式会导致效率很低的内存访问模式。确实，您只使用了8的1倍，所以浪费了7/8的内存使用量(因为std::vector可以分配更多的内存，所以可能会浪费更多)。此外，由于增加了填充，一个读/写操作是不连续的，并且处理器很难预取数据并有效地使用内存(因为每个缓存行都执行读/写操作，即多个 double 标量)。考虑应用矩阵的填充行或块(非标量)。

最后，使用大小为8的块的计划似乎太小了。在此处，对于parallel-for和reduce，仅指定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，在此非常小)。