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

Question

我写了这个SOR求解器代码.不要太费心这个算法做什么,这不是关注点.但仅仅是为了完整性:它可以解决线性方程组,这取决于系统的条件有多好.

我用一个病态的2097152行sparce矩阵(从不收敛)运行它,每行最多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个字节组的更新(数组中有8个双打sol).所以我设置work_line为8:

然后我认为8太低了,以至于无法避免NUMA档位并设置work_line为16:

在运行上述操作时,我想"我是谁来预测什么work_line是好的？让我们看看......",并计划每次运行work_line8到2048,步长为8(即每个缓存行的倍数,从1到1) 256).20和40个线程的结果(中间for循环分割的秒数x大小,在线程之间划分):

我认为低的情况work_line会因缓存同步而受到严重影响,而较大的情况work_line除了一定数量的线程之外没有任何好处(我假设因为内存路径是瓶颈).令人遗憾的是,似乎100%并行的问题在真实机器上呈现出这种不良行为.所以,在我确信多核系统是一个非常畅销的谎言之前,我先问你这里:

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

更新

根据建议,我测试了static和dynamic调度,但删除了数组上的原子读/写sol.作为参考,蓝色和橙色线与上图相同(仅适用于work_line = 248;).黄色和绿色线是新的.对于我所看到的:static对低值产生重大影响work_line,但在96之后,其好处dynamic超过了其开销,使其更快.原子操作完全没有区别.

Answer 1

稀疏矩阵向量乘法是内存绑定(见此处),它可以用简单的车顶线模型显示.内存限制问题受益于多插槽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原语仍然不是一个好主意.

Answer 2

尝试运行IPCM(英特尔性能计数器监视器).您可以观察内存带宽,并查看是否有更多内核.我的直觉是你的内存带宽有限.

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