c++ - Cómo optimizar el código de multiplicación matricial(matmul) para que se ejecute rápidamente en un solo núcleo de procesador

optimization parallel-processing

Question

Estoy trabajando en conceptos de programación paralela y tratando de optimizar el ejemplo de multiplicación de matrices en un solo núcleo. La implementación más rápida que se me ocurrió hasta ahora es la siguiente:

/* This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values. */    
void square_dgemm (int n, double* A, double* B, double* C)
{
  /* For each row i of A */
  for (int i = 0; i < n; ++i)
    /* For each column j of B */
    for (int j = 0; j < n; ++j) 
    {
      /* Compute C(i,j) */
      double cij = C[i+j*n];
      for( int k = 0; k < n; k++ )
          cij += A[i+k*n] * B[k+j*n];
      C[i+j*n] = cij;
    }
}

Los resultados son como a continuación. Cómo reducir los bucles y aumentar el rendimiento

login4.stampede(72)$ tail -f job-naive.stdout
Size: 480       Mflop/s:  1818.89       Percentage: 18.95
Size: 511       Mflop/s:  2291.73       Percentage: 23.87
Size: 512       Mflop/s:  937.061       Percentage:  9.76
Size: 639       Mflop/s:  293.434       Percentage:  3.06
Size: 640       Mflop/s:  270.238       Percentage:  2.81
Size: 767       Mflop/s:  240.209       Percentage:  2.50
Size: 768       Mflop/s:  242.118       Percentage:  2.52
Size: 769       Mflop/s:  240.173       Percentage:  2.50
Average percentage of Peak = 22.0802
Grade = 33.1204

Answer 1

En lugar de optimizar, puede ofuscar el código para que parezca que está optimizado.

Aquí hay una multiplicación de matriz con un solo bucle for cuerpo nulo (!):

/* This routine performs a dgemm operation
 *  C := C + A * B
 * where A, B, and C are lda-by-lda matrices stored in column-major format.
 * On exit, A and B maintain their input values. 
 * This implementation uses a single for loop: it has been optimised for space,
 * namely vertical space in the source file! */    
void square_dgemm(int n, const double *A, const double *B, double *C) {
    for (int i = 0, j = 0, k = -1;
         ++k < n || ++j < n + (k = 0) || ++i < n + (j = 0);
         C[i+j*n] += A[i+k*n] * B[k+j*n]) {}
}

Answer 2

Hay muchas formas de mejoras directas. La optimización básica es lo que escribió Rick James. Además, puede reorganizar la primera matriz por filas y la segunda por columnas. Luego, en tus bucles for () siempre harás ++ y nunca harás + = n. Los bucles donde saltas por n son mucho más lentos en comparación con ++.

Pero la mayoría de esas optimizaciones aguantan el golpe porque un buen compilador las hará por usted cuando use las banderas -O3 u -O4. Desenrollará los bucles, reutilizará registros, realizará operaciones lógicas en lugar de multiplicaciones, etc. Incluso cambiará el orden de sus bucles for i y for j si es necesario.

El problema principal con su código es que cuando tiene matrices NxN, utiliza 3 bucles que lo obligan a realizar operaciones O(N^3) . Esto es muy lento Creo que los algoritmos de última generación solo realizan operaciones ~ O(N^2.37) ( enlace aquí ). Para matrices grandes (digamos N = 5000), esta es una excelente optimización. Puede implementar el algoritmo Strassen fácilmente, lo que le dará una mejora de ~ N ^ 2.87 o usarlo en combinación con el algoritmo Karatsuba , que puede acelerar las cosas incluso para optimizaciones escalares regulares. No implemente nada por su cuenta. Descargue una implementación de código abierto. Multiplicar matrices como un gran tema con mucha investigación y algoritmos muy rápidos. El uso de 3 bucles no se considera una forma válida de hacer este trabajo de manera eficiente. Buena suerte

Answer 3

La implementación de vanguardia de la multiplicación de matrices en las CPU utiliza el algoritmo GotoBLAS . Básicamente, los bucles están organizados en el siguiente orden:

Loop5 for jc = 0 to N-1 in steps of NC
Loop4   for kc = 0 to K-1 in steps of KC
          //Pack KCxNC block of B
Loop3     for ic = 0 to M-1 in steps of MC
            //Pack MCxKC block of A
//--------------------Macro Kernel------------
Loop2       for jr = 0 to NC-1 in steps of NR
Loop1         for ir = 0 to MC-1 in steps of MR
//--------------------Micro Kernel------------
Loop0           for k = 0 to KC-1 in steps of 1
                //update MRxNR block of C matrix

Una idea clave que subyace a las implementaciones modernas de alto rendimiento de la multiplicación de matrices es organizar los cálculos dividiendo los operandos en bloques para la localidad temporal (3 bucles más externos), y empacar (copiar) dichos bloques en búferes contiguos que encajan en varios niveles de memoria para localidad espacial (3 bucles más internos).

La figura anterior (originalmente de este documento , utilizada directamente en este tutorial ) ilustra el algoritmo GotoBLAS implementado en BLIS . Los parámetros de bloqueo de caché {MC, NC, KC} determinan los tamaños de las submatrices de Bp (KC × NC) y Ai (MC × KC), de modo que encajan en varios cachés. Durante el cálculo, los paneles de filas Bp se empaquetan contiguamente en el búfer Bp para caber en el caché L3. Los bloques Ai se empaquetan de manera similar en el búfer Ai para caber en el caché L2. Los tamaños de bloque de registro {MR, NR} se relacionan con submatrices en registros que contribuyen a C. En el micro-núcleo (el bucle más interno), un pequeño micro-mosaico MR × NR de C se actualiza mediante un par de MR × KC y KC × NR astillas de Ai y Bp.

Para el algoritmo de Strassen con complejidad O (N ^ 2.87), es posible que le interese leer este documento . Otros algoritmos rápidos de multiplicación de matrices con complejidad asintótica menor que O (N ^ 3) pueden extenderse fácilmente en este documento . Existe una tesis reciente sobre los prácticos algoritmos de multiplicación rápida de matrices.

Los siguientes tutoriales pueden ser útiles si desea obtener más información sobre cómo optimizar la multiplicación de matrices en las CPU:

Cómo optimizar Wiki GEMM

GEMM: de Pure C a SSE Micro Kernels optimizados

BLISlab: un sandbox para optimizar GEMM para CPU y ARM

El documento más actualizado sobre cómo optimizar GEMM en las CPU (con AVX2 / FMA) paso a paso se puede descargar aquí: https://github.com/ULAFF/LAFF-On-HPC/blob/master/LAFF-On-PfHP.pdf

Un curso masivo abierto en línea que se ofrecerá en edX a partir de junio de 2019 (programación LAFF-On para alto rendimiento): https://github.com/ULAFF/LAFF-On-HPC http://www.cs.utexas.edu/users/flame/laff/pfhp/LAFF-On-PfHP.html

Answer 4

Mi C está bastante oxidado, y no sé qué de lo siguiente ya está haciendo el optimizador, pero aquí va ...

Como prácticamente todo el tiempo se dedica a hacer un producto de puntos, permítanme optimizarlo; puedes construir desde allí.

double* pa = &A[i];
double* pb = &B[j*n];
double* pc = &C[i+j*n];
for( int k = 0; k < n; k++ )
{
    *pc += *pa++ * *pb;
    pb += n;
}

Es probable que su código pase más tiempo en aritmética de subíndices que cualquier otra cosa. Mi código usa +=8 y +=(n<<3) , que es mucho más eficiente. (Nota: un double toma 8 bytes).

Otras optimizaciones:

Si conoce el valor de n , podría "desenrollar" al menos el bucle más interno. Esto elimina la sobrecarga del for .

Incluso si supiera que n es par, puede repetir n / 2 veces, duplicando el código en cada iteración. Esto reduciría la sobrecarga a la mitad (aprox.).

No verifiqué si la multiplicación de la matriz podría hacerse mejor en el orden de fila mayor versus columna mayor. +=8 es más rápido que +=(n<<3) ; Esto sería una pequeña mejora en los bucles externos.

Otra forma de "desenrollar" sería hacer dos productos punto en el mismo bucle interno. (Supongo que me estoy volviendo demasiado complejo como para explicarlo).

Las CPU son "hiperescalares" en estos días. Esto significa que pueden, hasta cierto punto, hacer varias cosas al mismo tiempo. Pero eso no significa que las cosas que deben hacerse consecutivamente puedan optimizarse de esa manera. Hacer dos productos punto independientes en el mismo bucle puede proporcionar más oportunidades para el hiperscaling.