我尝试应用循环展开来通过Gauss Jorda方法找到矩阵的逆矩阵,当矩阵的大小非常大,无法放入缓存时,减少内存访问的次数(瓶颈)。我让它运行得更快,但我得到的结果是错误的,我不知道为什么。
for(k=0;k<size;k++)
{
pivot=original[k][k];
for(j=0;j<size;j++)
{
original[k][j]/=pivot;
inverse[k][j]/=pivot;
}
for(i=0;i<size;i++)
{
if(i!=k)
{
pivot = original[i][k];
for(j=0;j<size;j++)
{
original[i][j] -= original[k][j]*pivot;
inverse[i][j] -= inverse[k][j]*pivot;
}
}
}
}
我希望通过接收内存访问次数来加快问题的执行速度。
我的循环展开版本如下:
for(k=0;k<sizeOfMatrix;k += 2)
{
pivot=original[k][k];
for(j=0;j<sizeOfMatrix;j++)
{
original[k][j]/=pivot;
inverse[k][j]/=pivot;
}
for(i=0;i<sizeOfMatrix;i++)
{
if(i!=k && i!=k+1)
{
pivot = original[i][k];
for(j=0;j<sizeOfMatrix;j++)
{
original[i][j] -= original[k][j]*pivot;
inverse[i][j] -= inverse[k][j]*pivot;
}
}
if(i!=k+1)
{
pivot=original[k][k];
for(j=0;j<sizeOfMatrix;j++)
{
original[k+1][j]/=pivot;
inverse[k+1][j]/=pivot;
}
pivot = original[i][k+1];
for(j=0;j<sizeOfMatrix;j++)
{
original[i][j] -= original[k+1][j]*pivot;
inverse[i][j] -= inverse[k+1][j]*pivot;
}
}
}
}
2条答案
按热度按时间watbbzwu1#
现在,让我们先不考虑原始代码的数值适应性问题,以及展开是否是您应该尝试的练习。您展示了通过重写 middle 循环体来展开 outer 循环的尝试。尽管可以沿着与转换结果类似的方式构建与原始代码等效的计算。它不会构成一个展开,因为循环展开不会改变在循环迭代范围内循环体中执行的操作的顺序。2你所做的转换确实会产生这样的改变。
但在您的例子中,转换后的代码就是错误的,至少因为
sizeOfMatrix
为奇数时,它在k
循环迭代期间执行越界数组访问,其中k == sizeOfMatrix - 1
。if(i!=k+1)
块的开始处使用了错误的主元值。在不进一步考虑您所尝试的转换是否有助于性能改进的情况下,我观察到,如果您想执行真正的展开,那么尝试展开的循环是错误的。如果真正的展开有任何好处,那么它将来自于展开最内层的循环。例如:
但是正如我在注解中所说的,如果你的编译器期望这样的展开来加速程序,那么当你在编译时启用优化时,它会自动这样做,而且编译器对这些事情的判断是相当不错的。
实际上,手动优化很容易使性能恶化,而不是提高性能,这可能是固有的,也可能是因为手动优化会干扰编译器自己执行更有效的优化。
如果最佳速度是你的目标,那么你可能应该使用一个调优的第三方实现,比如lapack或者atlas。
如果您坚持自己滚动,那么您可能只需改进算法(甚至坚持使用G-J)就可以获得更高的加速比,但您可能应该给予其中一些返回以获得更高的数值稳定性。
1.在每次外循环迭代时,选择当前消除列中具有最大量值的可用行,并在必要时将其与当前行交换。这将最小化大数除以非常小的数(尤其是零)的发生率。
1.(在上面的(1)之后),仅计算当前透视表右侧列中元素的更新。原始矩阵的其他元素已经对结果产生了所有影响。(但您不能使用逆矩阵来简化这种方法。)
如果你不对(1)做一些改变,那么你的逆矩阵将在一些实际上可逆的矩阵上失败。
沿着如此,例如:
这是从您的第一个代码中得到的灵感,改进了名称,并实现了算法改进(1)和(2)。
pbpqsu0x2#
你已经找到了问题的解决方案。因为目标是学术性的。目标是减少执行时间,减少内存访问的次数,在我的情况下减少一半。
1.我修改了第k行,我用它的主元划分了该行的元素
1.用新的行k修改矩阵的行k+1
1.修改行k+1,将行的元素除以其透视。
1.我修改了值为k和k + 1的其他行