c - 使用 open mp 的慢速稀疏矩阵 vector 积 (CSR)

Question

我正在尝试使用 open mp 加速稀疏矩阵 vector 乘积，代码如下:

void zAx(double * z, double * data, long * colind, long * row_ptr, double * x, int M){

long i, j, ckey;
int chunk = 1000;
//int * counts[8]={0};
#pragma omp parallel num_threads(8)
{ 
  #pragma omp for private(ckey,j,i) schedule(static,chunk)
  for (i=0; i<M; i++ ){ 
    z[i]=0;
    for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
      j = colind[ckey];
      z[i] += data[ckey]*x[j];
    }              
  }
}
}

现在，这段代码运行良好，并产生了正确的结果，但它只使我的速度提高了约 30%。我检查过线程都获得了相同数量的非零元素(它们是)，并且矩阵相当大(300,000 x 300,000)，所以我希望开销不是唯一的问题.我也尝试过使用不同的 block 大小和线程数运行，我得到了类似的性能。

还有什么我可以尝试从中获得额外速度的东西吗？或者我显然做错了什么？

干杯。

编辑:刚刚注释掉'//int * counts[8]={0}'，因为它是计算工作分配的剩余部分。不需要

Edit2(更多细节):

好吧，我计算了一个调用这个 5000 次的循环并得到了平均次数:

• seq:0.0036(秒？)
• 2 个线程:0.002613
• 4 个线程:0.002308
• 8 个线程:0.002384

矩阵的大小为:303544x303544，并且有:2122980 个非零元素。

使用更小的矩阵 30000x30000 我得到的时间更像

• 序列 0.000303
• 8 个线程 0.000078

看来大尺寸可能是我的问题。

Answer 1

欢迎来到内存限制问题的奇妙世界。为了减轻您的痛苦，我想告诉您，稀疏矩阵 vector 乘法是许多无法在单个多核芯片上有效并行化甚至矢量化的事情之一，除非所有 数据可以放入最后一级缓存或内存总线真的很宽。

为什么？仅仅是因为计算与内存访问的比率非常低。对于内循环的每次迭代，您获取一次列索引到 j(8 字节)，矩阵元素到 data(8 字节)， vector 元素的值(8 字节)和结果的先前值(因为编译器很少优化对共享变量的访问)(8 字节)。然后执行 2 个非常快的浮点运算 (FLOP) 并执行存储(尽管 += 运算符被转换为单个指令，但它仍然是一个“获取-修改-写入”指令)。您总共加载了 32 个字节，并对它们执行了 2 次 FLOP。这使得每字节 1/16 FLOPs。

现代支持 SSE 的 CPU 内核可以执行 4 个 double FLOP/周期，这通常会导致每个 CPU 内核大约 8 GFLOPS(假设基本频率为 2 GHz)。使用 AVX，这个数字翻了一番，因此您可以在 2 GHz Intel Sandy/Ivy Bridge 或 AMD 同等处理器上获得高达 16 GFLOPS 的每个内核。为了使数据处理能力达到饱和，给定 1/16 FLOPs/字节，您至少需要 128 GiB/s 的内存带宽。

Xeon X7560 等高端 Nehalem-EX 处理器 运行频率为 2.26 GHz(9.04 GFLOPS/内核)及其共享的 L3 缓存(L1 和 L2 缓存是每个内核的)提供大约 275 GiB/s。在 9,04 GFLOPS/核心下，您需要每个核心 144,64 GiB/s 才能满足 zAx 例程的内部循环。这意味着在理想情况下，该 CPU 的 L3 缓存不能提供超过 2 个完全向量化的乘法内核。

在没有 SSE 向量化的情况下， double 的 FLOPS 速率要低两倍，因此可以预期问题会扩展到 4 个线程。一旦您的问题变得比 L3 高速缓存大，事情就会变得非常糟糕，因为内存总线提供的带宽少了大约十倍。

尝试使用以下版本的内部循环，看看编译器是否足够智能以遵循 OpenMP 的宽松内存 View :

#pragma omp for private(ckey,j) schedule(static,chunk)
for (i=0; i<M; i++){
  double zi = 0.0;
  for (ckey=row_ptr[i]; ckey<row_ptr[i+1]; ckey++) {
    j = colind[ckey];
    zi += data[ckey]*x[j];
  }
  z[i] = zi;
}

不幸的是，您无能为力。稀疏矩阵 vector 乘法与 CPU 插槽的数量成比例，而不是与 CPU 内核的数量成比例。您需要一个带有独立内存 Controller 的多路系统，例如具有多个(后)Nehalem 或 AMD64 处理器的任何系统。

编辑:优化提示。您真的需要 long 来存储列索引和行指针吗？对于 2122980 个非零元素，使用 int 会很好。它将在内循环中为每个元素节省加载 4 个字节，在外循环中为每行节省另外 4 个字节。