多次元行列加算を高速化するには？

Question

私は、サイズ（m * l * 4）の行列Aを持ち、mのサイズは約100,000、l = 100です。リストのサイズは常にnと等しく、n < = mです。私は与えられたインデックスのリストを行列に加えたいと思っていました。 私はこの関数を書いており、この関数を何度も何度も呼び出す必要があります。私は、全体のコード内の関数の各部分を取り、私はMatrixAddition機能で撮影した時の私の60％を発見したどのくらいの時間を計算するためには、gprofを使用

void MatrixAddition(int l, int n, vector<int>& list, int ***A,int ***C,int cluster) 
{ 
    for(int i=0;i<l;i++) 
    { 
     for(int j=0;j<4;j++) 
       C[cluster][i][j]=0; 
    } 

for (int i = 0; i < l; i++) 
{ 
     for(int j=0;j<n;++j) 
    { 
     for(int k=0;k<4;k++) 
      C[cluster][i][k]+=A[list[j]][i][k]; 
    } 
} 

} 

。私の実行時間を減らすために、この関数を書く別の方法がありますか？

時間seconds秒
52.00 7.85 7.85 20 392.60 405.49 MatrixAddition（int型、int型、のstd ::ベクトル> &、int型***、int型***、int型）

名を呼ぶ/ MS/MSを呼び出すを呼び出します

Answer 1

iでループし、2番目の部分でjでループします。これにより、関数はよりキャッシュに適しています。

for(int j=0;j<n;++j) 
{ 
    for (int i = 0; i < l; i++) 
    { 
     for(int k=0;k<4;k++) 
      C[cluster][i][k]+=A[list[j]][i][k]; 
    } 
} 

また、-O3フラグを忘れないでください。

Answer 2

（更新：以前のバージョンではインデックスが間違っていました。このバージョンはかなり簡単に自動ベクトル化されます）。

ポインタ（ポインタへのポインタの配列ではなく）のC多次元配列、またはi*cols + jでインデックス付けされたフラット1D配列を使用して、メモリが連続しているようにします。これは、メモリ帯域幅を有効に活用するために、ハードウェアプリフェッチの有効性に大きな違いをもたらします。負荷が実際に別の負荷から来ている場合は、パフォーマンスが低下し、逆に予測可能なアドレスがあらかじめわかっていれば、負荷が必要になる前に開始することができます。

また、@ user31264の答えは正しいです。jのループが最も外側になるように、ループを交換する必要があります。これは良いですが、それだけでは十分ではありません。

これにより、コンパイラは自動的にベクトル化することもできます。実際、gccに自動ベクトル化をさせるのは驚くほど苦労しました。 （私はインデックスが間違っていましたので、私はコードのみで初めて見えたので、しかし、それは、おそらくです。コンパイラは、我々は連続したメモリをループし、知りませんでしたので。）

---

は私がそれで遊んGodbolt compiler explorer。

私は最終的にAとCのようなフラット1Dアレイを取り、インデックス自体を行い、このバージョンから良い素敵なコンパイラの出力を、得た：

void MatrixAddition_contiguous(int rows, int n, const vector<int>& list, 
           const int *__restrict__ A, int *__restrict__ C, int cluster) 
    // still auto-vectorizes without __restrict__, but checks for overlap and 
    // runs a scalar loop in that case 
{ 
    const int cols = 4; // or global constexpr or something 
    int *__restrict__ Ccluster = C + ((long)cluster) * rows * cols; 

    for(int i=0;i<rows;i++) 
    //#pragma omp simd 
    for(int k=0;k<4;k++) 
     Ccluster[cols*i + k]=0; 

    for(int j=0;j < cols;++j) { // loop over clusters in A in the outer-most loop 
    const int *__restrict__ Alistj = A + ((long)list[j]) * rows * cols; 
    // #pragma omp simd // Doesn't work: only auto-vectorizes with -O3 
    // probably only -O3 lets gcc see through the k=0..3 loop and treat it like one big loop 
    for (int i = 0; i < rows; i++) { 
     long row_offset = cols*i; 
     //#pragma omp simd // forces vectorization with 16B vectors, so it hurts AVX2 
     for(int k=0;k<4;k++) 
     Ccluster[row_offset + k] += Alistj[row_offset + k]; 
    } 
    } 
} 

手動list[j]は間違いなくCに保存することを実現し、コンパイラを助けた巻き上げ読み込まれるインデックスには影響しませんlist[j]。おそらく他のものを手動で持ち上げる必要はなかったでしょう。

ホイストA[list[j]]ではなく、list[j]は、previous approach where I had the indexing wrongのアーチファクトです。可能であれば、list[j]から負荷を持ち上げる限り、listがCと重複していないことがわからなくても、コンパイラはうまくいく可能性があります。GCC 5.3は、x86-64の-O3 -Wall -march=haswell -fopenmp（および-fverbose-asm）を標的と

内側のループは、次のとおりです。

.L26: 
    vmovdqu ymm0, YMMWORD PTR [r8+rax]  # MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B] 
    vpaddd ymm0, ymm0, YMMWORD PTR [rdx+rax] # vect__71.75, MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], MEM[base: vectp.73_90, index: ivtmp.91_26, offset: 0B] 
    add  r12d, 1 # ivtmp.88, 
    vmovdqu YMMWORD PTR [r8+rax], ymm0  # MEM[base: Ccluster_20, index: ivtmp.91_26, offset: 0B], vect__71.75 
    add  rax, 32 # ivtmp.91, 
    cmp  r12d, r9d # ivtmp.88, bnd.66 
    jb  .L26  #, 

だから、8をやっている非整列ロードおよび背面に整列していない店で、AVX2 vpadddで、一度に追加C.

自動ベクトル化されているので、ARM NEON、PPC Altivec、またはパック32ビット追加を行うことができるもので良いコードにする必要があります。

-ftree-vectorizer-verbose=2で何かを教えてもらえませんでしたが、clangの-Rpass-analysis=loop-vectorizeは少し役に立ちました。