OpenMP並列化（Block Matrix Mult）

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はまだ組み込み関数が必要です）。次のコードは、4コア/ 8ハードウェアスレッドSkylakeシステムのピークFLOPSの60％を取得します。ブロック行列の乗算を使用します。

ハイパースレッディングはパフォーマンスが低下しますので、スレッドの移行を防ぐために、コア数とスレッド数を同じにしてください。

export OMP_PROC_BIND=true 
export OMP_NUM_THREADS=4 

そしてこれは無限ループの反復ごとGEMM 10回行い、低値、中央値、およびプの高い比を出力し、この

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); 
    } 
} 

をコンパイルpeak_FLOPSまで、そして最後に中央値FLOPSまで。

（有効な場合ターボで）あなたは物理コアの数は、すべてのコアの周波数に次の行を

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

を調整する必要があります

、及びため16でコアあたりの浮動ポインタ操作の数Core2-Ivy Bridge、Haswell-Kaby Lakeの場合は32、Xeon Phi Knights Landingの場合は64です。

このコードは、NUMAシステムでは効率が悪い場合があります。それはKnight Landingとほとんど同じではありません（私はこれを調べ始めました）。