OpenMP並行化（塊矩陣多）

Question

我正在嘗試實現塊矩陣乘法並使其更加並行化。

這是我的代碼：

int i,j,jj,k,kk; 
float sum; 
int en = 4 * (2048/4); 
    #pragma omp parallel for collapse(2) 
for(i=0;i<2048;i++) { 
    for(j=0;j<2048;j++) { 
     C[i][j]=0; 
    } 
} 
for (kk=0;kk<en;kk+=4) { 
    for(jj=0;jj<en;jj+=4) { 
     for(i=0;i<2048;i++) { 
      for(j=jj;j<jj+4;j++) { 
       sum = C[i][j]; 
       for(k=kk;k<kk+4;k++) { 
        sum+=A[i][k]*B[k][j]; 
       } 
       C[i][j] = sum; 
      } 
     } 
    } 
} 

我一直在玩弄OpenMP的，但仍然有在盤算什麼有最好的方式在最短的時間內量完成這件事沒有運氣。

Answer 1

從矩陣乘法中獲得良好的性能是一項很大的工作。由於「最好的代碼是我不必寫的代碼」，更好地利用你的時間將是瞭解如何使用BLAS庫。

如果您使用的是X86處理器，則可以免費獲得英特爾數學核心函數庫（MKL），幷包含優化的並行化矩陣乘法運算。 https://software.intel.com/en-us/articles/free-mkl

（FWIW，我爲英特爾工作，但不是在MKL :-)）

Answer 2

最近，我開始尋找到稠密矩陣乘法（GEMM）一次。事實證明，Clang編譯器非常擅長優化GEMM，而不需要任何內部函數（GCC仍然需要內部函數）。以下代碼獲得了我的四核/八硬件線程Skylake系統的峯值FLOPS的60％。它使用塊矩陣乘法。

超線程的性能更差，因此您確保只使用等於核心數的線程並綁定線程以防止線程遷移。

export OMP_PROC_BIND=true 
export OMP_NUM_THREADS=4 

然後編譯這樣

clang -Ofast -march=native -fopenmp -Wall gemm_so.c 

代碼

#include <stdlib.h> 
#include <string.h> 
#include <stdio.h> 
#include <omp.h> 
#include <x86intrin.h> 

#define SM 80 

typedef __attribute((aligned(64))) float * restrict fast_float; 

static void reorder2(fast_float a, fast_float b, int n) { 
    for(int i=0; i<SM; i++) memcpy(&b[i*SM], &a[i*n], sizeof(float)*SM); 
} 

static void kernel(fast_float a, fast_float b, fast_float c, int n) { 
    for(int i=0; i<SM; i++) { 
    for(int k=0; k<SM; k++) { 
     for(int j=0; j<SM; j++) { 
     c[i*n + j] += a[i*n + k]*b[k*SM + j]; 
     } 
    } 
    } 
} 

void gemm(fast_float a, fast_float b, fast_float c, int n) { 
    int bk = n/SM; 

    #pragma omp parallel 
    { 
    float *b2 = _mm_malloc(sizeof(float)*SM*SM, 64); 
    #pragma omp for collapse(3) 
    for(int i=0; i<bk; i++) { 
     for(int j=0; j<bk; j++) { 
     for(int k=0; k<bk; k++) { 
      reorder2(&b[SM*(k*n + j)], b2, n); 
      kernel(&a[SM*(i*n+k)], b2, &c[SM*(i*n+j)], n); 
     } 
     } 
    } 
    _mm_free(b2); 
    } 
} 

static int doublecmp(const void *x, const void *y) { return *(double*)x < *(double*)y ? -1 : *(double*)x > *(double*)y; } 

double median(double *x, int n) { 
    qsort(x, n, sizeof(double), doublecmp); 
    return 0.5f*(x[n/2] + x[(n-1)/2]); 
} 

int main(void) { 
    int cores = 4; 
    double frequency = 3.1; // i7-6700HQ turbo 4 cores 
    double peak = 32*cores*frequency; 

    int n = SM*10*2; 

    int mem = sizeof(float) * n * n; 
    float *a = _mm_malloc(mem, 64); 
    float *b = _mm_malloc(mem, 64); 
    float *c = _mm_malloc(mem, 64); 

    memset(a, 1, mem), memset(b, 1, mem); 

    printf("%dx%d matrix\n", n, n); 
    printf("memory of matrices: %.2f MB\n", 3.0*mem*1E-6); 
    printf("peak SP GFLOPS %.2f\n", peak); 
    puts(""); 

    while(1) { 
    int r = 10; 
    double times[r]; 
    for(int j=0; j<r; j++) { 
     times[j] = -omp_get_wtime(); 
     gemm(a, b, c, n); 
     times[j] += omp_get_wtime(); 
    } 

    double flop = 2.0*1E-9*n*n*n; //GFLOP 
    double time_mid = median(times, r); 
    double flops_low = flop/times[r-1], flops_mid = flop/time_mid, flops_high = flop/times[0]; 
    printf("%.2f %.2f %.2f %.2f\n", 100*flops_low/peak, 100*flops_mid/peak, 100*flops_high/peak, flops_high); 
    } 
} 

這確實GEMM每一個無限循環的迭代10次並打印低，中值和FLOPS的高比到peak_FLOPS，最後是中值FLOPS。

您將需要調整以下行

int cores = 4; 
double frequency = 3.1; // i7-6700HQ turbo 4 cores 
double peak = 32*cores*frequency; 

物理內核的所有核心數量，頻率（與如果啓用渦輪增壓）和浮動每個核心的指針操作的數量是16爲Core2-Ivy Bridge，32爲Haswell-Kaby Lake，64爲Xeon Phi Knights Landing。

對於NUMA系統，此代碼可能效率較低。這跟Knight Landing差不多（我剛開始研究這個）。