Làm nhiều phép nhân ma trận-ma trận trong một hoạt động

Question

Tôi đang thực hiện một thuật toán đó, về bản chất, là một loạt các phép nhân ma trận ma trận như thế này:

 
Res = M1.M2.M3. ... .Mn 

ma trận của tôi là nổi 100x100 thực sự nhỏ, nhưng trình tự thực sự dài, theo thứ tự của hàng tỷ.

Tôi đã thử sử dụng CUBLAS cho các phép nhân nhưng điều này chậm, tuy nhiên tôi đã chú ý điều gì đó thú vị.

nhân 100x100 với ma trận 100x100 chậm, nhưng nhân 1.000.000x100 với 100x100 tương đối nhanh, điều này khiến tôi nghĩ .Nếu thay vì quét từ trái sang phải có 10.000 lần quét song song. Điều này nên được khá nhanh, và nếu tôi nhân ma trận của tôi khi tôi đã được thực hiện với điều này, tôi sẽ nhận được kết quả tương tự - chỉ nhanh hơn.

 
Res1 = M1.M2.M3. ... .Mn/1000-1 
Res1 = M1+n/1000.M2+n/1000.M3+n/1000. ... .M2(n/1000)-1 
... 
Res1 = M1+999*n/1000.M2+999*n/1000.M3+999*n/1000. ... .M1000*(n/1000)-1 
Res = Res1*Res2* ... *Res999 

giá trị của nó không có gì mà M_1 ... M_n nằm trong một nhóm khoảng 100 ma trận khác nhau, do tiêu thụ không gian là không thực sự là một vấn đề, tất cả tôi cần đến là được làm nhiều nhân lên trong một hoạt động .

Bây giờ, đây là vấn đề của tôi. Tôi đã thực hiện một ma trận ma trận (sgemm) thực hiện lấy cảm hứng từ một nvidia thể hiện trong tài liệu của họ, nhưng nó là một thứ tự khoảng 4 lần chậm như cublas. Có ai biết CUBLAS hoạt động như thế nào không? Và nếu mã có sẵn ở đâu đó?

Answer 1

Bạn đã xem CUBLAS (version 4.1) mới nhất chưa? Nó bao gồm một chế độ GEMM theo đợt mới được thiết kế riêng cho các lô lớn các nhân ma trận ma trận nhỏ.Tôi sẽ đề nghị làm một cây nhân đôi như Jonathan Dursi gợi ý trong câu trả lời của mình, bằng cách sử dụng API theo lô CUBLAS để tăng tốc nó, thay vì viết kernel tùy chỉnh của riêng bạn như anh ta gợi ý.

CUBLAS 4.1 được bao gồm trong CUDA Toolkit v4.1.

Answer 2

Vấn đề là cublas, vv được thiết kế để sử dụng tất cả các SM để nhân các ma trận lớn. Đó không phải là điều bạn muốn; bạn muốn thực hiện rất nhiều phép nhân ma trận nhỏ.

Có thể có cách nào đó để truyền nội dung này vào thứ gì đó CUBLAS có thể làm tốt cho bạn, nhưng tôi không nhìn thấy nó. Đề xuất của tôi sẽ như sau:

Viết hạt nhân sử dụng một khối chuỗi để nhân hai ma trận nhỏ của bạn và xuất kết quả.

Sau đó khởi động log kernel N với tấn khối và giải quyết những cặp nhân:

• Bước 1: nhân M x M , M x M ... M _{N - 2} x M _N-1, xuất ra M ', M' .. M ' _N/2
• Bước 2: nhân M' x M ', M' x M ' ... M' _{N/2 - 2} x M _N/2-1, xuất ra M '' , M '' .. M '' _N/4 ...

, vv

Sẽ có hệ số 50% phí trên bộ nhớ, nhưng tôi nghĩ bạn sẽ tận dụng tốt hơn lõi của bạn theo cách này.

Cập nhật

Ok, nếu bạn thực sự không muốn làm điều này trong giai đoạn, bạn có thể làm theo cách này nhưng nó sẽ đòi hỏi mã hóa hơn, và hiệu suất có lẽ sẽ tồi tệ hơn so với những gì bạn có thể nhận được một cái gì đó như cuBLAS và chuyển giao không đồng bộ. Tôi giả sử bạn đang sử dụng một Fermi, và bạn đã tắt bộ nhớ cache L1 để bạn có 48K chia sẻ mem cho mỗi SM.

Lưu trữ 100 ma trận ở dạng khối 2x2, với mỗi khối có liên quan trong bộ nhớ. Vì vậy, matrix[matrixnum,i,j] bắt đầu tại matricies[matrixnum*100*100 + i*100*50 + j*50*50]. Lưu ý rằng mỗi khối là 50 * 50 * 4 byte ~ 10K, do đó, 4 thoải mái phù hợp trong bộ nhớ chia sẻ.

Chỉ định mỗi chuỗi 4 chuỗi một (Nmatricies/Nblocks) chuỗi dài của ma trận để nhân lên, với một trong bốn chuỗi chịu trách nhiệm cho mỗi khối phép nhân.

Giả sử bạn đang tạo chuỗi 1 trong 4 và phần đầu tiên của ma trận bạn nhân với nhau là AxB. Bạn có trách nhiệm (1,1) của kết quả - (AB) _1,1 = A _1,1 B _1,1 + A _1,2 * B _2,1. Bạn tải trước trong A _1,1 vào khoanh vùng [0] trong bộ nhớ dùng chung.

• tải trong myblock [1] = B _1,1 từ bộ nhớ toàn cầu
• myblock [3] = myblock [0] * myblock [1] (ma trận mult, tất cả trong bộ nhớ chia sẻ)
• tải trong myblock [1] = A _1,2 từ toàn cầu
• tải trong myblock [2] = B _2,1 từ toàn cầu
• myblock [0] = myblock [3] + (myblock [ 1] * myblock [2]) (ma trận đa và bổ sung, tất cả trong bộ nhớ chia sẻ).

Bây giờ bạn có thể lặp lại điều này cho phần còn lại của chuỗi ma trận trong phần của chuỗi, chỉ xuất hiện khi hoàn tất.

Khi bạn hoàn thành, bạn sẽ kết thúc với (#SMs) ma trận trong bộ nhớ toàn cầu, vẫn phải nhân lên, nhưng sẽ không có thêm bộ nhớ tạm thời nào trong bộ nhớ toàn cầu và bạn đã thắng ' t đã phải sao chép dữ liệu vào bộ nhớ toàn cầu khác với các ma trận gốc và danh sách của những cái để giải quyết.

Một lần nữa, không có lý do thực sự để làm điều này ngoại trừ việc bạn không thể bị làm phiền khi truyền dữ liệu đến GPU theo từng giai đoạn và hiệu suất gần như chắc chắn sẽ tồi tệ hơn; có ít bộ nhớ toàn cục hơn, nhưng bạn có thể sẽ trả tiền cho nó bằng một GEMM được cuộn bằng tay. Tin tốt là 50 không phải là bội số của 8, vì vậy bạn có thể sẽ không có quá nhiều trong cách thức xung đột ngân hàng bộ nhớ chia sẻ.

Một lần nữa, đối với điểm thưởng, bạn có thể tính toán trước tất cả các khối của tất cả các sản phẩm ma trận ghép nối trước và sau đó bằng một nửa độ dài của danh sách.

Answer 3

LIBXSMM - một thư viện nhắm mục tiêu Kiến trúc Intel cho phép nhân ma trận nhỏ, dày đặc hay thưa thớt, và nhiều nếp cuộn nhỏ được chính xác có nghĩa là để khai thác hiệu quả tốt nhất cho phép nhân ma trận nhỏ.

Ngược lại với NVidia CUBLAS (hoặc Intel MKL), LIBXSMM không dựa trên giao diện hàng loạt. Thay vào đó, người ta có thể sắp xếp cho các cuộc gọi cá nhân và cũng cung cấp "các vị trí tiếp theo" tức là, nơi các toán hạng/ma trận của các phép nhân tiếp theo được đặt (trong bộ nhớ). Ưu điểm là không cần cấu trúc dữ liệu rõ ràng hoặc định dạng chỉ mục mô tả lô.

#include <libxsmm.h> 

int main() 
{ 
    const libxsmm_gemm_prefetch_type prefetch = LIBXSMM_PREFETCH_AUTO; 
    const double alpha = 1.0, beta = 1.0; /* accumulate C-matrix */ 
    const int m = 23, n = 23, k = 23;  /* some problem size */ 
    libxsmm_dmmfunction xmm = NULL;  /* function pointer */ 

    xmm = libxsmm_dmmdispatch(23, 23, 23, /* some problem size */ 
      NULL/*lda*/, NULL/*ldb*/, NULL/*ldc*/, 
      &alpha, &beta, NULL/*flags*/, 
      NULL/*&prefetch*/); 

    if (xmm) { /* JiT'ted code has been generated */ 
# pragma omp parallel for 
    for (int i = 0; i < nbatch; ++i) { 
     const double *const ai = a + i * asize; 
     const double *const bi = b + i * bsize; 
     /* e.g., matrix C is accumulated (instead of streamed) */ 
     double *const ci = c /*+ i * csize*/; 
     /* optionally provide "next locations" */ 
     xmm(ai, bi, ci/*, 
      ai + 1 * asize, 
      bi + 1 * bsize, 
      ci + 0 * csize 
     */); 
    } 
    } 
} 

LIBXSMM sản xuất tối ưu hóa cao và mã chuyên ngành (JIT), trong đó khai thác mới nhất phần mở rộng tập lệnh (SSE3, AVX, AVX2, và AVX-512). LIBXSMM là available theo giấy phép không cho phép (mệnh đề BSD-3).

LƯU Ý: Đây không phải về CUBLAS (như được yêu cầu ban đầu).