c++ - 优化四重嵌套 "for"循环

Question

我正在用 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];
}

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

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

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

由于模拟 3D 阵列所需的许多存储 block 可能彼此相距甚远，因此 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 .

但是到处重复这个数学是很疯狂的，迟早你会搞砸的，即使可读性的下降并没有让你先想死，所以你真的想把一维数组包装在一个帮助你处理数学的类(class)。一旦你这样做了，你还不如让那个类处理分配和释放，这样你就可以利用 all that RAII goodness .不再到处都是 new 和 delete 的 for 循环。全部包裹起来并系上蝴蝶结。

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