Vượt hiệu năng nhân ma trận NumPy chỉ với 150 dòng mã C

 (salykova.github.io)
1 điểm bởi GN⁺ 2024-07-05 | 1 bình luận | Chia sẻ qua WhatsApp
  • Phép nhân ma trận của NumPy dựa vào thư viện BLAS bên ngoài, nhưng cách triển khai này nhắm tới việc đưa hiệu năng đơn luồng và đa luồng lên mức BLAS chỉ bằng C thuần cùng FMA3 và AVX2
  • Trọng tâm hiệu năng là chia $C$ thành các khối nhỏ, trong đó microkernel 16×6 lặp lại rank-1 update bên trong các thanh ghi YMM để giảm truy cập bộ nhớ
  • Với ma trận kích thước bất kỳ, xử lý biên dễ trở thành nút thắt cổ chai, nên giải pháp kết hợp mask store và bộ đệm padding 0 để tránh suy giảm hiệu năng của mask load
  • Khả năng tái sử dụng cache được đảm bảo bằng blocking k_c, m_c, n_c, và hiệu năng đỉnh thực tế phụ thuộc đáng kể vào việc tinh chỉnh số luồng, kích thước kernel và kích thước tile
  • AVX-512 bị loại trừ để hỗ trợ CPU rộng hơn, nên trên CPU có AVX-512 thì BLAS có thể nhanh hơn; việc so sánh với OpenBLAS cũng được thực hiện trong điều kiện tắt AVX-512

Mục tiêu triển khai và đối tượng so sánh

  • Mã triển khai được công khai tại sgemm.c và tối ưu hóa phép nhân ma trận FP32 đa luồng trên các bộ xử lý hiện đại
  • NumPy dựa vào các thư viện BLAS bên ngoài cho các phép toán đại số tuyến tính như nhân ma trận
    • Ví dụ gồm Intel MKL, Accelerate, BLIS, GotoBLAS, OpenBLAS
    • OpenBLAS, GotoBLAS, BLIS được viết bằng C/FORTRAN/Assembly và bao gồm các triển khai nhân ma trận tối ưu thủ công theo từng vi kiến trúc CPU
  • Mục tiêu là một triển khai nhân ma trận được viết bằng C thuần không cần assembly mức thấp nhưng vẫn đáp ứng các điều kiện sau
    • Hoạt động với ma trận có kích thước bất kỳ
    • Chạy được trên các bộ xử lý x86-64 hiện đại
    • Cạnh tranh được với các thư viện BLAS hiện có
    • Mã đơn giản và dễ mở rộng
  • Tài liệu tham khảo gồm Fast Multidimensional Matrix Multiplication on CPU from Scratch của Simon Boehm, Matrix Multiplication của Sergey Slotin, Can you multiply a matrix? của Geohot, cùng các bài báo về GotoBLAS và BLIS

Điều kiện benchmark và cách tính FLOPS

  • Môi trường thử nghiệm là AMD Ryzen 7 9700X, 32GB DDR5 6000 MHz CL36, OpenBLAS 0.3.26, GCC 13.3, Ubuntu 24.04.1 LTS
  • Cờ biên dịch sử dụng là -O3 -march=native -mno-avx512f -fopenmp
  • Để so sánh công bằng, khi cài OpenBLAS cần đặt TARGET phù hợp và vô hiệu hóa tập lệnh AVX-512
    • Với bộ xử lý Zen4/5, biên dịch bằng make TARGET=ZEN
    • Nếu không, OpenBLAS sẽ mặc định dùng lệnh AVX-512
  • Phép nhân ma trận FP32 của OpenBLAS được chạy qua API cblas_sgemm
  • Benchmark được thực hiện trên ma trận vuông
    • Đánh giá từ m=n=k=200 đến m=n=k=10000 với bước nhảy 200
    • Phép nhân ma trận được lặp n_iter lần và thời gian chạy trung vị được dùng để đo hiệu năng
  • Khi nhân ma trận $A$ kích thước $M \times K$ với ma trận $B$ kích thước $K \times N$, tổng lượng tính toán là $2MNK$ FLOP
    • Hiệu năng được tính bằng FLOPS=(2*m*n*k)/exec_time

Giới hạn lý thuyết và nền tảng SIMD

  • CPU x86-64 hiện đại xử lý song song nhiều dữ liệu nhờ mở rộng SIMD
  • Các tập lệnh chính là AVX2 và FMA
    • Cả hai đều dùng thanh ghi YMM 256 bit
    • Mỗi thanh ghi YMM chứa được 8 số thực float 32 bit
  • Lệnh FMA VFMADD231PS thực hiện phép toán packed single theo dạng YMM1 = YMM2 * YMM3 + YMM1
  • Trên Ryzen 9700X, throughput của fused multiply-add là 0.5 cycles/instruction, tức 2 lệnh mỗi chu kỳ
  • Về lý thuyết, Ryzen 9700X có thể thực hiện 32 FLOP mỗi chu kỳ trên một lõi
    • Công thức là 8 floats × 2(add+mul) × 2(1/TP)
    • Với 8 lõi và xung duy trì 4.7GHz, đỉnh lý thuyết đa luồng được ước tính là 1203 FLOPS

Triển khai cơ bản và microkernel

  • Ma trận được lưu theo thứ tự column-major
    • A[row][col] được truy cập trong con trỏ C dưới dạng ptr[col*M + row]
  • Triển khai đơn giản nhất là duyệt mọi hàng và cột của $C$, rồi tính tích vô hướng giữa hàng của $A$ và cột của $B$ cho từng phần tử
  • Cốt lõi của triển khai hiệu năng cao là microkernel, chia $C$ thành các ma trận con $m_R \times n_R$ và tính từng phần một cách hiệu quả
  • Kernel khởi tạo $\bar{C}$ bằng 0 trong thanh ghi rồi lặp theo chiều $K$
    • Nạp vector cột của $\bar{A}$ và vector hàng của $\bar{B}$ vào thanh ghi
    • Tính tích ngoài của hai vector rồi cộng vào bộ tích lũy $\bar{C}$
    • Mỗi bước là một rank-1 update
  • So với lượng truy cập bộ nhớ $2K m_R n_R$ của cách naive, phương pháp này giảm số phần tử phải nạp vào thanh ghi xuống còn $(m_R+n_R)K$
  • CPU AVX có 16 thanh ghi YMM, nên kích thước kernel phải thỏa các ràng buộc sau
    • $(m_R/8) \cdot n_R + m_R/8 + 1 \le 16$
    • $m_R$ phải là bội số của 8
  • Về lý thuyết, $m_R$ và $n_R$ càng lớn và càng gần nhau thì càng giảm được truy cập bộ nhớ, nhưng trên Ryzen 9700X thực tế kernel 16×6 cho hiệu năng tốt nhất
  • Triển khai dùng intrinsic trong immintrin.h
    • __m256 là kiểu vector 256 bit, biểu diễn nội dung thanh ghi YMM
    • _mm256_loadu_ps dùng để nạp vector cột của A
    • _mm256_broadcast_ss dùng để broadcast một giá trị vô hướng của B thành vector 8 float
    • _mm256_fmadd_ps dùng để cập nhật bộ tích lũy
    • _mm256_storeu_ps dùng để ghi kết quả ra bộ nhớ
  • Assembly sinh ra có các lệnh SIMD FMA như vfmadd231ps và vbroadcastss

Padding cho ma trận kích thước bất kỳ

  • Kernel 16×6 cơ bản hoạt động trực tiếp khi $M$ và $N$ lần lượt là bội số của 16 và 6
  • Ở vùng biên, nếu số cột $n$ nhỏ hơn 6 thì vòng lặp ghi chỉ chạy đến j < n
  • Khi số hàng $m$ nhỏ hơn 16, _mm256_storeu_ps ghi 8 phần tử cùng lúc nên cần mask store
    • _mm256_maskstore_ps chỉ ghi các phần tử có bit mask được bật vào bộ nhớ
    • Mask được tạo theo số hàng chồng lấp $m$
  • Nếu cả phần nạp ở vùng biên cũng xử lý bằng _mm256_maskload_ps thì hiệu năng kernel có thể giảm mạnh
    • Lệnh bổ sung để tính mask tạo ra overhead
    • Vì $n$ không phải hằng số tại thời điểm biên dịch, compiler khó unroll vòng lặp hiệu quả
  • Thay vào đó, nếu $m \neq m_R$ thì sao chép $\bar{A}$ vào bộ đệm và padding bằng 0; nếu $n \neq n_R$ thì cũng sao chép $\bar{B}$ vào bộ đệm và điền 0
  • Phần triển khai liên quan nằm tại matmul_pad.h

Cache blocking và tái sử dụng dữ liệu

  • Giữa thanh ghi và DRAM có hệ thống phân cấp CPU cache; CPU desktop hiện đại thường dùng cache L1, L2 và L3
  • Cache nhanh hơn DRAM nhưng dung lượng có hạn, nên không thể chứa toàn bộ $A$, $B$, $C$ trong cache cùng lúc
  • Chia ma trận thành các khối nhỏ để đưa vào cache và tái sử dụng cùng dữ liệu cho nhiều rank-1 update được gọi là cache blocking hoặc tiling
  • Cache blocking đơn luồng có dạng 5 vòng lặp tương tự cấu trúc của BLIS
    • Vòng ngoài cùng tạo các khối $C_j$ và $B_j$ theo chiều $N$
    • Vòng tiếp theo tạo các khối $A_j$ và $B_p$ theo chiều $K$
    • $B_p$ được packing thành $\tilde{B}_p$ và, khi cần, được padding bằng 0 để tận dụng tái sử dụng trong cache L3
    • Vòng tiếp theo tạo các khối $C_i$ và $A_j$ theo chiều $M$, trong đó $A_j$ được packing thành $\tilde{A}_j$
    • Hai vòng cuối chia các khối cache thành các panel $m_R \times k_c$ và $k_c \times n_R$ để đưa vào kernel
  • Các khối đã packing $\tilde{A}_j$ và $\tilde{B}_p$ được lưu theo cách khác nhau
    • Các panel bên trong $\tilde{A}_j$ được lưu theo column-major
    • Các panel bên trong $\tilde{B}_p$ được lưu theo row-major
  • Các tham số cache blocking cần được điều chỉnh theo kích thước cache của từng mẫu CPU
    • $k_c \times n_c$ là điểm khởi đầu để lấp đầy cache L3
    • $m_c \times k_c$ là điểm khởi đầu để lấp đầy cache L2
    • $k_c \times n_R$ là điểm khởi đầu để lấp đầy cache L1
  • Trong thực tế, giá trị lớn hơn lý thuyết thường cho hiệu năng tốt hơn, và CPU tự quản lý bố trí cache, nên ở mức thuật toán cần thiết kế vòng lặp và mẫu truy cập phù hợp
  • Phần triển khai nằm tại matmul_cache.h

Tối ưu vi mô cho kernel

  • Thay vì định nghĩa bộ tích lũy bằng mảng như __m256 C_buffer[6][2], mã khai báo tường minh từng biến tích lũy
  • Cách này giúp GCC tối ưu mã tốt hơn và tránh register spilling
  • Phần tính mask cũng được đổi sang dùng lệnh vector
    • Dùng mảng tĩnh mask[32] cùng _mm256_cvtepi8_epi32 và _mm_loadu_si64
  • Phần triển khai này nằm tại matmul_micro.h

Chiến lược đa luồng

  • Cả phần tính toán số học lẫn packing đều được song song hóa
  • Vòng lặp thứ 5, 4 và 3 bên ngoài microkernel lặp theo đơn vị kích thước khối cache
    • Để giữ mọi luồng luôn bận, số lần lặp cần lớn hơn hoặc bằng số luồng
    • Kích thước ma trận đầu vào cần xấp xỉ từ số luồng × kích thước khối cache trở lên
  • Trên một luồng của Ryzen 9700X, kích thước khối cache cho hiệu năng tốt là $n_c=1535$, $m_c=1024$
    • Để tận dụng đủ 8 lõi cần ít nhất chiều kích thước khoảng $\max(m_c,n_c) \times 8 = 1535 \times 8 = 12280$
  • Ngược lại, hai vòng cuối lặp trên các khối nhỏ $m_R$, $n_R$ nên rất phù hợp để song song hóa
    • Thông thường $m_R$, $n_R$ nhỏ hơn 20
    • Nếu chọn $m_c$, $n_c$ là bội số của số lõi thì có thể chia việc đều hơn
  • Trên Ryzen 9700X, cách cho hiệu năng tốt nhất là song song hóa đồng thời hai vòng lặp bên trong bằng #pragma omp parallel for collapse(2) num_threads(NTHREADS)
  • Với bộ xử lý nhiều lõi, đặc biệt hơn 16 lõi, có thể cân nhắc song song lồng nhau và song song hóa 2 đến 3 vòng lặp
  • Việc packing $\tilde{A}$ và $\tilde{B}$ cũng được song song hóa bằng OpenMP
    • pack_blockA được song song hóa khi duyệt mc theo đơn vị MR
    • pack_blockB được song song hóa khi duyệt nc theo đơn vị NR
  • Trong triển khai đa luồng, các tham số cho hiệu năng tốt trên Ryzen 9700X là
    • $m_c = m_R \times \text{number of threads} \times 5$
    • $n_c = n_R \times \text{number of threads} \times 50$
  • Triển khai đa luồng hoàn chỉnh nằm tại matmul_parallel.h

Bài viết liên quan

1 bình luận

 
GN⁺ 2024-07-05
Ý kiến trên Hacker News

  • Nếu ý chính của bài viết này là thường thì vẫn còn dư địa hiệu năng, thì thực ra nó còn đánh giá thấp mức cải thiện. Dù công sức đổ vào các thư viện nhân ma trận lớn hơn rất nhiều so với phần lớn phần mềm khác
    Nếu không phải là mã đã được tối ưu hóa mạnh, việc cải thiện 10~1000 lần trở lên từ mã hiện có mà không cần quá nhiều công sức là chuyện thường gặp. Xét đại khái theo thứ tự quan trọng, quan trọng nhất là lựa chọn thuật toán có phù hợp không và có thể loại bỏ hẳn công việc đó không; việc giảm các thao tác nặng như đi-về kernel hay malloc cũng có tác động lớn.
    Về vector hóa, các intrinsic vector tường minh cũng tốt, nhưng nhiều khi chỉ cần tái cấu trúc dữ liệu từ mảng struct sang mảng/các mảng của struct là đã sinh ra cùng mã máy. Hiệu quả cache cũng quan trọng, và trong mã song song thì phức tạp hơn khi dữ liệu theo từng thread không được cô lập, như hiện tượng false sharing. Cuối cùng cũng có thể tối ưu theo phần cứng, chẳng hạn intrinsic hoặc assembly viết tay

    • Cũng không thể bỏ qua tác động của mạng. Có lần tôi phát hiện một truy vấn phân tán kéo khoảng 1 triệu hàng qua mạng rồi join để cuối cùng chỉ còn 5~10 hàng, và việc đó đã cải thiện hiệu năng hàng trăm lần
      Chỉ cần đổi truy vấn để phép join diễn ra trên máy chủ từ xa và qua mạng chỉ gửi 5~10 hàng là nhanh lên ngay. Overhead cố định và độ trễ thì luôn có, nhưng nếu gửi qua kết nối mạng nhiều dữ liệu hơn hẳn mức cần thiết thì cuối cùng hiệu năng sẽ hỏng. Bài “It's the latency, stupid” nói về tác động của độ trễ cũng đáng đọc: http://www.stuartcheshire.org/rants/latency.html
      Nhìn chung tôi đồng ý với các điểm cân nhắc ở trên và thứ tự tương đối của chúng
    • “Lựa chọn thuật toán có phù hợp không” trên thực tế phần nào đã trở thành cargo cult. Nhiều thuật toán “nhanh hơn” có hằng số thực tế tệ đến mức, phía làm nhiều việc hơn lại có hiệu năng tốt hơn
      Nhiều buổi phỏng vấn, thay vì xem cách suy luận vì sao triển khai chậm, benchmark và sửa nó, lại biến thành các câu đố học thuộc thuật toán obscure kiểu “vì Google làm thế”

  • Các mẫu lập trình phổ biến không đủ đặc thù hóa cho phần cứng nên để lại rất nhiều hiệu năng chưa khai thác. Bài này là một ví dụ thú vị, và một minh họa kinh điển khác là “There's plenty of room at the top”
    https://www.science.org/doi/10.1126/science.aam9744

  • Để hiểu nội dung này, các bài báo trong kho BLIS gần như là nguồn chuẩn. Tôi không hiểu vì sao lại cho rằng BLAS đã tối ưu mà hiệu năng không tốt; với ma trận đủ lớn thì nên kỳ vọng trên 90% đỉnh CPU
    Lần cuối tôi xem, OpenBLAS chạy tuần tự nhìn chung tương tự MKL, và BLAS triển khai GEMM chứ không phải matmul như khối tuyến tính đại số cơ bản. Tôi cũng không hiểu vì sao lại dùng numpy thay vì một framework benchmark thông thường; trên Zen thì nên so với BLAS của AMD, tức triển khai dựa trên BLIS. Trước đây BLIS có câu chuyện về song song hóa tốt hơn OpenBLAS, và AMD BLIS còn có chuyển đổi triển khai cho các chiều “nhỏ”; tôi không rõ hiện OpenBLAS có chưa
    Vector hóa micro-kernel không nhất thiết cần intrinsic SIMD; một trình biên dịch C khá tốt sẽ vector hóa hoàn toàn và cũng unroll vòng lặp. Micro-kernel C thuần của BLIS đạt hơn 80% hiệu năng so với triển khai tối ưu thủ công cho Haswell ở kích thước khối phù hợp. Khác biệt có lẽ là do prefetch, nhưng tôi không hiểu chính xác

    • Intrinsic SIMD và unroll vòng lặp thủ công chắc chắn là cần thiết. Đó là lý do mọi thư viện BLAS đều vector hóa và unroll vòng lặp thủ công
      Các trình biên dịch hiện đại cũng không thể tự động vector hóa và unroll vòng lặp đúng với tỷ lệ thành công 100%

  • Bài viết và phần triển khai trông ổn, nhưng tôi tò mò “bí quyết” là gì. OpenBLAS đã được tối ưu bằng assembly+C suốt hàng chục năm cho đúng bài toán này, vậy làm sao có thể thắng được?
    Bài viết nói khá chi tiết về caching, v.v., nên tôi thắc mắc liệu BLAS không tận dụng những thứ này, hay là phần này được tinh chỉnh tốt hơn cho một bộ xử lý cụ thể

    • OpenBLAS không được tối ưu đến mức đó cho một số kiến trúc mới nhất định. Ma trận cũng không lớn đến vậy, và numpy có overhead cffi
      Chênh lệch hiệu năng nổi bật hơn nhiều ở thông lượng đỉnh so với thông lượng trung bình, mà rất ít ứng dụng thật sự cần thông lượng đỉnh. Mã benchmark được hiển thị có vẻ như phía numpy đi qua allocator của Python còn bản triển khai C thì không đi qua allocator, nên đây là chỗ đầu tiên cần kiểm tra lỗi hoặc sự không nhất quán của microbenchmark. Nhiều routine của numpy hỗ trợ thao tác tại chỗ, vì vậy có lẽ cần xem rõ benchmark phiên bản in-place cho cả hai bên
      numpy cũng có kiểm tra biên và xử lý lỗi chạy độc lập với phần triển khai bên dưới, khiến nó rất chậm, thậm chí chậm hơn cả list Python thuần với các ma trận nhỏ. Nếu cộng thêm vài nghìn cycle overhead thuần thì rất khó làm cho nhanh
      Bản triển khai này là một cách tiếp cận khá có nguyên tắc nhằm bão hòa các cache liên quan, và theo một nghĩa nào đó là hiển nhiên, nhưng những cải tiến kỹ thuật rõ ràng như vậy đáng được nhấn mạnh trong các thảo luận kiểu này. OpenBLAS cũng đã đầu tư nhiều nhân lực, nhưng khó có khả năng họ đã nghĩ hết mọi thứ. Muốn giải thích đúng thì cần phân tích sâu mã của cả hai bên
    • Việc đánh bại OpenBLAS không hề đáng ngạc nhiên và cũng không phải chưa từng có tiền lệ. Chẳng hạn thư viện đại số tuyến tính Mir của ngôn ngữ D cũng đã làm được điều đó vài năm trước [1]
      Với các bản triển khai C++ và C, hãy xem cách tiếp cận metaprogramming [2], [3]. Điều thật sự đáng ngạc nhiên là nhiều ngôn ngữ hiện đại như Matlab, Julia, Mojo vẫn còn dựa vào OpenBLAS, dù tất nhiên mỗi bên đều có lý do riêng
      [1] Numeric age for D: Mir GLAS is faster than OpenBLAS and Eigen (2016):
      http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...
      [2] Vastly outperforming LAPACK with C++ metaprogramming (2018):
      https://wordsandbuttons.online/vastly_outperforming_lapack_w...
      [3] Outperforming LAPACK with C metaprogramming (2018):
      https://wordsandbuttons.online/outperforming_lapack_with_c_m...
    • -march=native có thể là một lợi thế vì nó biên dịch đúng theo model CPU cụ thể. numpy nhiều khả năng được biên dịch cho mục tiêu x86-64 cũ hơn và tổng quát hơn
      Trên CPU Ryzen, -march=native có lẽ dùng v4, còn numpy có vẻ nhắm tới v1 hoặc v2
      https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...
    • numpy 2.0 tích hợp Google highway để dùng SIMD tốt hơn trên nhiều vi kiến trúc, nên phép so sánh phía numpy sẽ khá hơn

  • Bài viết hay, và việc làm cho benchmark dễ tái hiện cũng rất tuyệt. Trên Xeon W-2245 16 nhân 3.90GHz của tôi, matmul.c thực hiện phép nhân ma trận 8192x8192 trong 1,41 giây với gcc -O3, 1,47 giây với clang -O2, còn NumPy là 1,07 giây
    Tôi nghĩ nếu dùng kernel AVX-512 thì sẽ nhanh hơn nhiều. Một lý do khác khiến hiệu năng chưa như ý có thể là OpenMP; theo kinh nghiệm, nếu quản lý thread pool rõ ràng bằng pthreads thì có thể giảm overhead. Thay vì hardcode số CPU, cũng nên dùng sysconf(_SC_NPROCESSORS_ONLN)

  • Không có lý do gì để một bên là Python còn bên kia là C, khiến gánh nặng khác nhau như vậy. Có thể viết cả hai bên bằng C, một bên gọi thư viện BLAS, bên kia gọi bản triển khai này, như vậy mới là so sánh ngang bằng

    • Ở đây so sánh với Python là đúng. Vì cách phổ biến nhất hiện nay để thực hiện các phép tính như thế này là dùng Python với numpy
      Overhead không quá lớn, nhưng như đã nói ở những chỗ khác trong thread này, gọi cho đúng là rất quan trọng. Đem mã numpy ngây thơ đấu với mã C đã được tinh chỉnh rõ ràng không phải là một so sánh công bằng

  • Không phải là hot path, nhưng sự kém hiệu quả trong việc tạo mask, cụ thể là việc dùng bit_mask, khiến tôi thấy hơi khó chịu. Cách hiệu quả hơn là tạo một mảng hằng toàn cục dạng {-1,-1,...,0,0,...} rồi load từ offset phần tử 16-m, 8-m, hoặc so sánh vector hằng {0,1,2,3,4,...} với m và m-8 đã được broadcast
    Tuy nhiên, nó chỉ tương ứng với một cột của ma trận, và vòng lặp maskload/maskstore phía sau tốn thời gian hơn nhiều, nên đây chỉ là một lời bắt bẻ rất nhỏ. Đặc biệt thao tác store vẫn chậm ngay cả trên Zen 4[1], còn lệnh AVX-512 thì nhanh hơn 6 lần dù khác biệt chỉ là nhận mask từ thanh ghi mask. Dù sao clang cũng tự động vector hóa phép shift, nên có lẽ chỉ chậm hơn đề xuất của tôi khoảng 2–3 lần
    [1]: https://uops.info/table.html?search=vmaskmovps&cb_lat=on&cb_...

    • Tôi là tác giả bài viết. Đây thực sự là lần đầu tôi tối ưu mã C và dùng intrinsics, nên tôi không phải chuyên gia trong lĩnh vực này, nhưng rất muốn học thêm
      Tôi rất biết ơn phản hồi đem lại góc nhìn mới. “Tạo mảng hằng toàn cục rồi load” theo trí nhớ của tôi khi thử thì có vẻ chậm hơn một chút so với shift bit mask, nhưng tôi sẽ thử lại để chắc chắn. Cách “so sánh vector hằng {0, 1, 2, 3, 4, ...} với m và m-8 đã được broadcast” là một ý tưởng hay, tôi sẽ thử xem
    • Khi tạo mảng hằng toàn cục, có thể để các phần tử là int8_t, rồi khi load thì mở rộng dấu byte thành int32_t. Tổ hợp _mm_loadu_si64 / _mm256_cvtepi8_epi32 sẽ được biên dịch thành một lệnh vpmovsxbd duy nhất dùng toán hạng bộ nhớ
      Như vậy, nếu căn chỉnh đúng bằng alignas(32), toàn bộ mảng hằng sẽ nằm gọn trong một cache line. Với use case trong bài gốc cần hai mask, nên lệnh vpmovsxbd thứ hai chắc chắn sẽ trúng cache L1D, khá phù hợp

  • tinyBLAS của jart thì sao
    https://hacks.mozilla.org/2024/04/llamafiles-progress-four-m...
    Và https://justine.lol/matmul/

    • Hôm qua tôi đã trao đổi khá sôi nổi với Justine, và trên workstation đó có vẻ implementation này nhanh hơn tinyBLAS ít nhất 2 lần. Toàn bộ cuộc thảo luận nằm trên Mozilla AI Discord: https://discord.com/invite/NSnjHmT5xY

  • Ngoài benchmark ra, lý do để đa luồng hóa chính phép nhân ma trận là gì nhỉ. Trên thực tế chẳng phải sẽ có lợi hơn nếu dùng đa luồng ở phía thuật toán sử dụng phép nhân sao?

    • Trong HPC thì thực tế thường làm như vậy. Tuy nhiên, chỉ cần thay bằng BLAS song song cũng có thể dễ dàng giúp ích cho một số loại mã R nhất định
      Nhưng trong mã HPC, thường thì GEMM không phải nút thắt cổ chai

  • Tôi mới chỉ đọc lướt, nhưng bài này có nhiều chi tiết và giải thích. Nó trông như một bài viết khá hay giải thích cách phép nhân ma trận nhanh được triển khai ra sao khi phản ánh các cân nhắc về kiến trúc, nên tôi đã đưa vào danh sách đọc