OpenMP 阵列性能不佳
OpenMP poor performance with arrays
我遇到以下问题：我正在尝试使用 openMP 在 C++ 中并行化一个非常简单的 PDE 求解器，但是如果我增加线程数，性能不会提高。该方程是一个包含对流的简单一维热方程。由于我在每个时间步都需要解决方案，因此我决定使用二维数组
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
分析
首先，您的工作粒度太小，多线程效率不高。
实际上，您的顺序时间是 9.6 毫秒，并且有 999 个时间步长。
因此，每个时间步大约需要 9.6 us，这是相当小的。
此外，内存访问效率不高：
- 一方面，使用
std::vector<std::vector<std::vector<double>>> 会在内部生成一个数组，该数组包含指向数组的指针，该数组包含指向双精度数组（全部动态分配）的指针。数组在内存中可能不连续，也可能对齐不当。这会显着减慢执行速度，因为处理器更难从内存中预取数据。考虑分配一个大数组而不是许多小数组（例如，一个大的展平 std::vector）。
- 另一方面，按照您的方式使用填充会导致非常低效的内存访问模式。实际上，您只使用 1 倍于 8，所以 7/8 的内存使用量被浪费了（可能更多，因为 std::vector 可以分配额外的内存）。此外，由于添加了填充， read/written 不是连续的，处理器很难预取数据并有效地使用内存（因为 reads/writes 是按缓存行执行的，即多个双精度标量）。考虑应用矩阵的填充行或块（不是标量）。
最后，使用块大小为 8 的时间表似乎太小了。对于 parallel-for 和 reduction，只指定 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，这是非常小的）。
我遇到以下问题：我正在尝试使用 openMP 在 C++ 中并行化一个非常简单的 PDE 求解器，但是如果我增加线程数，性能不会提高。该方程是一个包含对流的简单一维热方程。由于我在每个时间步都需要解决方案，因此我决定使用二维数组
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
分析
首先，您的工作粒度太小，多线程效率不高。 实际上，您的顺序时间是 9.6 毫秒，并且有 999 个时间步长。 因此，每个时间步大约需要 9.6 us，这是相当小的。
此外，内存访问效率不高：
- 一方面，使用
std::vector<std::vector<std::vector<double>>>会在内部生成一个数组，该数组包含指向数组的指针，该数组包含指向双精度数组（全部动态分配）的指针。数组在内存中可能不连续，也可能对齐不当。这会显着减慢执行速度，因为处理器更难从内存中预取数据。考虑分配一个大数组而不是许多小数组（例如，一个大的展平 std::vector）。 - 另一方面，按照您的方式使用填充会导致非常低效的内存访问模式。实际上，您只使用 1 倍于 8，所以 7/8 的内存使用量被浪费了（可能更多，因为 std::vector 可以分配额外的内存）。此外，由于添加了填充， read/written 不是连续的，处理器很难预取数据并有效地使用内存（因为 reads/writes 是按缓存行执行的，即多个双精度标量）。考虑应用矩阵的填充行或块（不是标量）。
最后，使用块大小为 8 的时间表似乎太小了。对于 parallel-for 和 reduction，只指定 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，这是非常小的）。