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 系统的 60% 的峰值 FLOPS。它使用块矩阵乘法。

超线程的性能更差，因此请确保只使用与内核数量相等的线程并绑定线程以防止线程迁移。

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;

物理核心的数量，所有核心的频率（如果启用涡轮增压），以及每个核心的浮动指针操作数，Core2-Ivy Bridge 为16，Haswell-Kaby Lake 为32，和64 Xeon Phi Knights Landing。

此代码在 NUMA 系统中可能效率较低。与 Knight Landing 相比，它的效果几乎没有（我刚开始研究这个）。