ループタイル/ブロックによる高密度行列の乗算

Question

大規模な高密度行列の乗算に効果的にループタイリング/ループブロッキングを使用する方法を教えてもらえないかと疑問に思っていました。私はやっているC = AB 1000x1000行列です。私はループタイリングのためのWikipediaの例に従ってきましたが、タイル張りを使用するよりも悪い結果が得られます。

http://en.wikipedia.org/wiki/Loop_tiling

http://software.intel.com/en-us/articles/how-to-use-loop-blocking-to-optimize-memory-use-on-32-bit-intel-architecture

私は以下のいくつかのコードを提供しています。 naiveメソッドはキャッシュミスのために非常に遅いです。転置方法は、Bの転置をバッファーに作成します。この方法では、行列乗算はO（n^3）となり、転置はO（n^2）となるため、転置は少なくとも1000倍高速です。ブロッキングのないwikiメソッドも高速で、バッファーは必要ありません。ブロック方法は遅いです。ブロッキングに関するもう1つの問題は、ブロックを数回更新する必要があることです。これはスレッド/ OpenMPにとって挑戦です。なぜなら、それが慎重でなければ競合状態を引き起こす可能性があるからです。

私は、転置方法の変更でAVXを使用すると、結果がEigenよりも速くなることを指摘しておきます。しかし、SSEの私の結果はEigenより少し遅いので、私はより良いキャッシュを使うことができると思います。

編集： 私は自分がしたいことを考えていると思います。これは、1969年
私はブロック行列の乗算を行い、各スレッドがCではなく、とBのサブ行列を扱う持っている必要があり http://en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms

からキャノンのアルゴリズムから来ています。例えば、私の行列を4つのブロックに分割するとします。私はそうするでしょう：

[C1 C2  [A1 A2 [B1 B2 
C3 C4] = A3 A4] x B3 B4] 
thread 1: C1 = A1B1 + A2B3 
thread 2: C2 = A1B2 + A2B4 
thread 3: C3 = A3B1 + A4B3 
thread 4: C4 = A3B2 + A4B4 

それは競合状態を取り除きます。私はこれについて考える必要があります。

void matrix_mult_naive(const float*A , const float* B, float* C, const int N, const int M, const int K) { 
    #pragma omp parallel for 
    for(int i=0; i<N; i++) { 
     for(int j=0; j<K; j++) { 
      float tmp = 0; 
      for(int l=0; l<M; l++) { 
       tmp += A[M*i+l]*B[K*l+j]; 
      } 
      C[K*i + j] = tmp; 
     } 
    } 
} 
void matrix_mult_transpose(const float*A , const float* B, float* C, const int N, const int M, const int K) { 
    float *B2 = (float*)aligned_malloc(M*K*sizeof(float), 32); 
    transpose(B, B2, M, K, 1); 
    #pragma omp parallel for 
    for(int i=0; i<N; i++) { 
     for(int j=0; j<K; j++) { 
      float tmp = 0; 
      for(int l=0; l<M; l++) { 
       tmp += A[M*i+l]*B2[M*j+l]; 
      } 
      C[K*i + j] = tmp; 
     } 
    } 
    aligned_free(B2); 
} 

void matrix_mult_wiki(const float*A , const float* B, float* C, const int N, const int M, const int K) { 
    for(int i=0; i<N; i++) { 
     for(int j=0; j<K; j++) { 
      C[K*i + j] = 0; 
     } 
    } 
    #pragma omp parallel for 
    for(int i=0; i<N; i++) { 
     for(int l=0; l<M; l++) { 
      float a = A[M*i+l]; 
      for(int j=0; j<K; j++) { 
       C[K*i + j] += a*B[K*l+j]; 
      } 
     } 
    } 
} 

void matrix_mult_wiki_block(const float*A , const float* B, float* C, const int N, const int M, const int K) { 
    const int block_size = 8; //I have tried several different block sizes 
    for(int i=0; i<N; i++) { 
     for(int j=0; j<K; j++) { 
      C[K*i + j] = 0; 
     } 
    } 
    for(int l2=0; l2<M; l2+=block_size) { 
     for(int j2=0; j2<K; j2+=block_size) { 
     #pragma omp parallel for 
      for(int i=0; i<N; i++) { 
       for(int l=l2; l<min(M, l2+block_size); l++) { 
        for(int j=j2; j<min(K, j2+block_size); j++) { 
         C[K*i + j] += A[M*i+l]*B[K*l+j]; 
        } 
       } 
      } 
     } 
    } 
} 

Answer 1

私が得ている最良の結果は、あなたのN以上のブロック、およびループを並べ替えることによって1つの以上forループを追加することです。私もループ不変のコードを持ち上げましたが、コンパイラのオプティマイザはうまくいけばこれを自動的に行うはずです。ブロックサイズは、キャッシュラインサイズをsizeof(float)で割ったものでなければなりません。これは、転置されたアプローチより約50％速くなっています。

AVXまたはブロッキングのいずれかを選択する必要がある場合は、AVX拡張子（vfmadd###psおよびhaddps）を使用するとまだかなり高速です。ブロックサイズが64/sizeof(float) == 16 floats == 2つの256ビットAVXレジスタの倍数である場合には、すでにテストしているので、両方を使うのが最も簡単です。1816522は、

タイルをマダニ：892431（51％速い）

AVX +タイリング：230512（87％速い）

タイル：

void matrix_mult_wiki_block(const float*A , const float* B, float* C, 
          const int N, const int M, const int K) { 
    const int block_size = 64/sizeof(float); // 64 = common cache line size 
    for(int i=0; i<N; i++) { 
     for(int j=0; j<K; j++) { 
      C[K*i + j] = 0; 
     } 
    } 
    for (int i0 = 0; i0 < N; i0 += block_size) { 
     int imax = i0 + block_size > N ? N : i0 + block_size; 

     for (int j0 = 0; j0 < M; j0 += block_size) { 
      int jmax = j0 + block_size > M ? M : j0 + block_size; 

      for (int k0 = 0; k0 < K; k0 += block_size) { 
       int kmax = k0 + block_size > K ? K : k0 + block_size; 

       for (int j1 = j0; j1 < jmax; ++j1) { 
        int sj = M * j1; 

        for (int i1 = i0; i1 < imax; ++i1) { 
         int mi = M * i1; 
         int ki = K * i1; 
         int kij = ki + j1; 

         for (int k1 = k0; k1 < kmax; ++k1) { 
          C[kij] += A[mi + k1] * B[sj + k1]; 
         } 
        } 
       } 
      } 
     } 
    } 
} 

---

転置

• キャノンの参考文献については、SUMMA algorithm従う方が良いです。

  ---

  誰が（私はここに終わった方法{〜1E9×50} X {50×50}、）トールスキニー乗算を最適化する場合には、転置アプローチがブロックアプローチに性能がほぼ同一でありますn = 18（浮動小数点数）までです。 n = 18は病理学的なケースです（17または19よりも悪い）、これを引き起こすキャッシュアクセスパターンはあまりありません。ブロックされたアプローチでは、全てが大きくなります。n