我正在尝试使用open mp加速稀疏矩阵向量产品,代码如下:
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),所以我希望开销不是唯一的问题.我也尝试使用不同的块大小和线程数运行,我得到了类似的性能.
还有什么我可以试着从中获得一点额外的速度吗？或者我显然做错了什么？
干杯.
编辑:刚刚注释掉'// int*counts [8] = {0}',因为它是计算工作分配的剩余部分.不需要
Edit2(更多细节):
好的,所以我定时拨打5000次循环并得到平均时间:
矩阵的大小为:303544x303544,具有:2122980非零元素.
使用更小的矩阵30000x30000,我得到更多的时间
所以似乎大尺寸可能是我的问题.
Hri*_*iev 13
欢迎来到记忆受限问题的精彩世界.为了减轻您的痛苦,我想告诉您,稀疏矩阵向量乘法是在单个多核芯片上无法有效并行化甚至矢量化的众多事物之一,除非所有数据都适合最后一个级别缓存或内存总线真的很宽.
为什么？仅仅因为计算与内存访问的比率非常低.对于内部循环的每次迭代,您将获取一次列索引为j(8字节),矩阵元素为data(8字节),向量元素的值(8字节)和结果的先前值(因为编译器很少优化)访问共享变量)(8个字节).然后执行2个非常快速的浮点运算(FLOP)并执行存储(虽然+=运算符被转换为单个指令,但它仍然是"fetch-modify-write").总共加载32个字节,然后对它们执行2个FLOP.这使得每字节1/16 FLOP.
具有现代SSE功能的CPU内核可以执行4个双精度FLOP /周期,这通常会导致每个CPU内核8 GFLOPS(假设基本频率为2 GHz).使用AVX,数量增加了一倍,因此在2 GHz Intel Sandy/Ivy Bridge或AMD等效产品上,每个内核最多可获得16 GFLOPS.为了使数据处理能力饱和,给定1/16 FLOP /字节,您需要至少128 GiB/s的内存带宽.
甲高端的Nehalem-EX处理器等至强X7560在2,26千兆赫(9,04 GFLOPS /芯)和它的共享的L3高速缓存运行(L1和L2高速缓存是每核)开约275吉布/秒.在9,04 GFLOPS/core时,每个核心需要144,64 GiB/s才能满足zAx日常工作的内循环.这意味着在理想情况下,此CPU的L3缓存无法提供超过2个完全向量化的乘法内核.
如果没有SSE矢量化,FLOPS速率是双精度的两倍,因此人们可以预期问题可以扩展到4个线程.一旦您的问题变得比L3缓存大,事情变得非常糟糕,因为内存总线的带宽减少了大约十倍.
尝试使用以下版本的内部循环来查看编译器是否足够智能以遵循OpenMP的宽松内存视图:
#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;
}
不幸的是,你无能为力.稀疏矩阵向量乘法与CPU插槽的数量成比例,而不是CPU核心的数量.您需要一个带有独立内存控制器的多插槽系统,例如任何具有多个(后)Nehalem或AMD64处理器的系统.
编辑:优化提示.你真的需要long存储列索引和行指针吗？使用2122980非零元素,您可以使用相当好的方法int.它将在内部循环中保存每个元素加载4个字节,在外部循环中每行加载另外4个字节.