我想优化这个短循环答案

【问题标题】：I want to optimize this short loop我想优化这个短循环
【发布时间】：2013-08-15 21:55:11
【问题描述】：

我想优化这个简单的循环：

unsigned int i;
while(j-- != 0){ //j is an unsigned int with a start value of about N = 36.000.000
   float sub = 0;
   i=1;
   unsigned int c = j+s[1];
   while(c < N) {
       sub += d[i][j]*x[c];//d[][] and x[] are arrays of float
       i++;
       c = j+s[i];// s[] is an array of unsigned int with 6 entries.
   }
   x[j] -= sub;                        // only one memory-write per j
}

对于 4000 MHz AMD Bulldozer，循环的执行时间约为一秒。我考虑过 SIMD 和 OpenMP（我通常使用它们来提高速度），但这个循环是递归的。

有什么建议吗？

【问题讨论】：

你能包含更多上下文吗？ “d[][] 和 x[] 是浮点数组”——它们是否声明为 __restrict__？ i 和 j 究竟有多大？ “这个循环是递归的……”——请提供更多上下文。
优化后的编译器输出是什么样的？
另外，从分析器开始查看热点在哪里。
这不是递归循环（没有这样的东西。）这是一个嵌套在另一个循环中的循环。
您可以进行各种调整，但要真正提升，请尝试在更大的级别上优化代码（您没有发布足够的周边代码让我提供任何建议）

标签： c++ performance algorithm

【解决方案1】：

认为您可能想要转置矩阵 d——意味着以可以交换索引的方式存储它——使 i 成为外部索引：

    sub += d[j][i]*x[c];

而不是

    sub += d[i][j]*x[c];

这应该会带来更好的缓存性能。

【讨论】：

很好，我认为这实际上是导致每个循环花费大量时间的原因
我尝试切换循环，但随后我获得了更多的内存写入，并且速度变慢了。我还尝试了使用切换循环的 OpenMP，但由于缓存冲突，它变得更慢。矩阵的转置是可能的，但不适合该算法的其他部分。
好的！也许我应该给你更多的信息。该循环是此文件中 CholeskyBackSolve 的一部分code.google.com/p/rawtherapee/source/browse/rtengine/… 我在上面发布的版本已经进行了一些优化，并且是从 Ingo 335 行开始的 while 的替代品
转置，使用良好的缓存局部性进行计算，然后再次转置可能仍然更快。
这是虚拟内存计算系统中最古老的问题。总是问自己你的数组是按行优先还是列优先存储的，并处理它们，以便访问相邻的元素而不是相距很远的元素。

【解决方案2】：

我同意转置以获得更好的缓存（但最后请参阅我的 cmets），还有更多工作要做，所以让我们看看我们可以用完整的功能做什么......

原始函数，供参考（为了我的理智而进行了一些整理）：

void MultiDiagonalSymmetricMatrix::CholeskyBackSolve(float *x, float *b){
    //We want to solve L D Lt x = b where D is a diagonal matrix described by Diagonals[0] and L is a unit lower triagular matrix described by the rest of the diagonals.
    //Let D Lt x = y. Then, first solve L y = b.

    float *y = new float[n];
    float **d = IncompleteCholeskyFactorization->Diagonals;
    unsigned int *s = IncompleteCholeskyFactorization->StartRows;
    unsigned int M = IncompleteCholeskyFactorization->m;
    unsigned int N = IncompleteCholeskyFactorization->n;
    unsigned int i, j;
    for(j = 0; j != N; j++){
        float sub = 0;
        for(i = 1; i != M; i++){
            int c = (int)j - (int)s[i];
            if(c < 0) break;
            if(c==j) {
                sub += d[i][c]*b[c];
            } else {
                sub += d[i][c]*y[c];
            }
        }
        y[j] = b[j] - sub;
    }

    //Now, solve x from D Lt x = y -> Lt x = D^-1 y
    // Took this one out of the while, so it can be parallelized now, which speeds up, because division is expensive
#pragma omp parallel for
    for(j = 0; j < N; j++){
        x[j] = y[j]/d[0][j];
    }

    while(j-- != 0){
        float sub = 0;
        for(i = 1; i != M; i++){
            if(j + s[i] >= N) break;
            sub += d[i][j]*x[j + s[i]];
        }
        x[j] -= sub;
    }
    delete[] y;
}

由于有关并行除法的评论可以提高速度（尽管只有 O(N)），我假设函数本身被调用了很多次。那么为什么要分配内存呢？只需将x 标记为__restrict__ 并将y 更改为x（__restrict__ 是一个 GCC 扩展，取自 C99。您可能想使用 define 来代替它。也许图书馆已经有一个)。

同样，虽然我猜你不能改变签名，但你可以让函数只接受一个参数并修改它。当 x 或 y 已设置时，永远不会使用 b。这也意味着您可以摆脱运行 ~N*M 次的第一个循环中的分支。如果必须有 2 个参数，请在开头使用 memcpy。

为什么d 是一个指针数组？一定是吗？这在原始代码中似乎太深了，所以我不会触及它，但是如果有任何展平存储数组的可能性，即使您无法转置它也会提高速度（乘法，加法，取消引用更快）而不是取消引用、添加、取消引用）。

所以，新代码：

void MultiDiagonalSymmetricMatrix::CholeskyBackSolve(float *__restrict__ x){
    // comments removed so that suggestions are more visible. Don't remove them in the real code!
    // these definitions got long. Feel free to remove const; it does nothing for the optimiser
    const float *const __restrict__ *const __restrict__ d = IncompleteCholeskyFactorization->Diagonals;
    const unsigned int *const __restrict__ s = IncompleteCholeskyFactorization->StartRows;
    const unsigned int M = IncompleteCholeskyFactorization->m;
    const unsigned int N = IncompleteCholeskyFactorization->n;
    unsigned int i;
    unsigned int j;
    for(j = 0; j < N; j++){ // don't use != as an optimisation; compilers can do more with <
        float sub = 0;
        for(i = 1; i < M && j >= s[i]; i++){
            const unsigned int c = j - s[i];
            sub += d[i][c]*x[c];
        }
        x[j] -= sub;
    }

    // Consider using processor-specific optimisations for this
#pragma omp parallel for
    for(j = 0; j < N; j++){
        x[j] /= d[0][j];
    }

    for( j = N; (j --) > 0; ){ // changed for clarity
        float sub = 0;
        for(i = 1; i < M && j + s[i] < N; i++){
            sub += d[i][j]*x[j + s[i]];
        }
        x[j] -= sub;
    }
}

嗯，它看起来更整洁，缺少内存分配和减少分支，如果没有别的，是一种提升。如果您可以更改 s 以在末尾包含一个额外的 UINT_MAX 值，则可以删除更多分支（i<M 检查，再次运行 ~N*M 次）。

现在我们不能让更多的循环并行，也不能合并循环。正如另一个答案中所建议的那样，现在的提升将是重新排列d。除了……重新排列d 所需的工作与执行循环的工作具有完全相同的缓存问题。它需要分配内存。不好。唯一需要进一步优化的选项是：更改 IncompleteCholeskyFactorization->Diagonals 本身的结构，这可能意味着很多更改，或者找到一种更适合按此顺序处理数据的不同算法。

如果您想走得更远，您的优化将需要影响相当多的代码（这不是一件坏事；除非有充分的理由让 Diagonals 成为指针数组，否则它似乎可以与重构）。

【讨论】：

如果你做UINT_MAX的事情，记得重新排列第二个等式以避免溢出：N - j > s[i]
非常感谢您的详细回答！我按照您的建议消除了 y 并限制了 x。这给了一个小的加速。我通过完全不同的事情获得的最大加速。我更改了对角线的内存分配方式，以便每个对角线从不同的 64 字节边界开始。这使速度提高了 2 倍 :-) Bulldozer 的只有 4 路关联 L1 缓存真的是一团糟……现在我会更进一步，也许我会完全重新设计该代码（这不是从我这里，我只尝试优化它）

【解决方案3】：

我想回答我自己的问题：性能不佳是由于（至少）Win7 将大内存块对齐到同一边界而导致缓存冲突未命中造成的。就我而言，对于所有缓冲区，地址具有相同的对齐方式（所有缓冲区的缓冲区地址 % 4096 相同），因此它们属于 L1 缓存的相同缓存集。我更改了内存分配以将缓冲区对齐到不同的边界以避免缓存冲突未命中，并获得了 2 倍的加速。感谢所有答案，尤其是 Dave 的答案！

【讨论】：