układanie pętli/blokowanie dużych mnożenia zwartych matryc

Question

Zastanawiam się, czy ktoś może mi pokazać, jak skutecznie wykorzystać pętlowanie/blokowanie pętli w celu efektywnego zwielokrotnienia dużej gęstości macierzy. Wykonuję C = AB z macierzami 1000x1000. Podążałem za przykładem na Wikipedii w kwestii układania pętli, ale mam gorsze wyniki przy użyciu kafelkowania niż bez.

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

mam pod warunkiem jakiś poniższy kod. Metoda naiwna jest bardzo powolna z powodu chybienia pamięci podręcznej. Metoda transpozycji powoduje przetransponowanie w buforze wartości B B. Ta metoda daje najszybszy wynik (mnożenie macierzy idzie jako O (n^3) i transpozycja jako O (n^2), więc wykonanie transpozycji jest co najmniej 1000x szybsze). Metoda wiki bez blokowania jest również szybka i nie wymaga bufora. Metoda blokowania jest wolniejsza. Innym problemem z blokowaniem jest konieczność aktualizacji bloku kilka razy. Jest to wyzwanie dla wątków/OpenMP, ponieważ może powodować warunki wyścigu, jeśli nie jest ostrożny.

Należy zauważyć, że przy użyciu AVX z modyfikacją metody transpozycji uzyskuję wyniki szybciej niż Eigen. Jednak moje wyniki z SSE są nieco wolniejsze niż Eigen, więc myślę, że lepiej byłoby użyć buforowania.

Edytuj: Myślę, że mam pomysł, co chcę zrobić. To wynika z algorytmu Cannona od 1969
http://en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms

trzeba wykonać mnożenie macierzy bloku i ma każda nić uchwyt podpotok matrycy C zamiast i B . Na przykład, jeśli podzielę moje macierze na cztery bloki. Zrobiłbym:

[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 

Który usuwa stan wyścigu. Muszę o tym pomyśleć.

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

Najlepsze wyniki stałam się dodając jeszcze jedną pętlę, która blokuje for nad N, a poprzez zmianę pętle. Podniosłem również niezmienny kod, ale optymalizator kompilatora powinien mieć nadzieję, że zrobi to automatycznie. Rozmiar bloku powinien być wielkością linii pamięci podręcznej podzielonej przez sizeof(float). To było ~ 50% szybsze niż transponowane podejście.

Jeśli musisz wybrać tylko jeden z AVX lub blokowanie, korzystanie z rozszerzeń AVX (vfmadd###ps i haddps) jest nadal znacznie szybsze. Używanie obu jest najlepsze i łatwe do dodania, biorąc pod uwagę, że już testujesz, jeśli rozmiar bloku jest wielokrotnością 64/sizeof(float) == 16 zmiennych == dwa 256-bitowe rejestry AVX.

• Przeniesione: 1816522 kleszcze
• planszy: 892431 (51% szybciej)
• AVX + płytki: 230512 (87% szybciej)

Okładziny:

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]; 
         } 
        } 
       } 
      } 
     } 
    } 
} 

---

Jeśli chodzi o odniesienie do armaty, SUMMA algorithm jest lepszy do naśladowania.

---

W przypadku gdy ktokolwiek inny jest optymalizacja Tall-chudy mnożenia ({~ 1E9 x 50} x {50 x 50}, jak skończyło się tutaj), transpozycja podejście jest prawie identyczny w działaniu do zablokowanego podejścia do n = 18 (pływaki). n = 18 to przypadek patologiczny (znacznie gorszy niż 17 lub 19) i nie widzę w ogóle wzorców dostępu do pamięci podręcznej, które to powodują. Wszystkie większe n są ulepszone z zablokowanym podejściem.