Tối ưu hóa 2x2 ma trận nhân: lắp ráp chậm so với SIMD nhanh

Question

Vấn đề

tôi đang học phép nhân ma trận các thuật toán hiệu suất cao như OpenBLAS hoặc GotoBLAS và tôi đang cố gắng để tái tạo một số kết quả. Câu hỏi này đề cập đến hạt nhân bên trong của một thuật toán nhân ma trận. Cụ thể, tôi đang xem máy tính C += AB, trong đó A và B là ma trận 2x2 loại double ở tốc độ tối đa của CPU của tôi. Có hai cách để làm điều này. Một cách là sử dụng hướng dẫn SIMD. Cách thứ hai là mã trực tiếp vào assembly bằng cách sử dụng thanh ghi SIMD.

Những gì tôi đã xem xét cho đến nay

OpenBLAS Tất cả các giấy tờ có liên quan, các trang web Tất nhiên, nhiều nhiều SO Q & Như đối phó với vấn đề này (quá nhiều để liệt kê), tôi đã biên soạn trên máy tính của tôi, xem qua các mã nguồn OpenBLAS, GotoBLAS và BLIS, các sách hướng dẫn của Agner.

Phần cứng

CPU của tôi là một i5 Intel - 540M. Bạn có thể tìm thấy thông tin CPUID có liên quan trên cpu-world.com. Kiến trúc vi mô là Nehalem (westmere), do đó, về mặt lý thuyết, nó có thể tính toán 4 lần độ lệch chính xác gấp đôi trên mỗi lõi mỗi chu kỳ. Tôi sẽ chỉ sử dụng một lõi (không có OpenMP), vì vậy với siêu phân luồng và Intel Turbo Boost 4 bước, tôi sẽ thấy một đỉnh của (2.533 Ghz + 4*0.133 Ghz) * (4 DP flops/core/cycle) * (1 core) = 12.27 DP Gflops. Để tham khảo, với cả hai lõi chạy ở đỉnh, Intel Turbo Boost tăng tốc 2 bước và tôi sẽ nhận được một đỉnh lý thuyết là 22.4 DP Gflops.

Cài đặt

Tôi tuyên bố ma trận 2x2 của tôi như double và khởi tạo chúng với mục ngẫu nhiên như trong đoạn mã dưới đây.

srand(time(NULL)); 
const int n = 2; 
double A[n*n]; 
double B[n*n]; 
double C[n*n]; 
double T[n*n]; 
for(int i = 0; i < n*n; i++){ 
    A[i] = (double) rand()/RAND_MAX; 
    B[i] = (double) rand()/RAND_MAX; 
    C[i] = 0.0; 
} 

tôi tính toán một câu trả lời đúng bằng ngây thơ ma trận ma trận multiplcation (hình dưới đây) cho phép tôi để kiểm tra kết quả của tôi hoặc là trực quan hoặc bằng cách tính toán các chỉ tiêu L2 của tất cả các yếu tố

// "true" answer 
for(int i = 0; i < n; i++) 
    for(int j = 0; j < n; j++) 
     for(int k = 0; k < n; k++) 
      T[i*n + j] += A[i*n + k]*B[k*n + j]; 

Để chạy mã và nhận được một ước tính của Gflops, tôi gọi mỗi chức năng nhân một lần để khởi động, và sau đó thực hiện nó bên trong một vòng lặp for cho maxiter lần, đảm bảo không có ma trận C mỗi khi tôi đang tính toán C += AB. Vòng lặp for được đặt bên trong hai câu lệnh clock() và điều này được sử dụng để ước tính Gflops. Đoạn mã sẽ minh họa phần này.

C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 
mult2by2(A,B,C); //warmup 
time1 = clock(); 
for(int i = 0; i < maxiter; i++){ 
     mult2by2(A,B,C); 
     C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 
} 
time2 = clock() - time1; 
time3 = (double)(time2)/CLOCKS_PER_SEC; 
gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; 
mult2by2(A,B,C); // to compute the norm against T 
norm = L2norm(n,C,T); 

đang SIMD

CPU của tôi hỗ trợ vectơ 128-bit, vì vậy tôi có thể phù hợp 2 double s trong mỗi vector. Đây là lý do chính tại sao tôi làm phép nhân 2x2 ma trận trong nhân bên trong. Mã SIMD tính toàn bộ hàng C cùng một lúc.

inline void 
    __attribute__ ((gnu_inline))   
    __attribute__ ((aligned(16))) mult2by2B(  
      const double* restrict A, 
      const double* restrict B, 
      double* restrict C 
     ) 

    { 

    register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; 
    xmm0 = _mm_load_pd(C); 
    xmm1 = _mm_load1_pd(A); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 1); 
    xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C,xmm2); 

    xmm0 = _mm_load_pd(C + 2); 
    xmm1 = _mm_load1_pd(A + 2); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 3); 
    //xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C + 2,xmm2); 
} 

Assmebly (Intel Cú pháp)

nỗ lực đầu tiên của tôi là tạo ra một thói quen lắp ráp riêng cho phần này và gọi nó là từ thói quen main. Tuy nhiên, nó rất chậm vì tôi không thể nội tuyến các hàm extern. Tôi đã viết lắp ráp như lắp ráp nội tuyến như hình dưới đây. Nó là giống hệt nhau cho sản phẩm được sản xuất bởi gcc -S -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel. Từ những gì tôi hiểu về sơ đồ kiến ​​trúc vi mô Nehalem, bộ xử lý này có thể thực hiện song song SSE ADD, SSE MUL và SSE MOV, giải thích việc xen kẽ các hướng dẫn MUL, ADD, MOV. Bạn sẽ thấy các hướng dẫn SIMD ở trên theo thứ tự khác vì tôi đã có một sự hiểu biết khác với sách hướng dẫn của Agner Fog. Tuy nhiên, gcc là thông minh và mã SIMD ở trên biên dịch cho hội đồng được hiển thị trong phiên bản nội tuyến.

inline void 
__attribute__ ((gnu_inline))   
__attribute__ ((aligned(16))) mult2by2A 
    ( 
     const double* restrict A, 
     const double* restrict B, 
     double* restrict C 
    ) 
    { 
    __asm__ __volatile__ 
    (
    "mov  edx, %[A]     \n\t" 
    "mov  ecx, %[B]     \n\t" 
    "mov  eax, %[C]     \n\t" 
    "movapd  xmm3, XMMWORD PTR [ecx]  \n\t" 
    "movapd  xmm2, XMMWORD PTR [ecx+16] \n\t" 
    "movddup xmm1, QWORD PTR [edx]  \n\t" 
    "mulpd  xmm1, xmm3     \n\t" 
    "addpd  xmm1, XMMWORD PTR [eax]  \n\t" 
    "movddup xmm0, QWORD PTR [edx+8]  \n\t" 
    "mulpd  xmm0, xmm2     \n\t" 
    "addpd  xmm0, xmm1     \n\t" 
    "movapd  XMMWORD PTR [eax], xmm0  \n\t" 
    "movddup xmm4, QWORD PTR [edx+16] \n\t" 
    "mulpd  xmm4, xmm3     \n\t" 
    "addpd  xmm4, XMMWORD PTR [eax+16] \n\t" 
    "movddup xmm5, QWORD PTR [edx+24] \n\t" 
    "mulpd  xmm5, xmm2     \n\t" 
    "addpd  xmm5, xmm4     \n\t" 
    "movapd  XMMWORD PTR [eax+16], xmm5 \n\t" 
    : // no outputs 
    : // inputs 
    [A] "m" (A), 
    [B] "m" (B), 
    [C] "m" (C) 
    : //register clobber 
    "memory", 
    "edx","ecx","eax", 
    "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5" 
    ); 
} 

Kết quả

tôi biên dịch mã của tôi với những lá cờ sau:

gcc -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel 

Kết quả cho maxiter = 1000000000 dưới:

********** Inline ASM 
L2 norm: 0.000000e+000, Avg. CPU time: 9.563000, Avg. Gflops: 1.673115 

********** SIMD Version 
L2 norm: 0.000000e+000, Avg. CPU time: 0.359000, Avg. Gflops: 44.568245 

Nếu tôi buộc các phiên bản SIMD không được gạch chân với __attribute__ ((noinline)), kết quả là:

********** Inline ASM 
L2 norm: 0.000000e+000, Avg. CPU time: 11.155000, Avg. Gflops: 1.434334 

********** SIMD Version 
L2 norm: 0.000000e+000, Avg. CPU time: 11.264000, Avg. Gflops: 1.420455 

Câu hỏi

1. Nếu cả ASM inline và triển khai SIMD sản xuất sản lượng lắp ráp giống hệt nhau, tại sao là phiên bản lắp ráp nên chậm hơn? Đó là nếu lắp ráp nội tuyến đã không nhận được inlined, được thực hiện rõ ràng bởi các thiết lập thứ hai của kết quả cho thấy hiệu suất giống hệt nhau cho "nội tuyến" ASM so với "noinline" SIMD. Lời giải thích duy nhất tôi có thể tìm thấy là trong Agner Sương mù Tập 2 trang 6:

   mã biên dịch có thể nhanh hơn so với mã lắp ráp bởi vì trình biên dịch có thể làm cho tối ưu hóa liên thủ tục và tối ưu hóa toàn bộ chương trình. Các lập trình viên lắp ráp thường phải thực hiện các chức năng được xác định rõ ràng với giao diện được xác định rõ ràng tuân theo tất cả các quy ước gọi để làm cho mã có thể kiểm tra và có thể xác minh được. Điều này ngăn cản nhiều phương thức tối ưu hóa mà các trình biên dịch sử dụng, chẳng hạn như như chức năng nội tuyến, đăng ký phân bổ, liên tục truyền, phổ biến subexpression loại bỏ qua các chức năng, lập lịch trình trên các chức năng, vv. trong số mã lắp ráp.

   Nhưng đầu ra của bộ kết hợp cho cả hai phiên bản đều giống nhau.

2. Tại sao tôi thấy 44 Gflops trong tập hợp kết quả đầu tiên?Đây là cách trên 12 đỉnh Gflops tôi tính toán, và là những gì tôi mong đợi nếu tôi chạy cả hai lõi với các phép tính chính xác đơn.

EDIT 1 Các bình luận nói có thể có loại bỏ mã chết tôi có thể xác nhận rằng điều này xảy ra đối với các hướng dẫn SIMD. Đầu ra -S cho thấy vòng lặp for cho chỉ số SIMcủa SIMD. Tôi có thể vô hiệu hóa điều đó bằng cách tắt tối ưu hóa trình biên dịch với -O0. Trong trường hợp đó, SIMD chạy 3x chậm như ASM, nhưng ASM vẫn chạy với cùng tốc độ. Định mức cũng không phải là bây giờ, nhưng nó vẫn OK ở 10^-16. Tôi cũng thấy rằng phiên bản ASM nội tuyến đang được inline với các thẻ APP và NO_APP, nhưng nó cũng bị hủy bỏ 8 lần trong vòng lặp for. Tôi nghĩ rằng unrolling nhiều lần sẽ ảnh hưởng đến hiệu suất rất nhiều, như tôi thường unroll vòng 4 lần. Bất cứ điều gì nhiều hơn, theo kinh nghiệm của tôi, dường như làm suy giảm hiệu suất.

Answer 1

GCC là tối ưu hóa đi chức năng nội tuyến của bạn sử dụng intrinsics, mult2by2B, do dòng

C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 

Nếu không có dòng đó phải mất 2,9 giây trên máy tính từ Coliru http://coliru.stacked-crooked.com/a/992304f5f672e257

Và với dòng nó chỉ cần 0.000001 http://coliru.stacked-crooked.com/a/9722c39bb6b8590a

Bạn cũng có thể thấy điều này trong hội đồng. Nếu bạn thả mã bên dưới vào http://gcc.godbolt.org/, bạn sẽ thấy rằng với dòng mã đó, nó bỏ qua toàn bộ hàm.

Tuy nhiên, khi bạn nội tuyến lắp ráp, GCC KHÔNG tối ưu hóa chức năng, mult2by2A, cách xa (mặc dù nó inline nó). Bạn cũng có thể thấy điều này trong hội đồng.

#include <stdio.h> 
#include <emmintrin.h>     // SSE2 
#include <omp.h> 

inline void 
    __attribute__ ((gnu_inline))   
    __attribute__ ((aligned(16))) mult2by2B(  
      const double* __restrict A, 
      const double* __restrict B, 
      double* __restrict C 
     ) 

    { 

    register __m128d xmm0, xmm1, xmm2, xmm3, xmm4; 
    xmm0 = _mm_load_pd(C); 
    xmm1 = _mm_load1_pd(A); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 1); 
    xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C,xmm2); 

    xmm0 = _mm_load_pd(C + 2); 
    xmm1 = _mm_load1_pd(A + 2); 
    xmm2 = _mm_load_pd(B); 
    xmm3 = _mm_load1_pd(A + 3); 
    //xmm4 = _mm_load_pd(B + 2); 
    xmm1 = _mm_mul_pd(xmm1,xmm2); 
    xmm2 = _mm_add_pd(xmm1,xmm0); 
    xmm1 = _mm_mul_pd(xmm3,xmm4); 
    xmm2 = _mm_add_pd(xmm1,xmm2); 
    _mm_store_pd(C + 2,xmm2); 
} 

int main() { 
    double A[4], B[4], C[4]; 
    int maxiter = 10000000; 
    //int maxiter = 1000000000; 
    double dtime; 
    dtime = omp_get_wtime(); 
    for(int i = 0; i < maxiter; i++){ 
     mult2by2B(A,B,C); 
     C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0; 
    } 
    dtime = omp_get_wtime() - dtime; 
    printf("%f %f %f %f\n", C[0], C[1], C[2], C[3]); 
    //gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter; 
    printf("time %f\n", dtime); 
}