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

Posted 2023-02-16

【中文标题】为啥这段代码不能线性扩展？

【英文标题】：Why doesn't this code scale linearly?为什么这段代码不能线性扩展？

【发布时间】：2014-12-01 03:17:50

【问题描述】：

我编写了这个 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 的好处超过了它的开销，使其更快。原子操作根本没有区别。

【问题讨论】：

我对 SOR/Gauss–Seidel 方法不太熟悉，但对于矩阵乘法或 Cholesky 分解，您获得良好缩放的唯一方法是使用循环平铺以便重用数据仍在缓存中。见***.com/questions/22479258/…。否则它会受到内存限制。

虽然我不熟悉该算法，但快速浏览一下该内部循环表明您可能有一些非常糟糕的空间内存局部性。 （通常是稀疏线性代数的情况）在这种情况下，您可能会受到内存访问的限制。

SOR 的时间复杂度是多少？ cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html#link_4 O(N^3/2)？使用 Matrix Mult，计算为 N^3，而读取为 N^2，这就是它可以很好地扩展的原因。因此，除非计算次数远大于读取次数，否则它将受到内存限制。如果忽略内核速度快而主内存慢的事实，许多基本算法似乎可以很好地扩展。 BLAS 2 级（例如 matrix*vec）可以很好地扩展而忽略慢速内存。只有 BLAS 级别 3 (O(N^3) 例如 GEMM、Choleksy、...) 才能很好地适应慢速内存。

Intel Linux 上的默认拓扑是分散的。这意味着在您的情况下，偶数线程对应于一个节点，奇数线程对应于另一个节点。我认为如果您尝试export GOMP_CPU_AFFINITY="0 2 4 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62" 和export OMP_NUM_THREADS=20，您的代码将在一个节点（一个套接字）上运行。

@Zboson，简称export GOMP_CPU_AFFINITY="0-62:2"。至于拓扑，核心编号由 Bios 设置，Linux 内核通过解析相应的 MP ACPI 表来找到它（MADT？不过我不会打赌）。 Bull 我们的大多数双插槽 Intel 机器都在单个封装中具有连续编号的内核。

【参考方案1】：

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

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

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

【讨论】：

我认为那里没有真正的软件锁。我没有看过程序集，但它们很可能是指令级别的原子读/写。无论如何，我将在没有原子读/写的情况下重新运行案例 3 的稀疏版本。对于更大的work_line，没有区别（我在另一台具有 4 个线程的机器上进行了测试），这是有道理的，因为不太可能发生冲突。对于较小的work_line，它可能是相关的。看到这个：gcc.gnu.org/onlinedocs/gcc-4.1.2/gcc/Atomic-Builtins.html

atomic read 和 atomic write 在 x86 上使用 lock 指令前缀实现，即没有重量级软件锁。

【参考方案2】：

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

【参考方案3】：

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

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

【讨论】：

我正在说服系统管理员允许我对/dev/cpu/*/msr...拥有读写权限...

这个算法其实众所周知是内存带宽有限的。

更不用说sol[col] 上的潜在缓存未命中只会让事情变得更糟。如果所有内核都已经停在内存上，那么这对 CPU 来说可能并不重要。但从带宽的角度来看，这样的缓存未命中会占用缓存行的带宽。

@VladimirF，我不怀疑该算法的 OP 实现受内存带宽限制，但您是否有关于该算法通常受内存限制的声明的来源？在cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html#link_4 有一些关于并行版本的讨论。我第一次实现 Cholesky 分解并没有很好地扩展，但经过深思熟虑后，我得到了很好的扩展。

嗯，我只从非常稀疏的矩阵中知道它，但是有很多关于循环平铺和 G.-S 的其他技巧的理论和文献。和 SOR 以提高缓存的重用性。使用它们是因为内存带宽限制。

【参考方案4】：

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

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

【参考方案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 原语仍然不是一个好主意。

【讨论】：

感谢您的提示。我只是在学习 OpenMP，无法理解减少的事情。

对numactl 产生了巨大影响。稍后我将使用 libnuma 在 NUMA 套接字之间正确拆分工作并相应地设置线程亲和性。

@lvella，您能否在使用numactl 后再次用结果更新您的问题？我很想看到结果。