OpenMP Parallelisierung (Block Matrix Mult)

Question

Ich versuche, Block Matrix Multiplikation zu implementieren und es parallelisierter zu machen.

Dies ist mein Code:

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

Ich habe um mit OpenMP spielen aber noch kein Glück gehabt in herauszufinden, was der beste Weg, dies in der geringsten Menge an Zeit getan zu haben.

Answer 1

Gute Leistung durch Matrixmultiplikation ist eine große Aufgabe. Da "der beste Code der Code ist, den ich nicht schreiben muss", wäre es viel besser, die Zeit zu nutzen, um zu verstehen, wie man eine BLAS-Bibliothek benutzt.

Wenn Sie X86-Prozessoren verwenden, ist die Intel Math Kernel Library (MKL) kostenlos verfügbar und enthält optimierte, parallelisierte Matrixmultiplikationsoperationen. https://software.intel.com/en-us/articles/free-mkl

(FWIW, ich für Intel arbeiten, aber nicht auf MKL :-))

Answer 2

Ich suche in dichte Matrixmultiplikation vor kurzem begonnen (GEMM) wieder. Es stellt sich heraus, dass der Clang-Compiler sehr gut in der Optimierung von GEMM ist, ohne irgendwelche intrinsischen Eigenschaften zu benötigen (GCC benötigt immer noch intrinsische Eigenschaften). Der folgende Code bekommt 60% der Peak FLOPS meines Skylake Systems mit vier Cores und acht Hardware Threads. Es verwendet die Blockmatrix-Multiplikation.

Hyper-Threading führt zu schlechterer Leistung. Daher stellen Sie sicher, dass Sie nur Threads verwenden, die der Anzahl der Kerne entsprechen, und binden Sie Threads, um eine Thread-Migration zu verhindern.

export OMP_PROC_BIND=true 
export OMP_NUM_THREADS=4 

kompilieren dann so

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

Der Code

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

Dies gilt GEMM 10 Mal pro Iteration von einer Endlos-Schleife und druckt die niedrige, mittlere und hohe Verhältnis von FLOPS zu peak_FLOPS und schließlich den medianen FLOPS.

Sie müssen die folgenden Zeilen

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

auf die Anzahl der physikalischen Kerne, die Frequenz für alle Kerne (mit Turbo falls aktiviert) einzustellen, und die Anzahl der Zeigeroperationen pro Kern floating welches 16 für Core2-Ivy Bridge, 32 für Haswell-Kaby Lake und 64 für die Xeon Phi Knights Landing.

Dieser Code ist bei NUMA-Systemen möglicherweise weniger effizient. Es geht nicht annähernd so gut mit Knight Landing (ich habe gerade angefangen, dies zu untersuchen).