use*_*675 4 c c++ parallel-processing openmp
我得到了一些使用OpenMP进行并行化的代码,并且在各种函数调用中,我注意到这个for循环对计算时间有一些好的负罪感.
double U[n][n];
double L[n][n];
double Aprime[n][n];
for(i=0; i<n; i++) {
for(j=0; j<n; j++) {
if (j <= i) {
double s;
s=0;
for(k=0; k<j; k++) {
s += L[j][k] * U[k][i];
}
U[j][i] = Aprime[j][i] - s;
} else if (j >= i) {
double s;
s=0;
for(k=0; k<i; k++) {
s += L[j][k] * U[k][i];
}
L[j][i] = (Aprime[j][i] - s) / U[i][i];
}
}
然而,在尝试并行化并在这里和那里应用一些信号量之后(没有运气),我开始认识到else if条件对早期的强烈依赖if(L[j][i]是一个处理过的数字U[i][i],可能在早期设置if)在我的意见中,由于竞争条件,使其不可并行化.
是否可以以这种方式并行化此代码,以便else if只有在早期if已经完成时才执行？
在尝试并行化之前,先尝试简化.
例如,if可以完全消除.
此外,代码以导致最差缓存性能的方式访问矩阵.这可能是真正的瓶颈.
注意:在下面的更新#3中,我做了基准测试fix5,从更新#2开始,缓存友好版本的性能优于原始版本3.9倍.
我已经分阶段清理了一些东西,所以你可以看到代码转换.
有了这个,应该可以omp成功添加指令.正如我在最高评论中提到的,变量的全局与函数范围会影响可能需要的更新类型(例如omp atomic update,等等)
作为参考,这是您的原始代码:
double U[n][n];
double L[n][n];
double Aprime[n][n];
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
if (j <= i) {
double s;
s = 0;
for (k = 0; k < j; k++) {
s += L[j][k] * U[k][i];
}
U[j][i] = Aprime[j][i] - s;
}
else if (j >= i) {
double s;
s = 0;
for (k = 0; k < i; k++) {
s += L[j][k] * U[k][i];
}
L[j][i] = (Aprime[j][i] - s) / U[i][i];
}
}
}
这else if (j >= i)是不必要的,可以直接替换else.但是,我们可以将j循环分成两个循环,这样就不需要if/else:
// fix2.c -- split up j's loop to eliminate if/else inside
double U[n][n];
double L[n][n];
double Aprime[n][n];
for (i = 0; i < n; i++) {
for (j = 0; j <= i; j++) {
double s = 0;
for (k = 0; k < j; k++)
s += L[j][k] * U[k][i];
U[j][i] = Aprime[j][i] - s;
}
for (; j < n; j++) {
double s = 0;
for (k = 0; k < i; k++)
s += L[j][k] * U[k][i];
L[j][i] = (Aprime[j][i] - s) / U[i][i];
}
}
U[i][i]在第二个j循环中是不变的,所以我们可以预先设定它:
// fix3.c -- save off value of U[i][i]
double U[n][n];
double L[n][n];
double Aprime[n][n];
for (i = 0; i < n; i++) {
for (j = 0; j <= i; j++) {
double s = 0;
for (k = 0; k < j; k++)
s += L[j][k] * U[k][i];
U[j][i] = Aprime[j][i] - s;
}
double Uii = U[i][i];
for (; j < n; j++) {
double s = 0;
for (k = 0; k < i; k++)
s += L[j][k] * U[k][i];
L[j][i] = (Aprime[j][i] - s) / Uii;
}
}
对矩阵的访问可能是缓存性能最差的方式.因此,如果可以翻转维度的分配,则可以实现内存访问的大量节省:
// fix4.c -- transpose matrix coordinates to get _much_ better memory/cache
// performance
double U[n][n];
double L[n][n];
double Aprime[n][n];
for (i = 0; i < n; i++) {
for (j = 0; j <= i; j++) {
double s = 0;
for (k = 0; k < j; k++)
s += L[k][j] * U[i][k];
U[i][j] = Aprime[i][j] - s;
}
double Uii = U[i][i];
for (; j < n; j++) {
double s = 0;
for (k = 0; k < i; k++)
s += L[k][j] * U[i][k];
L[i][j] = (Aprime[i][j] - s) / Uii;
}
}
更新:
在Op的第一个k-loop它
k<j和k<i你的第二个不是你必须修复它？
是的,我已修好了.这是一个太难看的变化fix1.c,所以我删除了它并将更改应用到fix2-fix4容易做的地方.
更新#2:
这些变量都是函数的本地变量.
如果你的意思是它们是函数作用域[ 没有 static ],这表示矩阵不能太大,因为除非代码增加了堆栈大小,否则它们仅限于堆栈大小限制(例如8 MB)
虽然矩阵似乎是VLA [因为n是小写],但我忽略了这一点.您可能想尝试使用固定维度数组的测试用例,因为我相信它们可能更快.
此外,如果矩阵是函数范围,并且想要并行化事物,那么您可能需要做(例如)#pragma omp shared(Aprime) shared(U) shared(L).
缓存的最大阻力是要计算的循环s.在fix4,我能够访问U缓存友好,但L访问很差.
如果我确实包含外部上下文,我需要发布更多内容
我猜到了,所以我推测性地进行了矩阵维度交换,不知道其他代码需要改变多少.
我创建了一个新版本,它将尺寸改L回原来的方式,但保持交换版本与其他版本保持一致.这为所有矩阵提供了最佳的缓存性能.也就是说,用于大多数矩阵访问的内循环使得每次迭代沿着高速缓存行递增.
事实上,试一试.它可以将事情改进到不需要并行的程度.我怀疑代码是内存绑定的,所以并行可能没那么多帮助.
// fix5.c -- further transpose to fix poor performance on s calc loops
//
// flip the U dimensions back to original
double U[n][n];
double L[n][n];
double Aprime[n][n];
double *Up;
double *Lp;
double *Ap;
for (i = 0; i < n; i++) {
Ap = Aprime[i];
Up = U[i];
for (j = 0; j <= i; j++) {
double s = 0;
Lp = L[j];
for (k = 0; k < j; k++)
s += Lp[k] * Up[k];
Up[j] = Ap[j] - s;
}
double Uii = Up[i];
for (; j < n; j++) {
double s = 0;
Lp = L[j];
for (k = 0; k < i; k++)
s += Lp[k] * Up[k];
Lp[i] = (Ap[j] - s) / Uii;
}
}
即使你真的需要原始尺寸,根据其他代码,你也可以进行转置并转换回去.这会使其他代码的内容保持不变,但是,如果此代码确实是一个瓶颈,额外的转置操作可能会小到足以值得这样做.
更新#3:
我在所有版本上都运行了基准测试.以下是相对于原始的经过时间和比率n等于1037:
orig: 1.780916929 1.000x
fix1: 3.730602026 0.477x
fix2: 1.743769884 1.021x
fix3: 1.765769482 1.009x
fix4: 1.762100697 1.011x
fix5: 0.452481270 3.936x
更高的比率更好.
无论如何,这是我能做的极限.祝你好运......