私はキャッシングを改善するために転置することに同意します (しかし、最後にそれについての私のコメントを参照してください)。
参照用の元の関数（私の正気のためにいくつかの整理を加えて）：
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) しかないにもかかわらず) のため、関数自体が頻繁に呼び出されると想定しています。では、なぜメモリを割り当てるのでしょうか。xasをマークしてどこでも__restrict__変更yするだけです (は C99 から取得した GCC 拡張機能です。これには a を使用することをお勧めします。ライブラリには既に 1 つある可能性があります)。x__restrict__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最後に余分な値を含めるように変更できる場合UINT_MAXは、より多くのブランチを削除できます (両方のi<Mチェックで、これも ~N*M 回実行されます)。
これ以上ループを並列にすることはできず、ループを結合することもできません。ブーストは、他の回答で示唆されているように、再配置することになりdます。ただし、再配置dに必要な作業には、ループを実行する作業とまったく同じキャッシュの問題があります。そして、メモリを割り当てる必要があります。良くない。さらに最適化する唯一のオプションは、IncompleteCholeskyFactorization->Diagonalsそれ自体の構造を変更することです。これはおそらく多くの変更を意味します。または、この順序でデータをより適切に処理する別のアルゴリズムを見つけることです。
さらに先に進みたい場合は、最適化がかなり多くのコードに影響を与える必要があります (悪いことではありません。ポインタの配列である正当な理由がない限りDiagonals、リファクタリングでできるようです)。