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

Question

我编写了这个SOR求解器代码。不要过多关注算法的具体作用，这不是我们关心的重点。但为了完整起见，它可以解决一个线性方程组，取决于系统的条件。

我在一个病态的2097152行稀疏矩阵上运行它（它永远不会收敛），每行最多有7列非零元素。

简单来说，外部的do-while循环将执行10000次迭代（即我传递给它的max_iters值），中间的for循环将执行2097152次迭代，按照work_line分割成若干块，并在OpenMP线程之间进行划分。最内层的for循环将执行7次迭代，除了极少数情况（不到1%）外，它可能会少于7次。

在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个线程的结果（以中间的for循环的分裂大小为单位，平均分配给线程）：

我认为具有较小work_line的情况会因缓存同步而受到严重影响，而较大的work_line在一定数量的线程之外不提供任何好处（我假设是因为内存通道是瓶颈）。很遗憾，一个似乎100％并行的问题在实际机器上呈现出如此糟糕的行为。因此，在被说服多核系统是一个非常出色的谎言之前，我首先问你：

如何使此代码按线性比例扩展到核心数？我缺少什么？是否有什么东西使它不像一开始看起来那么好？

更新

根据建议，我尝试使用static和dynamic调度，但是去除了sol数组上的原子读/写操作。供参考，蓝色和橙色线与上一个图表相同（只到 work_line = 248;）。黄色和绿色线是新的。就我所看到的，

- lvella

---

我对SOR/Gauss-Seidel方法不是很熟悉，但是通过矩阵乘法或Cholesky分解，唯一能够获得良好缩放的方法就是使用循环分块以便在数据仍在缓存中时重复使用它。请参见https://dev59.com/8n3aa4cB1Zd3GeqPhLp_。否则，它将受到内存限制。 - Z boson

尽管我不熟悉这个算法，但快速浏览内部循环表明您可能有非常差的空间存储器局部性（这通常是稀疏线性代数的情况）。在这种情况下，您可能会受到内存访问限制。 - Mysticial

SOR的时间复杂度是多少？http://www.cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html#link_4 O(N^3/2)？矩阵乘法的计算复杂度为N^3，而读取的复杂度为N^2，因此它可以很好地扩展。因此，除非计算量比读取量大得多，否则它将受到内存限制。许多基本算法似乎可以很好地扩展，如果您忽略核心速度快而主内存速度慢的事实。 BLAS 2级（例如matrix*vec）会忽略缓慢的内存而扩展良好。只有BLAS 3级（O（N^3）例如GEMM，Choleksy，...）可以很好地扩展缓慢的内存。 - Z boson

Linux系统在使用Intel时采用的默认拓扑结构是分散的。这意味着在你的情况下，偶数线程对应一个节点，奇数线程对应另一个节点。我认为，如果你尝试使用“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`，你的代码将在一个节点上运行（即一个套接字）。 - Z boson

@Zboson，这是export GOMP_CPU_AFFINITY="0-62:2"的简写。至于拓扑结构，核心编号由BIOS设置，Linux内核通过解析相应的MP ACPI表（MADT？虽然我不确定）来找到它。Bull的大多数双插槽Intel机器都将核心放在一个单独的包中按顺序编号。 - Hristo Iliev

Answer 1

稀疏矩阵向量乘法受到内存限制（见此处），可以用简单的屋顶线模型来展示。内存限制问题受益于多插槽NUMA系统的更高内存带宽，但只有在数据初始化以使数据分布在两个NUMA域的情况下才能从中受益。我有一些理由相信您正在串行加载矩阵，因此所有的内存都分配在单个NUMA节点上。在这种情况下，您将无法从双插槽系统上可用的双倍内存带宽中受益，而且使用schedule(dynamic)或schedule(static)并没有任何意义。您可以启用内存交错NUMA策略，以便将内存分配在两个NUMA节点之间。因此，每个线程将最终具有50％的本地内存访问和50％的远程内存访问，而不是所有线程都受到100％远程内存访问的第二个CPU的影响。启用该策略的最简单方法是使用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基元并不是一个好主意。

- Hristo Iliev

---

谢谢你的建议。我正在学习OpenMP，对于reduction这个概念还有些困惑。 - lvella

使用 numactl 工具确实起到了很大的作用。稍后我会使用 libnuma 来正确地在 NUMA 插槽之间分配工作，并相应地设置线程亲和性。 - lvella

@lvella，请问您能否再次更新您的问题，并分享使用numactl后的结果呢？我非常期待看到结果。 - Z boson

Answer 2

尝试运行IPCM（Intel Performance Counter Monitor）以监视内存带宽，并查看是否随着更多核心而达到最大值。我的直觉是您的限制在于内存带宽。

作为一个快速的估算，我发现Xeon的未缓存读取带宽约为10 GB/s。如果您的时钟频率为2.5 GHz，则每个时钟周期一个32位字。您的内部循环基本上只是一个多重加法操作，其周期可以用一只手数出，再加上几个周期用于循环开销。在使用10个线程后，不会获得任何性能提升并不令我感到意外。

- Mark Lakata

---

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

这个算法实际上被广泛认为是受内存带宽限制的。 - Vladimir F Героям слава

我只了解在稀疏矩阵中的应用，但是有很多关于循环分块和其他技巧的理论和文献可以用于G.-S.和SOR算法以提高缓存重用。它们被使用是因为存在内存带宽限制。 - Vladimir F Героям слава

我无法在同一台Xeon机器上运行它，但在我的桌面机器上，L3缓存未命中率为32％，内存读取速率为15.41 GBytes / 3111 Mticks。Mticks如何转换为秒？那些是CPU时钟周期吗？ - lvella

我不确定，但我认为那是兆时钟周期。这将为您提供12 GB/s的速度，这大约是内存限制的极限......所以您肯定受到内存带宽的限制。 - Mark Lakata

Answer 3

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

它甚至可能需要锁定整个 sol 数组，除非它可以想出哪些读操作可能与哪些写操作冲突，而实际上，OpenMP 处理器并不一定那么智能。

没有代码完全按比例缩放，但请放心，有许多代码的缩放比你的更接近线性。

- John Bollinger

---

我认为那里没有真正的软件锁。我没有查看汇编代码，但它们很可能是原子读/写指令级别可用的。无论如何，我将重新运行一个更简洁的第3个案例，而不使用原子读/写。对于更大的“work_line”，这没有任何区别（我在另一台具有4个线程的机器上进行了测试），这是有道理的，因为冲突的可能性非常小。对于较小的“work_line”，这可能是相关的。请参见：https://gcc.gnu.org/onlinedocs/gcc-4.1.2/gcc/Atomic-Builtins.html - lvella

在x86上实现原子读取和原子写入是使用lock指令前缀，即没有重量级软件锁。 - Hristo Iliev

Answer 4

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

- dohashi

Answer 5

即使您的代码中没有显式的互斥锁，您的进程之间有一个共享资源：内存及其总线。您在代码中看不到这一点，因为是硬件负责处理来自CPU的所有不同请求，但尽管如此，它仍然是一个共享资源。因此，每当您的一个进程写入内存时，使用该内存位置的所有其他进程都必须从主内存重新加载它，并且它们都必须使用相同的内存总线来执行此操作。内存总线饱和，您无法从仅用于恶化情况的其他CPU核心获得更多性能提升。

- cmaster - reinstate monica