por - multiplicacion de matrices online

Replicando el rendimiento de la multiplicación de matrices BLAS: ¿Puedo hacerlo coincidir? (1)

Fondo

Si ha estado siguiendo mis publicaciones, estoy intentando replicar los resultados encontrados en el artículo seminal de Kazushige Goto sobre la multiplicación de la matriz cuadrada C = AB . Mi último post sobre este tema se puede encontrar here . En esa versión de mi código, sigo la estrategia de empaquetado y almacenamiento de memoria de Goto con un kernel interno que computa bloques 2x8 de C usando intrínsecos SSE3 de 128 bits. Mi CPU es i5-540M con hyperthreading apagado. Información adicional sobre mi hardware se puede encontrar en otra post y se repite a continuación.

Mi hardware

Mi CPU es un Intel i5 - 540M. Puede encontrar la información relevante de CPUID en cpu-world.com. La microarquitectura es Nehalem (westmere), por lo que teóricamente puede calcular 4 fracasos de doble precisión por núcleo por ciclo. Estaré usando solo un núcleo (no OpenMP), por lo tanto, con el subprocesamiento desactivado y el Intel Turbo Boost de 4 pasos, debería ver un máximo de ( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops . Para referencia, con ambos núcleos funcionando al máximo, Intel Turbo Boost da una aceleración de 2 pasos y debería obtener un máximo teórico de 22.4 DP Gflops .

Mi software

Windows7 de 64 bits, pero MinGW / GCC de 32 bits debido a restricciones en mi computadora.

¿Qué hay de nuevo esta vez?

Yo 2x4 bloques de 2x4 de C Esto ofrece un mejor rendimiento y está en línea con lo que dice Goto (la mitad de los registros deben usarse para calcular C ). He probado muchos tamaños: 1x8 , 2x8 , 2x4 , 4x2 , 2x2 , 4x4 .
Mi núcleo interno es un ensamblado x86 codificado a mano, optimizado a mi mejor capacidad (coincide con algunos de los núcleos que escribió Goto), lo que proporciona un aumento de rendimiento bastante grande en SIMD. Este código se desenrolla 8 veces dentro del kernel interno (definido como una macro para mayor comodidad), lo que brinda el mejor rendimiento de las otras estrategias de enrollamiento que probé.
Uso los contadores de rendimiento de Windows para clock() los códigos, en lugar de clock()
Programo el núcleo interno independientemente del código total, para ver qué tan bien está funcionando mi ensamblaje codificado a mano.
Reporto el mejor resultado de una serie de ensayos, en lugar de un promedio de los ensayos.
No más OpenMP (solo rendimiento de un solo núcleo).
NOTA Recopilaré OpenBLAS esta noche para usar solo 1 núcleo para poder comparar.

Algunos resultados preliminares

N es la dimensión de las matrices cuadradas, Total Gflops/s es el Gflops / s de todo el código, y Kernel Gflops/s es el Gflops / s del núcleo interno. Puede ver que con un pico de 12.26 Gflops / s en un núcleo, el núcleo interno está obteniendo alrededor de un 75% de eficiencia, mientras que el código general es aproximadamente un 60% de eficiencia.

Me gustaría acercarme al 95% de eficiencia para el núcleo y al 80% para el código general. ¿Qué más puedo hacer para mejorar el rendimiento de, al menos, el núcleo interno?

N Total Gflops/s Kernel Gflops/s 256 7.778089 9.756284 512 7.308523 9.462700 768 7.223283 9.253639 1024 7.197375 9.132235 1280 7.142538 8.974122 1536 7.114665 8.967249 1792 7.060789 8.861958

Código fuente

Si te sientes particularmente magnánimo, prueba mi código en tu máquina. Compilado con gcc -std=c99 -O2 -m32 -msse3 -mincoming-stack-boundary=2 -masm=intel somatmul2.c -o somatmul2.exe . Siéntase libre de probar otras banderas, pero he encontrado que estas funcionan mejor en mi máquina.

#include <stdio.h> #include <time.h> #include <stdlib.h> #include <string.h> #include <xmmintrin.h> #include <math.h> #include <omp.h> #include <stdint.h> #include <windows.h> #ifndef min #define min( a, b ) ( ((a) < (b)) ? (a) : (b) ) #endif 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+=4) { *ptr++ = *(src + i*n + j ); *ptr++ = *(src + i*n + j+1); *ptr++ = *(src + i*n + j+2); *ptr++ = *(src + i*n + j+3); } 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); } 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) { //nc cols, kc rows, nr size ofregister strips 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+=4) { *ptr++ = *(src + i*n + j ); *ptr++ = *(src + i*n + j+1); *ptr++ = *(src + i*n + j+2); *ptr++ = *(src + i*n + j+3); } 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); } #define KERNEL0(add0,add1,add2,add3) / "mulpd xmm4, xmm6 /n/t" / "addpd xmm0, xmm4 /n/t" / "movapd xmm4, XMMWORD PTR [edx+"#add2"] /n/t" / "mulpd xmm7, xmm4 /n/t" / "addpd xmm1, xmm7 /n/t" / "movddup xmm5, QWORD PTR [eax+"#add0"] /n/t" / "mulpd xmm6, xmm5 /n/t" / "addpd xmm2, xmm6 /n/t" / "movddup xmm7, QWORD PTR [eax+"#add1"] /n/t" / "mulpd xmm4, xmm5 /n/t" / "movapd xmm6, XMMWORD PTR [edx+"#add3"] /n/t" / "addpd xmm3, xmm4 /n/t" / "movapd xmm4, xmm7 /n/t" / " /n/t" inline void __attribute__ ((gnu_inline)) __attribute__ ((aligned(64))) dgemm_2x4_asm_j ( const int mc, const int nc, const int kc, const double* restrict locA, const int cs_a, // mr const double* restrict locB, const int rs_b, // nr double* restrict C, const int rs_c ) { const double* restrict a1 = locA; for (int i = 0; i < mc ; i+=cs_a) { const double* restrict b1 = locB; double* restrict c11 = C + i*rs_c; for (int j = 0; j < nc ; j+=rs_b) { __asm__ __volatile__ ( "mov eax, %[a1] /n/t" "mov edx, %[b1] /n/t" "mov edi, %[c11] /n/t" "mov ecx, %[kc] /n/t" "pxor xmm0, xmm0 /n/t" "movddup xmm7, QWORD PTR [eax] /n/t" // a1 "pxor xmm1, xmm1 /n/t" "movapd xmm6, XMMWORD PTR [edx] /n/t" // b1 "pxor xmm2, xmm2 /n/t" "movapd xmm4, xmm7 /n/t" // a1 "pxor xmm3, xmm3 /n/t" "sar ecx, 3 /n/t" // divide by 2^num "L%=: /n/t" // start loop KERNEL0( 8, 16, 16, 32) KERNEL0( 24, 32, 48, 64) KERNEL0( 40, 48, 80, 96) KERNEL0( 56, 64, 112, 128) KERNEL0( 72, 80, 144, 160) KERNEL0( 88, 96, 176, 192) KERNEL0( 104, 112, 208, 224) KERNEL0( 120, 128, 240, 256) "add eax, 128 /n/t" "add edx, 256 /n/t" " /n/t" "dec ecx /n/t" "jne L%= /n/t" // end loop " /n/t" "mov esi, %[rs_c] /n/t" // don''t need cs_a anymore "sal esi, 3 /n/t" // times 8 "lea ebx, [edi+esi] /n/t" // don''t need b1 anymore "addpd xmm0, XMMWORD PTR [edi] /n/t" // c11 "addpd xmm1, XMMWORD PTR [edi+16] /n/t" // c11 + 2 "addpd xmm2, XMMWORD PTR [ebx] /n/t" // c11 "addpd xmm3, XMMWORD PTR [ebx+16] /n/t" // c11 + 2 "movapd XMMWORD PTR [edi], xmm0 /n/t" "movapd XMMWORD PTR [edi+16], xmm1 /n/t" "movapd XMMWORD PTR [ebx], xmm2 /n/t" "movapd XMMWORD PTR [ebx+16], xmm3 /n/t" : // no outputs : // inputs [kc] "m" (kc), [a1] "m" (a1), [cs_a] "m" (cs_a), [b1] "m" (b1), [rs_b] "m" (rs_b), [c11] "m" (c11), [rs_c] "m" (rs_c) : //register clobber "memory", "eax","ebx","ecx","edx","esi","edi", "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5","xmm6","xmm7" ); b1 += rs_b*kc; c11 += rs_b; } a1 += cs_a*kc; } } double 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 = 4; double locA[mc*kc] __attribute__ ((aligned(64))); double locB[kc*nc] __attribute__ ((aligned(64))); LARGE_INTEGER frequency, time1, time2; double time3 = 0.0; QueryPerformanceFrequency(&frequency); // zero C memset(C, 0, n*n*sizeof(double)); int ii,jj,kk; //#pragma omp parallel num_threads(2) shared(A,B,C) private(ii,jj,kk,locA,locB) {//use all threads in parallel //#pragma omp for for ( jj = 0; jj < n; jj+=nc) { for ( kk = 0; kk < n; kk+=kc) { cpack(locB, B + kk*n + jj, nc, kc, nr, n); for ( ii = 0; ii < n; ii+=mc) { rpack(locA, A + ii*n + kk, kc, mc, mr, n); QueryPerformanceCounter(&time1); dgemm_2x4_asm_j( mc, nc, kc, locA , mr, locB , nr, C + ii*n + jj, n ); QueryPerformanceCounter(&time2); time3 += (double) (time2.QuadPart - time1.QuadPart); } } } } return time3 / frequency.QuadPart; } 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); } void main() { LARGE_INTEGER frequency, time1, time2; double time3, best_time; double kernel_time, best_kernel_time; QueryPerformanceFrequency(&frequency); int best_flag; double gflops, kernel_gflops; const int trials = 100; int nmax = 4096; printf("%16s %16s %16s/n","N","Total Gflops/s","Kernel Gflops/s"); int mc = 256; int kc = 256; int nc = 128; for (int n = kc; n <= nmax; n+=kc) { double *A = NULL; double *B = NULL; double *C = NULL; A = _mm_malloc (n*n * sizeof(*A),64); if (!A) {printf("A failed/n"); return;} B = _mm_malloc (n*n * sizeof(*B),64); if (!B) {printf("B failed/n"); return;} C = _mm_malloc (n*n * sizeof(*C),64); if (!C) {printf("C failed/n"); return;} 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; } } // warmup blis_dgemm_ref(n,A,B,C,mc,nc,kc); for (int count = 0; count < trials; count++){ QueryPerformanceCounter(&time1); kernel_time = blis_dgemm_ref(n,A,B,C,mc,nc,kc); QueryPerformanceCounter(&time2); time3 = (double)(time2.QuadPart - time1.QuadPart) / frequency.QuadPart; if (count == 0) { best_time = time3; best_kernel_time = kernel_time; } else { best_flag = ( time3 < best_time ? 1 : 0 ); if (best_flag) { best_time = time3; best_kernel_time = kernel_time; } } } gflops = compute_gflops(best_time, n); kernel_gflops = compute_gflops(best_kernel_time, n); printf("%16d %16f %16f/n",n,gflops,kernel_gflops); _mm_free(A); _mm_free(B); _mm_free(C); } printf("tests are done/n"); }

EDITAR

Reemplace las funciones de embalaje con las siguientes versiones nuevas y más rápidas:

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) { for (int i = 0; i < mc/mr; ++i) for (int j = 0; j < kc; ++j) for (int k = 0; k < mr; ++k) *dst++ = *(src + i*n*mr + k*n + j); } 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) { for (int i = 0; i < nc/nr; ++i) for (int j = 0; j < kc; ++j) for (int k = 0; k < nr; ++k) *dst++ = *(src + i*nr + j*n + k); }

Resultados con nuevas funciones de embalaje

Buen impulso para el rendimiento general:

N Total Gflops/s Kernel Gflops/s 256 7.915617 8.794849 512 8.466467 9.350920 768 8.354890 9.135575 1024 8.168944 8.884611 1280 8.174249 8.825920 1536 8.285458 8.938712 1792 7.988038 8.581001

LINPACK de 32 bits con 1 hilo

CPU frequency: 2.792 GHz Number of CPUs: 1 Number of cores: 2 Number of threads: 1 Performance Summary (GFlops) Size LDA Align. Average Maximal 128 128 4 4.7488 5.0094 256 256 4 6.0747 6.9652 384 384 4 6.5208 7.2767 512 512 4 6.8329 7.5706 640 640 4 7.4278 7.8835 768 768 4 7.7622 8.0677 896 896 4 7.8860 8.4737 1024 1024 4 7.7292 8.1076 1152 1152 4 8.0411 8.4738 1280 1280 4 8.1429 8.4863 1408 1408 4 8.2284 8.7073 1536 1536 4 8.3753 8.6437 1664 1664 4 8.6993 8.9108 1792 1792 4 8.7576 8.9176 1920 1920 4 8.7945 9.0678 2048 2048 4 8.5490 8.8827 2176 2176 4 9.0138 9.1161 2304 2304 4 8.1402 9.1446 2432 2432 4 9.0003 9.2082 2560 2560 4 8.8560 9.1197 2688 2688 4 9.1008 9.3144 2816 2816 4 9.0876 9.3089 2944 2944 4 9.0771 9.4191 3072 3072 4 8.9402 9.2920 3200 3200 4 9.2259 9.3699 3328 3328 4 9.1224 9.3821 3456 3456 4 9.1354 9.4082 3584 3584 4 9.0489 9.3351 3712 3712 4 9.3093 9.5108 3840 3840 4 9.3307 9.5324 3968 3968 4 9.3895 9.5352 4096 4096 4 9.3269 9.3872

Editar 2

Aquí hay algunos resultados de OpenBLAS de un solo hilo, que tardó 4 horas en compilarse anoche. Como puede ver, se está acercando al 95% del uso de la CPU. El rendimiento máximo de un solo hilo con ambos núcleos en 11.2 Gflops (Intel Turbo Boost de 2 pasos). Necesito apagar el otro núcleo para obtener 12.26 Gflops (4 pasos Intel Turbo Boost). Supongamos que las funciones de embalaje en OpeBLAS no generan una sobrecarga adicional. Entonces, el núcleo de OpenBLAS debe estar ejecutándose al menos tan rápido como el código total de OpenBLAS. Así que necesito hacer que mi núcleo funcione a esa velocidad. Todavía tengo que descubrir cómo hacer mi montaje más rápido. Me centraré en esto en los próximos días.

Ejecutó las siguientes pruebas desde la línea de comandos de Windows con: start /realtime /affinity 1 Mi Código:

N Total Gflops/s Kernel Gflops/s 256 7.927740 8.832366 512 8.427591 9.347094 768 8.547722 9.352993 1024 8.597336 9.351794 1280 8.588663 9.296724 1536 8.589808 9.271710 1792 8.634201 9.306406 2048 8.527889 9.235653

OpenBLAS:

N Total Gflops/s 256 10.599065 512 10.622686 768 10.745133 1024 10.762757 1280 10.832540 1536 10.793132 1792 10.848356 2048 10.819986

Teóricamente es posible ver ese código y razonar si se podría organizar para hacer un mejor uso de los recursos de microarquitectura, pero incluso los arquitectos de rendimiento de Intel podrían no recomendar hacerlo de esa manera. Podría ser útil usar una herramienta como VTune o Intel Performance Counter Monitor para saber qué parte de su carga de trabajo es memoria frente a front-end versus back-end. Intel Architecture Code Analyzer también puede ser una fuente rápida de ayuda para reducir cuál de los posibles problemas que se enumeran a continuación para realizar el seguimiento primero.

Es probable que Nominal Animal esté en el camino correcto en los comentarios que hablan sobre las instrucciones de intercalación que acceden a la memoria y las que hacen cálculos. Algunas otras posibilidades:

El uso de otras instrucciones para algunos de los cálculos puede reducir la presión en uno de los puertos de ejecución (consulte la sección 3.3.4 de esta presentación ). En particular, mulpd siempre va a enviar al puerto 1 en Westmere. Tal vez si hay algún ciclo en el que el puerto 0 no se esté utilizando, podría colarse en un FP escalar allí.
Uno u otro de los prefetchers de hardware podrían saturar el bus antes o contaminar la memoria caché con líneas que no usas .
Por otro lado, existe una pequeña posibilidad de que el orden de las referencias de memoria o el diseño de memoria implícito en dgemm_2x4_asm_j esté falsificando los prefetchers.
Cambiar el orden relativo de pares de instrucciones que no tienen ninguna dependencia de datos podría llevar a un mejor uso de los recursos de front-end o back-end.