优化四重嵌套“for”循环

Question

我正在开发一个2D数值模型，并且我想加快特定成员函数的速度，这会减慢我的代码。该函数需要循环遍历模型中的每个网格点，然后在每个网格点上执行双重求和，范围为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]; 
} 

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

Answer 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]; 
    } 
} 

不创建3D阵列double s。它为指向数组double的指针数组创建一个指针数组。每个数组都是它自己的内存块，谁知道它将在哪里存储。这导致了一个称为“poor spacial locality”的数据结构。所有的结构件可能散落在各处。在三维阵列中，您可以跳到三个不同的地方，以查明您的价值在哪里。

由于模拟3D阵列所需的许多存储块可能远不及彼此，因此CPU可能无法提前有效加载高速缓存（高速存储器），必须停止有用的工作它正在等待访问速度较慢的存储，可能更频繁地访问RAM。这是一个很好的，高水平的表现。另一方面，如果整个数组位于一个内存块中，是“连续的”，则CPU可以读取更大的内存块，也许是所有内存块，它需要立即进入缓存。此外，如果编译器知道程序将使用的内存全部在一个大块中，它可以执行各种常规优化，这将使您的程序更快。

那么，我们如何获得一个3D存储器块？如果大小是静态的，这很容易

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

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

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

，做手工的数组索引数学

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

貌似这应该是所有额外的数学比较慢，对吧？原来编译器隐藏的数学看起来很像你每次使用[]时的情况。你几乎没有损失，当然不会像失去一个不必要的cache miss那样多。

但是，在整个地方重复这个数学是疯狂的，迟早你会搞砸了，即使可读性耗尽并不是你首先希望死亡的，所以你真的想把1D数组包装在一个班级帮手处理你的数学。一旦你这样做了，你可能会让这个类处理分配和释放，这样你就可以利用all that RAII goodness。没有更多for环路new s和delete s遍布各处。它全部包裹起来并用弓绑起来。

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

Answer 2

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

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

到（注意没有用+=现在，它在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编译指示可能是一种可以简化并行化循环的方法（ S）。