为什么这段代码不能线性扩展？

Question

我编写了这个 SOR 求解器代码。不要太在意这个算法做了什么，这不是这里的问题。但只是为了完整起见：它可能会求解线性方程组，具体取决于系统的条件。

我使用病态的 2097152 行空间矩阵（从不收敛）运行它，每行最多有 7 个非零列。

翻译：外部do-while 循环将执行10000 次迭代（我传递的值为max_iters），中间for 将执行2097152 次迭代，分成work_line 的块，在OpenMP 线程之间划分。最里面的for 循环将有 7 次迭代，除非在极少数情况下（小于 1%）可以更少。

sol 数组的值中的线程之间存在数据依赖关系。中间for 的每次迭代都会更新一个元素，但最多会读取数组的其他 6 个元素。由于 SOR 不是一个精确的算法，因此在读取时，它可以具有该位置上的任何先前或当前值（如果您熟悉求解器，这是一个 Gauss-Siedel，它在某些地方容忍 Jacobi 行为，以便并行性）。

typedef struct{
    size_t size;

    unsigned int *col_buffer;
    unsigned int *row_jumper;
    real *elements;
} Mat;

int work_line;

// Assumes there are no null elements on main diagonal
unsigned int solve(const Mat* matrix, const real *rhs, real *sol, real sor_omega, unsigned int max_iters, real tolerance)
{
    real *coefs = matrix->elements;
    unsigned int *cols = matrix->col_buffer;
    unsigned int *rows = matrix->row_jumper;
    int size = matrix->size;
    real compl_omega = 1.0 - sor_omega;
    unsigned int count = 0;
    bool done;

    do {
        done = true;
        #pragma omp parallel shared(done)
        {
            bool tdone = true;

            #pragma omp for nowait schedule(dynamic, work_line)
            for(int i = 0; i < size; ++i) {
                real new_val = rhs[i];
                real diagonal;
                real residual;
                unsigned int end = rows[i+1];

                for(int j = rows[i]; j < end; ++j) {
                    unsigned int col = cols[j];
                    if(col != i) {
                        real tmp;
                        #pragma omp atomic read
                        tmp = sol[col];

                        new_val -= coefs[j] * tmp;
                    } else {
                        diagonal = coefs[j];
                    }
                }

                residual = fabs(new_val - diagonal * sol[i]);
                if(residual > tolerance) {
                    tdone = false;
                }

                new_val = sor_omega * new_val / diagonal + compl_omega * sol[i];
                #pragma omp atomic write
                sol[i] = new_val;
            }

            #pragma omp atomic update
            done &= tdone;
        }
    } while(++count < max_iters && !done);

    return count;
}

正如你所看到的，并行区域内没有锁，所以，对于他们总是教我们的，这是一种 100% 并行的问题。这不是我在实践中看到的。

我的所有测试均在 Intel(R) Xeon(R) CPU E5-2670 v2 @ 2.50GHz、2 个处理器、每个 10 个内核、启用超线程、总计多达 40 个逻辑内核上运行。

在我的第一组运行中，work_line 固定为 2048，线程数从 1 到 40 不等（总共 40 次运行）。这是每次运行的执行时间（秒 x 线程数）的图表：

惊喜的是对数曲线，所以我想既然工作线这么大，共享缓存没有很好地使用，所以我挖了这个虚拟文件/sys/devices/system/cpu/cpu0/cache/index0/coherency_line_size告诉我这个处理器的L1缓存同步更新以 64 个字节为一组（数组 sol 中有 8 个双精度数）。所以我将work_line 设置为8：

然后我认为 8 太低，无法避免 NUMA 停顿并将work_line 设置为 16：

在运行上述程序时，我想“我是谁来预测 work_line 是好的？让我们看看...”，并计划从 8 到 2048 运行每个 work_line，步长为 8（即每缓存行的倍数，从 1 到 256）。 20 和 40 线程的结果（秒 x 中间for 循环的分割大小，在线程之间划分）：

我相信work_line 较低的情况会严重影响缓存同步，而较大的work_line 在一定数量的线程之外没有任何好处（我假设是因为内存路径是瓶颈）。很遗憾，一个看似 100% 并行的问题在真机上却表现出如此糟糕的行为。所以，在我确信多核系统是一个非常畅销的谎言之前，我先在这里问你：

如何使此代码与内核数量成线性关系？我错过了什么？问题中是否有什么东西使它不如最初看起来的那么好？

更新

根据建议，我使用static 和dynamic 调度进行了测试，但删除了数组sol 上的原子读/写。作为参考，蓝色和橙色线与上图相同（最多为work_line = 248;）。黄线和绿线是新线。对于我所看到的：static 对低 work_line 产生了显着影响，但在 96 之后，dynamic 的好处超过了它的开销，使其更快。原子操作根本没有区别。

Answer 1

您的内部循环有一个omp atomic read，而您的中间循环有一个omp atomic write，该位置可能与其中一个读取的位置相同。 OpenMP 有义务确保对同一位置的原子写入和读取进行序列化，因此实际上它可能确实需要引入一个锁，即使没有任何明确的锁。

它甚至可能需要锁定整个 sol 数组，除非它能够以某种方式确定哪些读取可能与哪些写入发生冲突，而且实际上，OpenMP 处理器并不一定那么聪明。

没有代码绝对线性扩展，但请放心，有许多代码确实比您的代码更接近线性扩展。

Answer 2

我怀疑您遇到了缓存问题。当一个线程更新sol 数组中的值时，它会使存储同一高速缓存行的其他 CPU 上的高速缓存无效。这会强制更新缓存，然后导致 CPU 停止。

Answer 3

尝试运行 IPCM (Intel Performance Counter Monitor)。您可以查看内存带宽，看看它是否会因更多内核而最大化。我的直觉是你的内存带宽有限。

作为信封计算的快速回溯，我发现 Xeon 上的未缓存读取带宽约为 10 GB/s。如果您的时钟为 2.5 GHz，则每个时钟周期一个 32 位字。您的内部循环基本上只是一个多重加法操作，您可以一方面计算其周期，再加上几个周期作为循环开销。在 10 个线程之后，您没有获得任何性能提升，这并不让我感到惊讶。

Answer 4

即使您的代码中没有明确的互斥锁，您的进程之间也有一个共享资源：内存及其总线。您在代码中看不到这一点，因为它是负责处理来自 CPU 的所有不同请求的硬件，但无论如何，它是一个共享资源。

因此，每当您的一个进程写入内存时，所有其他使用它的进程都必须从主内存中重新加载该内存位置，并且它们都必须使用相同的内存总线才能这样做。内存总线饱和，您无法从额外的 CPU 内核中获得更多性能提升，这只会使情况变得更糟。

Answer 5

稀疏矩阵向量乘法受内存限制（请参阅here），可以使用简单的屋顶线模型来显示。内存绑定问题受益于多插槽 NUMA 系统的更高内存带宽，但前提是数据初始化是以数据分布在两个 NUMA 域之间的方式完成的。我有一些理由相信您正在串行加载矩阵，因此它的所有内存都分配在单个 NUMA 节点上。在这种情况下，您将无法从双插槽系统上可用的双倍内存带宽中受益，并且使用schedule(dynamic) 或schedule(static) 真的没关系。您可以做的是启用内存交错 NUMA 策略，以便在两个 NUMA 节点之间分布内存分配。因此，每个线程最终将获得 50% 的本地内存访问和 50% 的远程内存访问，而不是让第二个 CPU 上的所有线程都受到 100% 的远程内存访问。启用该策略的最简单方法是使用numactl：

$ OMP_NUM_THREADS=... OMP_PROC_BIND=1 numactl --interleave=all ./program ...

OMP_PROC_BIND=1 启用线程锁定，应该会稍微提高性能。

我还想指出这一点：

done = true;
#pragma omp parallel shared(done)
{
    bool tdone = true;

    // ...

    #pragma omp atomic update
    done &= tdone;
}

可能是一个不太有效的重新实现：

done = true;
#pragma omp parallel reduction(&:done)
{
    // ...
        if(residual > tolerance) {
            done = false;
        }
    // ...
}

由于内部循环中完成的工作量，这两种实现之间不会有显着的性能差异，但为了可移植性和可读性重新实现现有的 OpenMP 原语仍然不是一个好主意。