在尝试并行化之前,先尝试简化。例如,可以完全消除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
:
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
循环中是不变的,因此我们可以预先保存它:
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;
}
}
矩阵的访问方式对缓存性能来说可能是最差的。因此,如果可以调整维度的分配顺序,就可以实现大量节省内存访问的目标:
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循环中,它是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
的维度改回原始方式但保留其他矩阵交换版本的新版本。这提供了所有矩阵的最佳缓存性能。也就是说,大多数矩阵访问的内部循环是沿着缓存行递增的。
实际上,请试试。它可能会将事情改进到不需要并行化的程度。我怀疑代码已经受限于内存,因此并行可能并没有那么有帮助。
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
比率越高越好。
无论如何,这是我能做到的极限。祝你好运...
U
,L
和Aprime
是全局范围还是函数范围?这会影响到openmp
。有多线程原子更新的考虑,但openmp
有处理这个问题的工具。所以,你能贴出更多的代码吗?此外,你有一些不变量,所以看起来可以通过一些重构使其更可并行化。通过拆分事物,你可以从j
循环中移除/移动if
。另外,n
的值是多少?在我看来,它是可以并行化的。 - Craig Esteytime
足够了),特别是涉及内存访问优化的“之前”和“之后”,那就太好了。我发现我的一个旧版本不能充分利用简单(仅有约10%的)并行化,而且我怀疑错误的内存访问是罪魁祸首。 - deamentiaemundi