loop tiling/blocking cho phép nhân số lớn dày đặc

Question

Tôi đã tự hỏi liệu có ai đó có thể chỉ cho tôi cách sử dụng vòng lặp/ngăn chặn vòng lặp cho phép nhân ma trận dày đặc lớn một cách hiệu quả hay không. Tôi đang làm C = AB với ma trận 1000x1000. Tôi đã làm theo ví dụ trên Wikipedia cho tiling vòng lặp nhưng tôi nhận được kết quả tồi tệ hơn bằng cách sử dụng ốp lát hơn là không có.

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

tôi đã cung cấp một số mã dưới đây. Phương pháp ngây thơ là rất chậm do nhớ cache. Phương pháp chuyển vị tạo ra việc chuyển đổi B vào bộ đệm. Phương pháp này cho kết quả nhanh nhất (phép nhân ma trận đi như O (n^3) và transpose là O (n^2) do đó làm transpose là ít nhất 1000x nhanh hơn). Phương pháp wiki mà không chặn cũng nhanh và không cần bộ đệm. Phương pháp chặn chậm hơn. Một vấn đề khác với việc chặn là nó phải cập nhật khối nhiều lần. Đây là một thách thức đối với luồng/OpenMP vì nó có thể gây ra các điều kiện chủng tộc nếu không cẩn thận.

Tôi nên chỉ ra rằng bằng cách sử dụng AVX với một sửa đổi phương pháp chuyển vị, tôi nhận được kết quả nhanh hơn Eigen. Tuy nhiên, kết quả của tôi với SSE là một chút chậm hơn so với Eigen vì vậy tôi nghĩ rằng tôi có thể sử dụng bộ nhớ đệm tốt hơn.

Chỉnh sửa: Tôi nghĩ mình có ý tưởng mình muốn làm. Nó xuất phát từ thuật toán Cannon từ 1969.
http://en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms

tôi cần phải làm ma trận khối nhân và có mỗi thread xử lý một tiểu ma trận của C hơn Một và B . Ví dụ nếu tôi chia ma trận thành bốn khối. Tôi sẽ làm:

[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 

Điều đó loại bỏ điều kiện chủng tộc. Tôi sẽ phải suy nghĩ về điều này.

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

Kết quả tốt nhất mà tôi đã nhận được bằng cách thêm một người nữa for vòng chặn trên N của bạn, và bằng cách sắp xếp lại các vòng. Tôi cũng hoisted loop-bất biến mã, nhưng trình tối ưu hóa của trình biên dịch hy vọng sẽ làm điều này tự động. Kích thước khối phải là kích thước đường bộ nhớ cache chia cho sizeof(float). Điều này đã nhận nó ~ 50% nhanh hơn so với cách tiếp cận transposed.

Nếu bạn chỉ phải chọn một trong AVX hoặc chặn, sử dụng các phần mở rộng AVX (vfmadd###ps và haddps) vẫn còn nhanh hơn đáng kể. Việc sử dụng cả hai cách tốt nhất và dễ hiểu để thêm rằng bạn đã thử nghiệm nếu kích thước khối là bội số của 64/sizeof(float) == 16 float == hai thanh ghi AVX 256 bit.

• Hoán vị: 1.816.522 ve
• Lát: 892.431 (51% nhanh hơn)
• AVX + ốp lát: 230.512 (87% nhanh hơn)

Lát:

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

---

Đối với tham chiếu Cannon, số SUMMA algorithm là một cách tốt hơn để theo dõi.

---

Trong trường hợp bất cứ ai khác là tối ưu hóa phép nhân cao-gầy ({~ 1e9 x 50} x {50 x 50}, làm thế nào tôi đã kết thúc ở đây), cách tiếp cận hoán là gần như giống hệt về hiệu suất với phương pháp chặn tối đa n = 18 (phao). n = 18 là trường hợp bệnh lý (cách tệ hơn 17 hoặc 19) và tôi không thấy các mẫu truy cập bộ nhớ cache gây ra điều này. Tất cả các số lớn hơn n được cải tiến với phương pháp chặn.