私はこの単純なループ最適化したいと思い、この短いループ

Question

最適化したい：

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のブルドーザーで約1秒の実行時間を持っています。私はSIMDとOpenMPについて考えましたが、これは通常より速くするために使用しますが、このループは再帰的です。

提案がありますか？

Answer 1

...

オリジナル機能、参照のための（私の正気のためのいくつかの以降の整理整頓を持つ）：

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を使用してもかまいません。

同様に、署名を変更することはできないと思いますが、関数には1つのパラメータしか取らずに変更することができます。 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の値を追加するように変更できる場合は、より多くのブランチ（〜N * M回実行するi<Mの両方のチェック）を削除することができます。

これ以上ループを並列化することはできません。ループを組み合わせることはできません。ブーストは、他の答えに示唆されているように、dを並べ替えることになります。例外を除いて... dを再配置するために必要な作業は、ループを実行する作業とまったく同じキャッシュ問題を持ちます。そして、メモリを割り当てる必要があります。良くない。さらに最適化するオプションは、IncompleteCholeskyFactorization->Diagonals自体の構造を変更することです。これはおそらく多くの変更を意味するか、データをこの順序でうまく処理する別のアルゴリズムを見つけることになります。

さらに進めたい場合は、コードのかなりの部分に最適化が必要です（悪いことではありません; Diagonalsがポインタの配列である理由がない限り、リファクタリング）。

Answer 2

を使用すると、行列Dを移調したいかもしれないと思う - の手段は、あなたがインデックスを交換することができ、このような方法でそれを保存する - 私アウターインデックスます

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

代わりの

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

これにより、キャッシュのパフォーマンスが向上します。それでは、私たちはフル機能で何ができるか見てみましょう、私は（しかし、最後にその上で私のコメントを参照してください）、より良いキャッシュのために転置に同意し、やるべきことがあります

Answer 3

私は自分の質問に答えたいと思います：悪いパフォーマンスは、（少なくとも）Win7が大きなメモリブロックを同じ境界に揃えているために、キャッシュ競合ミスが原因です。私のケースでは、すべてのバッファについて、アドレスは同じアラインメント（bufferadress％4096はすべてのバッファで同じでした）なので、L1キャッシュの同じキャッシュセットに入ります。キャッシュの競合ミスを避けるためにバッファを別の境界に合わせるためにメモリ割り当てを変更し、第2因子の高速化を得ました。すべての答え、特にDaveの回答に感謝します！