respuestas - estadística elemental johnson r

No se puede superar el 50% máximo. rendimiento teórico en matriz multiplicar

Question

Problema

Estoy aprendiendo sobre HPC y optimización de código. Intento replicar los resultados en el papel de multiplicación de la matriz seminal de Goto ( http://www.cs.utexas.edu/users/pingali/CS378/2008sp/papers/gotoPaper.pdf ). A pesar de mis mejores esfuerzos, no puedo obtener más del 50% del rendimiento teórico máximo de CPU.

Fondo

Vea los problemas relacionados aquí ( Multiplicación de matriz optimizada 2x2: ensamblaje lento versus SIMD rápido ), incluida información sobre mi hardware

Lo que he intentado

Este documento relacionado ( http://www.cs.utexas.edu/users/flame/pubs/blis3_ipdps14.pdf ) tiene una buena descripción de la estructura algorítmica de Goto. Proporciono mi código fuente a continuación.

Mi pregunta

Estoy pidiendo ayuda general. He estado trabajando en esto durante demasiado tiempo, he probado muchos algoritmos diferentes, ensamblaje en línea, núcleos internos de varios tamaños (2x2, 4x4, 2x8, ..., mxn con myn grandes), pero parece que no puedo romper 50% CPU Gflops . Esto es puramente para propósitos educativos y no una tarea.

Código fuente

Esperemos que sea comprensible. Por favor, pregunte si no. Configuré la estructura macro (para los bucles) como se describe en el segundo documento anterior. Embalé las matrices como se describe en cualquiera de los documentos y se muestra gráficamente en la Figura 11 aquí ( http://www.cs.utexas.edu/users/flame/pubs/BLISTOMSrev2.pdf ). Mi núcleo interno calcula bloques de 2x8, ya que parece ser el cálculo óptimo para la arquitectura Nehalem (consulte el código fuente de GotoBLAS - núcleos). El núcleo interno se basa en el concepto de calcular las actualizaciones de rango 1 como se describe aquí ( http://code.google.com/p/blis/source/browse/config/template/kernels/3/bli_gemm_opt_mxn.c )

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <x86intrin.h>
#include <math.h>
#include <omp.h>
#include <stdint.h>


// define some prefetch functions
#define PREFETCHNTA(addr,nrOfBytesAhead) /
        _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_NTA)

#define PREFETCHT0(addr,nrOfBytesAhead) /
        _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_T0)

#define PREFETCHT1(addr,nrOfBytesAhead) /
        _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_T1)

#define PREFETCHT2(addr,nrOfBytesAhead) /
        _mm_prefetch(((char *)(addr))+nrOfBytesAhead,_MM_HINT_T2)

// define a min function
#ifndef min
    #define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
#endif

// zero a matrix
void zeromat(double *C, int n)
{
    int i = n;
    while (i--) {
        int j = n;
        while (j--) {
            *(C + i*n + j) = 0.0;
        }
    }
}

// compute a 2x8 block from (2 x kc) x (kc x 8) matrices
inline void 
__attribute__ ((gnu_inline))        
__attribute__ ((aligned(64))) dgemm_2x8_sse(
                int k,
                const double* restrict a1, const int cs_a,
                const double* restrict b1, const int rs_b,
                      double* restrict c11, const int rs_c
                )
{

    register __m128d xmm1, xmm4, //
                    r8, r9, r10, r11, r12, r13, r14, r15; // accumulators

    // 10 registers declared here

    r8 = _mm_xor_pd(r8,r8); // ab
    r9 = _mm_xor_pd(r9,r9);
    r10 = _mm_xor_pd(r10,r10);
    r11 = _mm_xor_pd(r11,r11);

    r12 = _mm_xor_pd(r12,r12); // ab + 8
    r13 = _mm_xor_pd(r13,r13);
    r14 = _mm_xor_pd(r14,r14);
    r15 = _mm_xor_pd(r15,r15);

        // PREFETCHT2(b1,0);
        // PREFETCHT2(b1,64);



    //int l = k;
    while (k--) {

        //PREFETCHT0(a1,0); // fetch 64 bytes from a1

            // i = 0
            xmm1 = _mm_load1_pd(a1);

            xmm4 = _mm_load_pd(b1);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r8 = _mm_add_pd(r8,xmm4);

            xmm4 = _mm_load_pd(b1 + 2);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r9 = _mm_add_pd(r9,xmm4);

            xmm4 = _mm_load_pd(b1 + 4);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r10 = _mm_add_pd(r10,xmm4);

            xmm4 = _mm_load_pd(b1 + 6);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r11 = _mm_add_pd(r11,xmm4);

            //
            // i = 1
            xmm1 = _mm_load1_pd(a1 + 1);

            xmm4 = _mm_load_pd(b1);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r12 = _mm_add_pd(r12,xmm4);

            xmm4 = _mm_load_pd(b1 + 2);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r13 = _mm_add_pd(r13,xmm4);

            xmm4 = _mm_load_pd(b1 + 4);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r14 = _mm_add_pd(r14,xmm4);

            xmm4 = _mm_load_pd(b1 + 6);
            xmm4 = _mm_mul_pd(xmm1,xmm4);
            r15 = _mm_add_pd(r15,xmm4);

        a1 += cs_a;
        b1 += rs_b;

        //PREFETCHT2(b1,0);
        //PREFETCHT2(b1,64);

    }

        // copy result into C

        PREFETCHT0(c11,0);
        xmm1 = _mm_load_pd(c11);
        xmm1 = _mm_add_pd(xmm1,r8);
        _mm_store_pd(c11,xmm1);

        xmm1 = _mm_load_pd(c11 + 2);
        xmm1 = _mm_add_pd(xmm1,r9);
        _mm_store_pd(c11 + 2,xmm1);

        xmm1 = _mm_load_pd(c11 + 4);
        xmm1 = _mm_add_pd(xmm1,r10);
        _mm_store_pd(c11 + 4,xmm1);

        xmm1 = _mm_load_pd(c11 + 6);
        xmm1 = _mm_add_pd(xmm1,r11);
        _mm_store_pd(c11 + 6,xmm1);

        c11 += rs_c;

        PREFETCHT0(c11,0);
        xmm1 = _mm_load_pd(c11);
        xmm1 = _mm_add_pd(xmm1,r12);
        _mm_store_pd(c11,xmm1);

        xmm1 = _mm_load_pd(c11 + 2);
        xmm1 = _mm_add_pd(xmm1,r13);
        _mm_store_pd(c11 + 2,xmm1);

        xmm1 = _mm_load_pd(c11 + 4);
        xmm1 = _mm_add_pd(xmm1,r14);
        _mm_store_pd(c11 + 4,xmm1);

        xmm1 = _mm_load_pd(c11 + 6);
        xmm1 = _mm_add_pd(xmm1,r15);
        _mm_store_pd(c11 + 6,xmm1);

}

// packs a matrix into rows of slivers
inline void 
__attribute__ ((gnu_inline))        
__attribute__ ((aligned(64))) rpack(        double* restrict dst, 
          const double* restrict src, 
            const int kc, const int mc, const int mr, const int n)
{
    double tmp[mc*kc] __attribute__ ((aligned(64)));
    double* restrict ptr = &tmp[0];

    for (int i = 0; i < mc; ++i)
        for (int j = 0; j < kc; ++j)
            *ptr++ = *(src + i*n + j);

    ptr = &tmp[0];

    //const int inc_dst = mr*kc;
    for (int k = 0; k < mc; k+=mr)
        for (int j = 0; j < kc; ++j)
            for (int i = 0; i < mr*kc; i+=kc)
                *dst++ = *(ptr + k*kc + j + i);

}

// packs a matrix into columns of slivers
inline void 
__attribute__ ((gnu_inline))        
__attribute__ ((aligned(64)))  cpack(double* restrict dst, 
                const double* restrict src, 
                const int nc, 
                const int kc, 
                const int nr, 
                const int n)
{
    double tmp[kc*nc] __attribute__ ((aligned(64)));
    double* restrict ptr = &tmp[0];

    for (int i = 0; i < kc; ++i)
        for (int j = 0; j < nc; ++j)
            *ptr++ = *(src + i*n + j);

    ptr = &tmp[0];

    // const int inc_k = nc/nr;
    for (int k = 0; k < nc; k+=nr)
        for (int j = 0; j < kc*nc; j+=nc)
            for (int i = 0; i < nr; ++i)
                *dst++ = *(ptr + k + i + j);

}

void blis_dgemm_ref(
        const int n,
        const double* restrict A,
        const double* restrict B,
        double* restrict C,
        const int mc,
        const int nc,
        const int kc
    )
{
    int mr = 2;
    int nr = 8;
    double locA[mc*kc] __attribute__ ((aligned(64)));
    double locB[kc*nc] __attribute__ ((aligned(64)));
    int ii,jj,kk,i,j;
    #pragma omp parallel num_threads(4) shared(A,B,C) private(ii,jj,kk,i,j,locA,locB)
    {//use all threads in parallel
        #pragma omp for
        // partitions C and B into wide column panels
        for ( jj = 0; jj < n; jj+=nc) {
        // A and the current column of B are partitioned into col and row panels
            for ( kk = 0; kk < n; kk+=kc) {
                cpack(locB, B + kk*n + jj, nc, kc, nr, n);
                // partition current panel of A into blocks
                for ( ii = 0; ii < n; ii+=mc) {
                    rpack(locA, A + ii*n + kk, kc, mc, mr, n);
                    for ( i = 0; i < min(n-ii,mc); i+=mr) {
                        for ( j = 0; j < min(n-jj,nc); j+=nr) {
                            // inner kernel that compues 2 x 8 block
                            dgemm_2x8_sse( kc,
                                       locA + i*kc          ,  mr,
                                       locB + j*kc          ,  nr,
                                       C + (i+ii)*n + (j+jj),  n );
                        }
                    }
                }
            }
        }
    }
}

double compute_gflops(const double time, const int n)
{
    // computes the gigaflops for a square matrix-matrix multiplication
    double gflops;
    gflops = (double) (2.0*n*n*n)/time/1.0e9;
    return(gflops);
}

// ******* MAIN ********//
void main() {
    clock_t time1, time2;
    double time3;
    double gflops;
    const int trials = 10;

    int nmax = 4096;
    printf("%10s %10s/n","N","Gflops/s");

    int mc = 128;
    int kc = 256;
    int nc = 128;

    for (int n = kc; n <= nmax; n+=kc) { //assuming kc is the max dim
        double *A = NULL;
        double *B = NULL;
        double *C = NULL;

        A = _mm_malloc (n*n * sizeof(*A),64);
        B = _mm_malloc (n*n * sizeof(*B),64);
        C = _mm_malloc (n*n * sizeof(*C),64);

        srand(time(NULL));

        // Create the matrices
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                A[i*n + j] = (double) rand()/RAND_MAX;
                B[i*n + j] = (double) rand()/RAND_MAX;
                //D[j*n + i] = B[i*n + j]; // Transpose
                C[i*n + j] = 0.0;
            }
        }

            // warmup
            zeromat(C,n);
            blis_dgemm_ref(n,A,B,C,mc,nc,kc);
            zeromat(C,n);
            time2 = 0;
            for (int count = 0; count < trials; count++){// iterations per experiment here
                    time1 = clock();
                    blis_dgemm_ref(n,A,B,C,mc,nc,kc);
                    time2 += clock() - time1;
                    zeromat(C,n);
                }
            time3 = (double)(time2)/CLOCKS_PER_SEC/trials;
            gflops = compute_gflops(time3, n);
            printf("%10d %10f/n",n,gflops);

        _mm_free(A);
        _mm_free(B);
        _mm_free(C);

        }

    printf("tests are done/n");
}

EDITAR 1

OS = Win 7 64 bit

Compiler = gcc 4.8.1, pero 32 bit y mingw (32 bit también. Estoy trabajando para obtener una versión "no instalable" de mingw64, por lo que puedo generar código / trabajo más rápido con más registros XMM, etc. Si alguien tiene un enlace a una instalación de mingw64 que es similar a la de mingw-get por favor. Mi computadora de trabajo tiene demasiadas restricciones de administrador.

Answer 1

Embalaje

Pareces estar empaquetando el bloque de la matriz A demasiada frecuencia. Tú lo haces

rpack(locA, A + ii*n + kk, kc, mc, mr, n);

Pero esto solo depende de ii y kk y no de jj pero está dentro del bucle interno en jj por lo que reenvasa lo mismo para cada iteración de jj . No creo que sea necesario. En mi código hago el embalaje antes de la multiplicación de matrices. Probablemente es más eficiente empaquetar dentro de la multiplicación de la matriz mientras los valores aún están en el caché, pero es más difícil hacerlo. Pero el empaquetado es una operación O (n ^ 2) y la multiplicación de matrices es una operación O (n ^ 3), por lo que no es muy ineficiente empaquetar fuera de la multiplicación de matrices para matrices grandes (también lo sé por las pruebas, comentando El empaque solo cambia la eficiencia en un pequeño porcentaje. Sin embargo, al volver a rpack con rpack cada iteración jj se ha convertido en una operación O (n ^ 3).

Tiempo de pared

Quieres el tiempo de pared. En Unix, la función clock () no devuelve el tiempo de pared (aunque lo hace en Windows con MSVC). Devuelve el tiempo acumulado para cada hilo. Este es uno de los errores más comunes que he visto en SO para OpenMP.

Use omp_get_wtime() para obtener el tiempo de pared.

Tenga en cuenta que no sé cómo funciona la función clock() con MinGW o MinGW-w64 (son proyectos separados). MinGW se vincula con MSVCRT, por lo que supongo que clock() con MinGW devuelve el tiempo de pared como lo hace con MSVC. Sin embargo, MinGW-w64 no se vincula a MSVCRT (por lo que entiendo que se vincula a algo como glibc). Es posible que clock() en MinGW-w64 realice lo mismo que clock() con Unix.

Hyper Threading

Hyper threading funciona bien para el código que bloquea la CPU a menudo. Esa es la mayoría del código porque es muy difícil escribir código que no detenga la CPU. Es por eso que Intel inventó Hyper Threading. Es más fácil cambiar de tarea y darle a la CPU algo más que hacer que optimizar el código. Sin embargo, para el código que está altamente optimizado, el hyper-threading puede dar peores resultados. En mi propio código de multiplicación de matrices es ciertamente el caso. Establezca la cantidad de subprocesos en la cantidad de núcleos físicos que tiene (dos en su caso).

Mi código

A continuación se muestra mi código. No inner64 función inner64 aquí. Puede encontrarlo en Diferencia en el rendimiento entre MSVC y GCC para obtener un código de multplicación matricial altamente optimizado (con el nombre desagradable y engañoso de AddDot4x4_vec_block_8wide )

Escribí este código antes de leer el documento de Goto y también antes de leer los manuales de optimización de Agner Fog. Pareces reordenar / empaquetar las matrices en el bucle principal. Eso probablemente tenga más sentido. No creo que los reordene de la misma manera que usted y también reordené solo una de las matrices de entrada (B) y no tanto como usted.

El rendimiento de este código en mi sistema (Xeon [email protected]) con Linux y GCC es aproximadamente el 75% del pico para este tamaño de matriz (4096x4096). El MKL de Intel recibe aproximadamente el 94% del pico en mi sistema para este tamaño de matriz, por lo que claramente hay espacio para mejorar.

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <immintrin.h>

extern "C" void inner64(const float *a, const float *b, float *c);
void (*fp)(const float *a, const float *b, float *c) = inner64;

void reorder(float * __restrict a, float * __restrict b, int n, int bs) {
    int nb = n/bs;
    #pragma omp parallel for
    for(int i=0; i<nb; i++) {
        for(int j=0; j<nb; j++) {
            for(int i2=0; i2<bs; i2++) {
                for(int j2=0; j2<bs; j2++) {
                    b[bs*bs*(nb*i+j) + bs*i2+j2]= a[bs*(i*n+j) + i2*n + j2];    
                }
            }
        }
    }
}

inline void gemm_block(float * __restrict a, float * __restrict b, float * __restrict c, int n, int n2) {
    for(int i=0; i<n2; i++) {
        fp(&a[i*n], b, &c[i*n]);
    }
}

void gemm(float * __restrict a, float * __restrict b, float * __restrict c, int n, int bs) {
    int nb = n/bs;
    float *b2 = (float*)_mm_malloc(sizeof(float)*n*n,64);
    reorder(b,b2,n,bs);
    #pragma omp parallel for
    for(int i=0; i<nb; i++) {
        for(int j=0; j<nb; j++) {
            for(int k=0; k<nb; k++) {
                gemm_block(&a[bs*(i*n+k)],&b2[bs*bs*(k*nb+j)],&c[bs*(i*n+j)], n, bs);
            }
        }
    }
    _mm_free(b2);
}

int main() {
    float peak = 1.0f*8*4*2*3.69f;
    const int n = 4096;
    float flop = 2.0f*n*n*n*1E-9f;
    omp_set_num_threads(4);

    float *a = (float*)_mm_malloc(sizeof(float)*n*n,64);
    float *b = (float*)_mm_malloc(sizeof(float)*n*n,64);
    float *c = (float*)_mm_malloc(sizeof(float)*n*n,64);
    for(int i=0; i<n*n; i++) {
        a[i] = 1.0f*rand()/RAND_MAX;
        b[i] = 1.0f*rand()/RAND_MAX;
    }

    gemm(a,b,c,n,64); //warm OpenMP up
    while(1) {
        for(int i=0; i<n*n; i++) c[i] = 0;
        double dtime = omp_get_wtime();
        gemm(a,b,c,n,64);   
        dtime = omp_get_wtime() - dtime;
        printf("time %.2f s, efficiency %.2f%%/n", dtime, 100*flop/dtime/peak);
    }
}