我希望优化这个简单的循环：
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 -- 这意味着将其存储在可以交换索引的方式中 -- 将 i 作为外部索引：
sub += d[j][i]*x[c];
替代
sub += d[i][j]*x[c];
这应该会导致更好的缓存性能。
我同意为了更好的缓存而进行转置 (但请看我的评论)，还有更多工作要做，所以让我们看看我们可以用完整函数做些什么...
参考原始函数（已经进行了一些整理）：
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;
}
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]。 - Dave
__restrict__？i和j实际上有多大？“这个循环是递归的……”--请提供更多上下文信息。 - Brian Cain