12

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

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（我通常使用它们来提高速度），但这个循环是递归的。

有什么建议么？

4

3 回答 3

10

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

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

代替

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

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

于 2013-08-15T22:11:39.847 回答
6

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

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

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它。也许图书馆已经有一个）。

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

为什么是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包含一个额外的值，则可以删除更多分支（两个检查，再次运行 ~N*M 次）。UINT_MAXi<M

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

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

于 2013-08-16T00:12:40.327 回答
2

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

于 2013-08-16T20:26:00.453 回答