c++ - ¿Cómo consigue BLAS un rendimiento tan extremo?
fortran (8)
Con respecto al código original en MM multiplicado, la referencia de memoria para la mayoría de las operaciones es la causa principal del mal rendimiento. La memoria se ejecuta entre 100 y 1000 veces más lenta que la caché.
La mayor parte de la aceleración proviene del empleo de técnicas de optimización de bucle para esta función de triple bucle en MM multiplicado. Se utilizan dos técnicas principales de optimización de bucle; desenrollar y bloquear. Con respecto al desenrollado, desenrollamos los dos bucles más externos y lo bloqueamos para la reutilización de datos en caché. El desenrollado del bucle externo ayuda a optimizar temporalmente el acceso a los datos al reducir el número de referencias de memoria a los mismos datos en diferentes momentos durante toda la operación. Bloquear el índice de bucle en un número específico, ayuda a retener los datos en la memoria caché. Puede elegir optimizar para caché L2 o caché L3.
Por curiosidad, decidí comparar mi propia función de multiplicación de matrices con la implementación de BLAS ... Por lo menos, me sorprendió el resultado:
Implementación personalizada, 10 ensayos de multiplicación de matrices 1000x1000:
Took: 15.76542 seconds.Implementación BLAS, 10 ensayos de multiplicación de matrices 1000x1000:
Took: 1.32432 seconds.
Esto está usando números de coma flotante de precisión simple.
Mi implementación:
template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 )
throw std::runtime_error("Error sizes off");
memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}
Tengo dos preguntas:
- Dado que una multiplicación matriz-matriz dice: nxm * mxn requiere n * n * m multiplicaciones, por lo que en el caso anterior 1000 ^ 3 o 1e9 operaciones. ¿Cómo es posible en mi procesador de 2.6Ghz para BLAS hacer 10 * 1e9 operaciones en 1.32 segundos? Incluso si las multiplicaciones fueran una sola operación y no se hiciera nada más, debería tomar ~ 4 segundos.
- ¿Por qué mi implementación es mucho más lenta?
Entonces, en primer lugar, BLAS es solo una interfaz de aproximadamente 50 funciones. Hay muchas implementaciones competitivas de la interfaz.
En primer lugar mencionaré cosas que en gran parte no están relacionadas:
- Fortran vs C, no hace diferencia
- Algoritmos de matriz avanzada como Strassen, las implementaciones no los usan ya que no ayudan en la práctica
La mayoría de las implementaciones dividen cada operación en operaciones de matriz o vector de pequeña dimensión de la manera más o menos obvia. Por ejemplo, una gran multiplicación de matrices de 1000x1000 puede dividirse en una secuencia de multiplicaciones de matrices de 50x50.
Estas operaciones de pequeña dimensión de tamaño fijo (llamadas núcleos) están codificadas en código ensamblador específico de CPU usando varias características de CPU de su objetivo:
- Instrucciones estilo SIMD
- Paralelismo de nivel de instrucción
- Conocimiento de la caché
Además, estos núcleos se pueden ejecutar en paralelo entre sí utilizando múltiples hilos (núcleos de CPU), en el patrón típico de diseño de reducción de mapa.
Eche un vistazo a ATLAS que es la implementación BLAS de código abierto más comúnmente utilizada. Tiene muchos núcleos que compiten entre sí, y durante el proceso de compilación de la biblioteca ATLAS ejecuta una competencia entre ellos (algunos incluso están parametrizados, por lo que el mismo núcleo puede tener configuraciones diferentes). Prueba diferentes configuraciones y luego selecciona la mejor para el sistema de destino particular.
(Sugerencia: Es por eso que si está utilizando ATLAS, es mejor que construya y ajuste la biblioteca a mano para su máquina en particular y luego utilice una precompilada).
Esta es una aceleración realista. Para ver un ejemplo de lo que se puede hacer con el ensamblador SIMD sobre código C ++, vea algunas funciones de matriz de iPhone, por ejemplo, que eran 8 veces más rápidas que la versión C y ni siquiera están "optimizadas". Aún no hay revestimiento. es una operación de pila innecesaria.
Además, su código no es " restrictivo ": ¿cómo sabe el compilador que cuando modifica C, no modifica A y B?
La mayoría de los argumentos a la segunda pregunta: ensamblador, división en bloques, etc. (pero no menos que N ^ 3 algoritmos, están realmente sobredesarrollados): juegan un papel. Pero la baja velocidad de su algoritmo es causada esencialmente por el tamaño de la matriz y la desafortunada disposición de los tres bucles anidados. Sus matrices son tan grandes que no caben de una vez en la memoria caché. Puede reorganizar los bucles de manera que se haga todo lo posible en una fila en caché, reduciendo dramáticamente las actualizaciones de caché (por cierto, dividir en pequeños bloques tiene un efecto analógico, mejor si los bucles sobre los bloques se organizan de manera similar). Una implementación modelo para matrices cuadradas sigue. En mi computadora, su consumo de tiempo fue de aproximadamente 1:10 en comparación con la implementación estándar (como la suya). En otras palabras: nunca programe una multiplicación de la matriz a lo largo del esquema de "columna de filas de filas" que aprendimos en la escuela. Después de haber reorganizado los bucles, se obtienen más mejoras desenrollando bucles, código de ensamblador, etc.
void vector(int m, double ** a, double ** b, double ** c) {
int i, j, k;
for (i=0; i<m; i++) {
double * ci = c[i];
for (k=0; k<m; k++) ci[k] = 0.;
for (j=0; j<m; j++) {
double aij = a[i][j];
double * bj = b[j];
for (k=0; k<m; k++) ci[k] += aij*bj[k];
}
}
}
Una observación más: esta implementación es aún mejor en mi computadora que reemplazando todo por la rutina BLAS cblas_dgemm (pruébelo en su computadora). Pero mucho más rápido (1: 4) está llamando dgemm_ de la biblioteca Fortran directamente. Creo que esta rutina no es, en realidad, Fortran, sino un código ensamblador (no sé qué hay en la biblioteca, no tengo las fuentes). No tengo muy claro por qué cblas_dgemm no es tan rápido ya que, a mi entender, es simplemente un contenedor para dgemm_.
No sé específicamente sobre la implementación de BLAS, pero existen algoritmos más eficientes para la multiplicación de matrices que tienen una complejidad superior a O (n3). Una muy conocida es el algoritmo de Strassen
Por muchas razones.
Primero, los compiladores Fortran están altamente optimizados, y el lenguaje les permite ser como tales. C y C ++ son muy flexibles en términos de manejo de matrices (por ejemplo, el caso de punteros que se refieren a la misma área de memoria). Esto significa que el compilador no puede saber de antemano qué hacer y se ve obligado a crear un código genérico. En Fortran, sus casos son más optimizados, y el compilador tiene un mejor control de lo que sucede, lo que le permite optimizar más (por ejemplo, el uso de registros).
Otra cosa es que Fortran almacene cosas en columnas, mientras que C almacena datos en filas. No he revisado su código, pero tenga cuidado de cómo realiza el producto. En C debe escanear filas: de esta manera escanea su matriz a lo largo de la memoria contigua, reduciendo las fallas de caché. La falta de caché es la primera fuente de ineficiencia.
En tercer lugar, depende de la implementación de blas que esté utilizando. Algunas implementaciones pueden escribirse en ensamblador y optimizarse para el procesador específico que está utilizando. La versión de netlib está escrita en fortran 77.
Además, está realizando muchas operaciones, la mayoría repetidas y redundantes. Todas esas multiplicaciones para obtener el índice son perjudiciales para el rendimiento. Realmente no sé cómo se hace esto en BLAS, pero hay muchos trucos para evitar costosas operaciones.
Por ejemplo, podría volver a trabajar su código de esta manera
template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Error sizes off");
memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
a1 = cc2*ADim2;
a3 = cc2*BDim1
for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
a2=cc1*ADim1;
ValT b = B[a3+cc1];
for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
C[a1+cr1] += A[a2+cr1]*b;
}
}
}
}
Pruébalo, estoy seguro de que vas a guardar algo.
En su pregunta # 1, la razón es que la multiplicación de la matriz se escala como O (n ^ 3) si usa un algoritmo trivial. Hay algoritmos que escalan mucho mejor .
Primero, hay algoritmos más eficientes para la multiplicación de matrices que el que estás usando.
En segundo lugar, su CPU puede hacer mucho más que una instrucción a la vez.
Su CPU ejecuta 3-4 instrucciones por ciclo, y si se usan las unidades SIMD, cada instrucción procesa 4 flotantes o 2 dobles. (por supuesto, esta cifra tampoco es precisa, ya que la CPU normalmente solo puede procesar una instrucción SIMD por ciclo)
En tercer lugar, su código está lejos de ser óptimo:
- Está utilizando punteros crudos, lo que significa que el compilador tiene que asumir que pueden alias. Hay palabras clave o indicadores específicos del compilador que puede especificar para indicar al compilador que no tienen alias. Alternativamente, debe usar otros tipos de punteros sin procesar, que solucionan el problema.
- Estás agotando la memoria caché realizando un recorrido ingenuo de cada fila / columna de las matrices de entrada. Puede usar el bloqueo para realizar la mayor cantidad de trabajo posible en un bloque más pequeño de la matriz, que se ajusta en la memoria caché de la CPU, antes de pasar al siguiente bloque.
- Para tareas puramente numéricas, Fortran es bastante inmejorable, y C ++ requiere mucha persuasión para alcanzar una velocidad similar. Se puede hacer, y hay algunas bibliotecas que lo demuestran (normalmente usan plantillas de expresión), pero no es trivial, y no sucede solo .
Un buen punto de partida es el gran libro The Science of Programming Matrix Computations de Robert A. van de Geijn y Enrique S. Quintana-Ortí. Proporcionan una versión de descarga gratuita.
BLAS está dividido en tres niveles:
El nivel 1 define un conjunto de funciones de álgebra lineal que operan únicamente en vectores. Estas funciones se benefician de la vectorización (por ejemplo, del uso de SSE).
Las funciones de nivel 2 son operaciones matriz-vector, por ejemplo, algún producto matriz-vector. Estas funciones podrían implementarse en términos de funciones de Nivel1. Sin embargo, puede aumentar el rendimiento de estas funciones si puede proporcionar una implementación dedicada que utilice una arquitectura multiprocesador con memoria compartida.
Las funciones de nivel 3 son operaciones como el producto matriz-matriz. De nuevo, podría implementarlos en términos de funciones de Level2. Pero las funciones de Level3 realizan operaciones O (N ^ 3) en datos O (N ^ 2). Entonces, si su plataforma tiene una jerarquía de caché, puede aumentar el rendimiento si proporciona una implementación dedicada que sea compatible con la caché / la memoria caché . Esto está muy bien descrito en el libro. El impulso principal de las funciones de Level3 proviene de la optimización de caché. Este impulso excede significativamente el segundo impulso del paralelismo y otras optimizaciones de hardware.
Por cierto, la mayoría (o incluso todas) de las implementaciones BLAS de alto rendimiento NO se implementan en Fortran. ATLAS se implementa en C. GotoBLAS / OpenBLAS se implementa en C y sus partes críticas de rendimiento en Assembler. Solo la implementación de referencia de BLAS se implementa en Fortran. Sin embargo, todas estas implementaciones BLAS proporcionan una interfaz Fortran de modo que se puede vincular con LAPACK (LAPACK gana todo su rendimiento de BLAS).
Los compiladores optimizados desempeñan un papel menor a este respecto (y para GotoBLAS / OpenBLAS el compilador no importa en absoluto).
En mi humilde opinión, ninguna implementación de BLAS usa algoritmos como el algoritmo Coppersmith-Winograd o el algoritmo Strassen. No estoy exactamente seguro de la razón, pero esta es mi suposición:
- Tal vez no sea posible proporcionar una implementación optimizada de caché de estos algoritmos (es decir, perderías más de lo que ganarías)
- Estos algoritmos son numéricamente no estables. Como BLAS es el núcleo computacional de LAPACK, es un no-go.
Editar / actualizar:
El documento nuevo e innovador para este tema son los documentos BLIS . Están excepcionalmente bien escritos. Para mi conferencia "Fundamentos de software para computación de alto rendimiento" implementé el producto matriz-matriz siguiendo su papel. En realidad, implementé varias variantes del producto matriz-matriz. Las variantes más simples están escritas completamente en C simple y tienen menos de 450 líneas de código. Todas las otras variantes simplemente optimizan los bucles
for (l=0; l<MR*NR; ++l) {
AB[l] = 0;
}
for (l=0; l<kc; ++l) {
for (j=0; j<NR; ++j) {
for (i=0; i<MR; ++i) {
AB[i+j*MR] += A[i]*B[j];
}
}
A += MR;
B += NR;
}
El rendimiento general del producto matriz-matriz solo depende de estos bucles. Alrededor del 99.9% del tiempo lo pasamos aquí. En las otras variantes usé intrínsecos y código ensamblador para mejorar el rendimiento. Puedes ver el tutorial pasando por todas las variantes aquí:
ulmBLAS: Tutorial sobre GEMM (Producto Matrix-Matrix)
Junto con los documentos de BLIS, es bastante fácil entender cómo las bibliotecas como Intel MKL pueden obtener ese rendimiento. ¡Y por qué no importa si usa almacenamiento principal en filas o columnas!
Los puntos de referencia finales están aquí (llamamos a nuestro proyecto ulmBLAS):
Los puntos de referencia para ulmBLAS, BLIS, MKL, openBLAS y Eigen
Otra Edición / Actualización:
También escribí algunos tutoriales sobre cómo BLAS se usa para problemas numéricos de álgebra lineal, como resolver un sistema de ecuaciones lineales:
Factorización LU de alto rendimiento
(Esta factorización LU es utilizada, por ejemplo, por Matlab para resolver un sistema de ecuaciones lineales).
Espero encontrar tiempo para extender el tutorial para describir y demostrar cómo realizar una implementación paralela altamente escalable de la factorización LU como en PLASMA .
Ok, aquí tienes: Codificar una factorización de LU paralela optimizada en caché
PD: También realicé algunos experimentos para mejorar el rendimiento de uBLAS. En realidad, es bastante simple impulsar (sí, jugar con palabras :)) el rendimiento de uBLAS:
Aquí un proyecto similar con BLAZE :