优化四重嵌套“for”循环

Posted 2023-02-16

【中文标题】优化四重嵌套“for”循环

【英文标题】：Optimizing a quadruple nested "for" loop

【发布时间】：2017-03-02 19:37:51

【问题描述】：

我正在用 C++ 开发一个 2D 数值模型，我想加快一个会减慢我的代码速度的特定成员函数。该函数需要遍历模型中的每个i,j 网格点，然后在l 和m 上的每个网格点处执行双重求和。函数如下：

int Class::Function(void) 
    double loadingEta;
    int i,j,l,m;

    //etaLatLen=64, etaLonLen=2*64
    //l_max = 12

    for (i=0; i<etaLatLen; i++) 
        for (j=0; j < etaLonLen; j++) 
            loadingEta = 0.0;
            for (l=0; l<l_max+1; l++) 
                for (m=0; m<=l; m++) 
                    loadingEta += etaLegendreArray[i][l][m] * (SH_C[l][m]*etaCosMLon[j][m] + SH_S[l][m]*etaSinMLon[j][m]);
                
            
            etaNewArray[i][j] = loadingEta;
        
    

    return 1;

我一直在尝试更改循环顺序以加快速度，但无济于事。任何帮助将非常感激。谢谢！

编辑 1：

所有五个数组都在我的类的构造函数中分配如下：

etaLegendreArray = new double**[etaLatLen];
for (int i=0; i<etaLatLen; i++) 
    etaLegendreArray[i] = new double*[l_max+1];
    for (int l=0; l<l_max+1; l++) 
        etaLegendreArray[i][l] = new double[l_max+1];
    


SH_C = new double*[l_max+1];
SH_S = new double*[l_max+1];
for (int i=0; i<l_max+1; i++) 
    SH_C[i] = new double[l_max+1]; 
    SH_S[i] = new double[l_max+1];


etaCosMLon = new double*[etaLonLen];
etaSinMLon = new double*[etaLonLen];
for (int j=0; j<etaLonLen; j++) 
    etaCosMLon[j] = new double[l_max+1];
    etaSinMLon[j] = new double[l_max+1];

也许这些是一维数组而不是多维数组会更好？

【问题讨论】：

更改循环顺序不会降低复杂性。如果您想真正加快速度，您可能希望在多个进程或线程之间分配工作，但这也有开销。

你的数组是如何定义的？您或许可以提高数据的缓存能力。

听起来您正在通过 2D 网格传递 2D 过滤器。所以使用 KissFFT 转换到频域，卷积，然后再转换回空间域。

感谢您的回复。我添加了有关如何将数组分配给帖子的代码。

【参考方案1】：

在这里跳入 X-Y 领域。与其加快算法速度，不如尝试加快数据访问速度。

etaLegendreArray = new double**[etaLatLen];
for (int i=0; i<etaLatLen; i++) 
    etaLegendreArray[i] = new double*[l_max+1];
    for (int l=0; l<l_max+1; l++) 
        etaLegendreArray[i][l] = new double[l_max+1];
    

不创建 doubles 的 3D 数组。它创建一个指向数组的指针数组，该数组指向doubles 的数组。每个数组都是它自己的内存块，谁知道它在存储中的位置。这会产生一个具有所谓“poor spacial locality”的数据结构。结构的所有部分都可能散落在各处。在 3D 数组中，您跳到三个不同的地方只是为了找出您的价值所在。

因为模拟 3D 阵列所需的许多存储块可能彼此相距甚远，CPU 可能无法提前有效地加载缓存（高速内存）并不得不停止它的有用工作做并等待访问较慢的存储，可能更频繁地访问 RAM。这是一个很好的高级article on how much this can hurt 性能。

另一方面，如果整个数组在一个内存块中，并且是“连续的”，CPU 可以读取更大的内存块，也许是全部，它需要一次全部进入缓存。另外，如果编译器知道程序将使用的内存都在一个大块中，它可以执行各种常规优化，从而使您的程序更快。

那么我们如何获得一个全是一个内存块的 3D 数组呢？如果尺寸是静态的，这很容易

double etaLegendreArray[SIZE1][SIZE2][SIZE3];

这看起来不是你的情况，所以你要做的是分配一个一维数组，因为它将是一个连续的内存块。

double * etaLegendreArray= new double [SIZE1*SIZE2*SIZE3];

并手动进行数组索引数学运算

etaLegendreArray[(x * SIZE2 + y) * SIZE3 + z] = data;

看起来应该会因为所有额外的数学运算而变慢，对吧？事实证明，每次使用[] 时，编译器都会隐藏看起来很像的数学。你几乎没有损失任何东西，当然也没有你失去一个不必要的cache miss 那么多。

但是在所有地方重复这个数学是疯狂的，迟早你会搞砸即使可读性的消耗不会让你首先希望死亡，所以你真的想把一维数组包装在一个帮助您处理数学的课程。一旦你这样做了，你还不如让那个类处理分配和释放，这样你就可以利用all that RAII goodness。不再有for 循环news 和deletes 到处都是。全部用蝴蝶结包起来。

Here is an example of a 2D Matrix class easily extendable to 3D. 将以一种很好的可预测和缓存友好的方式处理您可能需要的基本功能。

【参考方案2】：

如果 CPU 支持它并且编译器进行了足够优化，您可能会从 the C99 fma（融合乘加）函数中获得一些小收益，将您的一些两步操作（乘，然后加）转换为一个步操作。它还可以提高准确性，因为融合运算只需要进行一次浮点舍入，而不是乘法和加法一次。

假设我没看错，您可以将最内层循环的表达式更改为：

loadingEta += etaLegendreArray[i][l][m] * (SH_C[l][m]*etaCosMLon[j][m] + SH_S[l][m]*etaSinMLon[j][m]);

to（注意+=现在不使用了，它被合并到fma中）：

loadingEta = fma(etaLegendreArray[i][l][m], fma(SH_C[l][m], etaCosMLon[j][m], SH_S[l][m]*etaSinMLon[j][m]), loadingEta);

我不希望有任何神奇的性能方面的东西，但它可能会有所帮助（同样，只有优化足以让编译器内联硬件指令来完成工作；如果它正在调用库函数，你'将失去对函数调用开销的任何改进）。同样，它应该会通过避免您发生的两个舍入步骤来提高准确性。

请注意，some compilers with appropriate compilation flags, they'll convert your original code to hardware FMA instructions for you;如果这是一个选项，我会选择，因为（如您所见）fma 函数往往会降低代码的可读性。

您的编译器也可能提供浮点指令的矢量化版本，这可能会显着提高性能（请参阅上一个关于自动转换为 FMA 的链接）。

大多数其他改进将需要有关目标、所使用的输入数组的性质等的更多信息。简单的线程可能会为您带来一些好处，OpenMP pragma 可能是一种简化并行化循环的方法（s )。

【讨论】：

可能还值得一提的是在对浮点数求和时所涉及的陷阱（即：排序，所以首先对最小值求和）。根据架构（即：嵌入式），如果整数运算明显更快，使用定点归一化/求和也可能是值得的。

有趣，我不知道 FMA。我在 g++6 中编译，所以我会看看 FMA 是否是内置优化的一部分。谢谢！