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差不多（我刚开始研究这个）。