¿Cómo logra BLAS un rendimiento tan extremo?

108

Por curiosidad, decidí comparar mi propia función de multiplicación de matrices con la implementación de BLAS ... Estaba por decir lo menos sorprendido con el resultado:

Implementación personalizada, 10 pruebas 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á utilizando números de punto 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:

  1. Dado que una multiplicación matriz-matriz dice: nxm * mxn requiere n * n * m multiplicaciones, entonces en el caso anterior 1000 ^ 3 o 1e9 operaciones. ¿Cómo es posible que BLAS en mi procesador de 2.6Ghz realice operaciones 10 * 1e9 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.
  2. ¿Por qué mi implementación es mucho más lenta?
DeusAduro
fuente
17
BLAS ha sido optimizado hacia arriba y hacia abajo por el otro especialista en la materia. Supongo que está aprovechando la unidad de punto flotante SIMD en su chip y jugando muchos trucos para mejorar el comportamiento de almacenamiento en caché también ...
dmckee --- ex-moderador gatito
3
Aún así, ¿cómo se realizan las operaciones 1E10 en un procesador de 2.63E9 ciclos / segundo en 1.3 segundos?
DeusAduro
9
Múltiples unidades de ejecución, revestimiento de tuberías y datos múltiples de instrucción única ((SIMD) que significa realizar la misma operación en más de un par de operandos al mismo tiempo). Algunos compiladores pueden apuntar a las unidades SIMD en chips comunes, pero casi siempre tienes que activarlos explícitamente, y es útil saber cómo funciona todo ( en.wikipedia.org/wiki/SIMD ). Asegurarse contra pérdidas de caché es casi con certeza la parte difícil.
dmckee --- exmoderador gatito
13
La suposición es incorrecta. Hay mejores algoritmos conocidos, consulte Wikipedia.
MSalters
2
@DeusAduro: En mi respuesta a ¿Cómo escribir un producto de matriz matricial que pueda competir con Eigen? Publiqué un pequeño ejemplo sobre cómo implementar un producto matriz-matriz eficiente en caché.
Michael Lehn

Respuestas:

141

Un buen punto de partida es el gran libro La ciencia de la programación de cálculos matriciales de Robert A. van de Geijn y Enrique S. Quintana-Ortí. Proporcionan una versión de descarga gratuita.

BLAS se divide en tres niveles:

  • El nivel 1 define un conjunto de funciones de álgebra lineal que operan solo 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 haga uso de alguna arquitectura de multiprocesador con memoria compartida.

  • Las funciones de nivel 3 son operaciones como el producto matriz-matriz. Nuevamente, podría implementarlos en términos de funciones de Nivel2. 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 optimizada para caché / compatible con caché . Esto está muy bien descrito en el libro. El principal impulso de las funciones de Level3 proviene de la optimización de la caché. Este impulso supera 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 de BLAS proporcionan una interfaz Fortran de modo que se puede vincular contra LAPACK (LAPACK obtiene todo su rendimiento de BLAS).

Los compiladores optimizados juegan un papel menor a este respecto (y para GotoBLAS / OpenBLAS el compilador no importa en absoluto).

En mi humilde opinión, la implementación de BLAS no utiliza 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ía más de lo que ganaría)
  • Estos algoritmos son numéricamente no estables. Como BLAS es el núcleo computacional de LAPACK, esto no es posible.

Editar / Actualizar:

El documento nuevo e innovador para este tema son los documentos BLIS . Están excepcionalmente bien escritos. Para mi conferencia "Conceptos básicos de software para computación de alto rendimiento", implementé el producto matriz-matriz después de su artículo. De hecho, 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 demás 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. Aproximadamente el 99,9% del tiempo se pasa 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, resulta bastante fácil comprender 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):

Puntos de referencia para ulmBLAS, BLIS, MKL, openBLAS y Eigen

Otra edición / actualización:

También escribí un tutorial sobre cómo se usa BLAS para problemas de álgebra lineal numérica, como resolver un sistema de ecuaciones lineales:

Factorización LU de alto rendimiento

(Esta factorización LU la utiliza, por ejemplo, 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: Codificación de una factorización de LU paralela optimizada de caché

PD: También hice algunos experimentos para mejorar el rendimiento de uBLAS. En realidad, es bastante simple impulsar (sí, jugar con las palabras :)) el rendimiento de uBLAS:

Experimentos sobre uBLAS .

Aquí un proyecto similar con BLAZE :

Experimentos en BLAZE .

Michael Lehn
fuente
3
Nuevo enlace a “Benchmarks for ulmBLAS, BLIS, MKL, openBLAS y Eigen”: apfel.mathematik.uni-ulm.de/~lehn/ulmBLAS/#toc3
Ahmed Fasih
Resulta que el ESSL de IBM utiliza una variación del algoritmo Strassen: ibm.com/support/knowledgecenter/en/SSFHY8/essl_welcome.html
ben-albrecht
2
la mayoría de los enlaces están muertos
Aurélien Pierre
Se puede encontrar un PDF de TSoPMC en la página del autor, en cs.utexas.edu/users/rvdg/tmp/TSoPMC.pdf
Alex Shpilkin
Aunque el algoritmo Coppersmith-Winograd tiene una gran complejidad de tiempo en papel, la notación Big O oculta una constante muy grande, por lo que solo comienza a ser viable para matrices ridículamente grandes.
DiehardThe Tryhard
26

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 medida no están relacionadas:

  • Fortran vs C, no importa
  • Los algoritmos matriciales avanzados como Strassen, las implementaciones no los utilizan porque no ayudan en la práctica.

La mayoría de las implementaciones dividen cada operación en operaciones matriciales o vectoriales de pequeñas dimensiones 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 kernels) están codificadas en código de ensamblaje específico de la CPU utilizando varias características de la CPU de su destino:

  • Instrucciones de estilo SIMD
  • Paralelismo a nivel de instrucción
  • Conciencia de caché

Además, estos núcleos se pueden ejecutar en paralelo entre sí utilizando varios subprocesos (núcleos de CPU), en el patrón de diseño de reducción de mapa típico.

Eche un vistazo a ATLAS, que es la implementación de BLAS de código abierto más utilizada. Tiene muchos kernels competidores diferentes, y durante el proceso de construcción de la biblioteca ATLAS se ejecuta una competencia entre ellos (algunos incluso están parametrizados, por lo que el mismo kernel puede tener configuraciones diferentes). Prueba diferentes configuraciones y luego selecciona la mejor para el sistema de destino en particular.

(Sugerencia: es por eso que si está usando ATLAS, es mejor construir y ajustar la biblioteca a mano para su máquina en particular y luego usar una prediseñada).

Andrés Tomazos
fuente
ATLAS ya no es la implementación BLAS de código abierto más utilizada. Ha sido superado por OpenBLAS (una bifurcación de GotoBLAS) y BLIS (una refactorización de GotoBLAS).
Robert van de Geijn
1
@ ulaff.net: Eso tal vez. Esto fue escrito hace 6 años. Creo que la implementación de BLAS más rápida actualmente (en Intel, por supuesto) es Intel MKL, pero no es de código abierto.
Andrew Tomazos
14

Primero, hay algoritmos más eficientes para la multiplicación de matrices que el que está usando.

En segundo lugar, su CPU puede realizar más de una instrucción a la vez.

Su CPU ejecuta 3-4 instrucciones por ciclo, y si se utilizan 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ás usando punteros sin procesar, lo que significa que el compilador debe asumir que pueden usar un alias. Hay palabras clave o indicadores específicos del compilador que puede especificar para decirle al compilador que no tienen alias. Alternativamente, debe usar otros tipos que no sean punteros sin procesar, que se encargan del problema.
  • Estás destruyendo la 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 cabe en la memoria caché de la CPU, antes de pasar al siguiente bloque.
  • Para tareas puramente numéricas, Fortran es prácticamente imbatible, y C ++ requiere mucha persuasión para alcanzar una velocidad similar. Se puede hacer, y hay algunas bibliotecas que lo demuestran (generalmente usando plantillas de expresión), pero no es trivial, y no sucede simplemente .
jalf
fuente
Gracias, agregué restringir el código correcto según la sugerencia de Justicle, no vi mucha mejora, me gusta la idea de bloque. Por curiosidad, sin conocer el tamaño de la caché de la CPU, ¿cómo sería correcto el código óptimo?
DeusAduro
2
Tu no Para obtener un código óptimo, necesita conocer el tamaño de la caché de la CPU. Por supuesto, la desventaja de esto es que está codificando efectivamente su código para obtener el mejor rendimiento en una familia de CPU.
jalf
2
Al menos el bucle interior aquí evita cargas escalonadas. Parece que esto está escrito para una matriz que ya se está transponiendo. ¡Por eso es "sólo" un orden de magnitud más lento que BLAS! Pero sí, todavía está mal por la falta de bloqueo de caché. ¿Estás seguro de que Fortran ayudaría mucho? Creo que todo lo que ganaría aquí es que restrict(sin alias) es el valor predeterminado, a diferencia de C / C ++. (Y desafortunadamente ISO C ++ no tiene una restrictpalabra clave, por lo que debe usarla __restrict__en compiladores que la proporcionen como una extensión).
Peter Cordes
11

No sé específicamente sobre la implementación de BLAS, pero hay algoritmos más eficientes para la multiplicación de matrices que tienen una complejidad mejor que O (n3). Uno bien conocido es el algoritmo de Strassen

softveda
fuente
8
El algoritmo de Strassen no se usa en numéricos por dos razones: 1) No es estable. 2) Ahorras algunos cálculos, pero eso tiene el precio de poder aprovechar las jerarquías de caché. En la práctica, incluso pierdes rendimiento.
Michael Lehn
4
Para la implementación práctica del Algoritmo Strassen estrechamente construido sobre el código fuente de la biblioteca BLAS, hay una publicación reciente: " Algoritmo Strassen Reloaded " en SC16, que logra un rendimiento más alto que BLAS, incluso para el tamaño del problema 1000x1000.
Jianyu Huang
4

La mayoría de los argumentos para la segunda pregunta - ensamblador, división en bloques, etc. (pero no menos de 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 a la vez en la memoria caché. Puede reorganizar los bucles de manera que se haga todo lo posible en una fila en la caché, de esta manera reduciendo drásticamente las actualizaciones de la caché (por cierto, la división en bloques pequeños tiene un efecto analógico, mejor si los bucles sobre los bloques se organizan de manera similar). A continuación, se muestra una implementación de modelo para matrices cuadradas. En mi computadora, su consumo de tiempo fue de 1:10 en comparación con la implementación estándar (como la suya). En otras palabras: nunca programe una multiplicación de matrices a lo largo del "

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

Un comentario más: esta implementación es incluso mejor en mi computadora que reemplazar todo por la rutina BLAS cblas_dgemm (¡pruébalo en tu computadora!). Pero mucho más rápido (1: 4) es llamar directamente a dgemm_ de la biblioteca Fortran. Creo que esta rutina no es de hecho Fortran sino código ensamblador (no sé qué hay en la biblioteca, no tengo las fuentes). No me queda muy claro por qué cblas_dgemm no es tan rápido, ya que, que yo sepa, es simplemente un contenedor para dgemm_.

Wolfgang Jansen
fuente
3

Esta es una aceleración realista. Para ver un ejemplo de lo que se puede hacer con el ensamblador SIMD sobre el código C ++, vea algunos ejemplos de funciones de matriz de iPhone : eran 8 veces más rápidas que la versión C y ni siquiera el ensamblaje "optimizado", todavía no hay revestimiento de tuberías y hay son operaciones de pila innecesarias.

Además, su código no es " restringido correcto " - ¿cómo sabe el compilador que cuando modifica C, no modifica A y B?

Justicle
fuente
Seguro que si llamaste a la función como mmult (A ..., A ..., A); ciertamente no obtendría el resultado esperado. Una vez más, aunque no estaba tratando de vencer / volver a implementar BLAS, solo viendo lo rápido que realmente es, por lo que la verificación de errores no estaba en mente, solo la funcionalidad básica.
DeusAduro
3
Lo siento, para que quede claro, lo que estoy diciendo es que si pones "restringir" en tus punteros, obtendrás un código mucho más rápido. Esto se debe a que cada vez que modifica C, el compilador no tiene que volver a cargar A y B, lo que acelera drásticamente el ciclo interno. Si no me cree, compruebe el desmontaje.
Justicle
@DeusAduro: esto no es una comprobación de errores; es posible que el compilador no pueda optimizar los accesos a la matriz B [] en el bucle interno porque es posible que no pueda averiguar que los punteros A y C nunca alias B formación. Si hubiera un alias, sería posible que el valor en la matriz B cambiara mientras se ejecuta el ciclo interno. Sacar el acceso al valor B [] del bucle interno y colocarlo en una variable local podría permitir al compilador evitar accesos continuos a B [].
Michael Burr
1
Hmmm, primero intenté usar la palabra clave '__restrict' en VS 2008, aplicada a A, B y C. Esto no mostró cambios en el resultado. Sin embargo, mover el acceso a B, desde el bucle más interno al bucle exterior mejoró el tiempo en un ~ 10%.
DeusAduro
1
Lo siento, no estoy seguro de VC, pero con GCC debes habilitarlo -fstrict-aliasing. También hay una mejor explicación de "restringir" aquí: cellperformance.beyond3d.com/articles/2006/05/…
Justicle del
2

Con respecto al código original en MM Multiply, la referencia de memoria para la mayoría de las operaciones es la principal causa del mal rendimiento. La memoria funciona a 100-1000 veces más lento que el 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 multiplicar. 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 los bloqueamos para la reutilización de datos en la caché. El desenrollado del bucle externo ayuda a optimizar el acceso a los datos temporalmente 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 caché. Puede optar por optimizar la caché L2 o la caché L3.

https://en.wikipedia.org/wiki/Loop_nest_optimization

Pari Rajaram
fuente
-24

Por muchas razones.

Primero, los compiladores de 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 están más simplificados y el compilador tiene un mejor control de lo que sucede, lo que le permite optimizar más (por ejemplo, utilizando registros).

Otra cosa es que Fortran almacena las cosas en columnas, mientras que C almacena los datos en filas. No he comprobado su código, pero tenga cuidado con el rendimiento del producto. En C, debe escanear por filas: de esta manera escanea su matriz a lo largo de la memoria contigua, reduciendo las pérdidas 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 estar escritas en ensamblador y optimizadas 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 de ellas 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 operaciones costosas.

Por ejemplo, puede 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ébelo, estoy seguro de que guardará algo.

En su pregunta # 1, la razón es que la multiplicación de matrices escala como O (n ^ 3) si usa un algoritmo trivial. Hay algoritmos que escalan mucho mejor .

Stefano Borini
fuente
36
Esta respuesta es completamente incorrecta, lo siento. Las implementaciones de BLAS no están escritas en fortran. El código crítico para el rendimiento está escrito en ensamblador, y los más comunes en estos días están escritos en C arriba. Además, BLAS especifica el orden de fila / columna como parte de la interfaz, y las implementaciones pueden manejar cualquier combinación.
Andrew Tomazos
10
Sí, esta respuesta es completamente incorrecta. Desafortunadamente, está lleno de tonterías comunes, por ejemplo, la afirmación de que BLAS fue más rápido debido a Fortran. Tener 20 (!) Valoraciones positivas es algo malo. ¡Ahora este sinsentido se extiende aún más debido a la popularidad de Stackoverflow!
Michael Lehn
12
Creo que está confundiendo la implementación de referencia no optimizada con las implementaciones de producción. La implementación de referencia es solo para especificar la interfaz y el comportamiento de la biblioteca, y fue escrita en Fortran por razones históricas. No es para uso de producción. En producción, la gente usa implementaciones optimizadas que exhiben el mismo comportamiento que la implementación de referencia. He estudiado los aspectos internos de ATLAS (que respalda a Octave - Linux "MATLAB") que puedo confirmar de primera mano que está escrito en C / ASM internamente. Es casi seguro que las implementaciones comerciales también lo sean.
Andrew Tomazos
5
@KyleKanos: Sí, aquí está la fuente de ATLAS: sourceforge.net/projects/math-atlas/files/Stable/3.10.1 Hasta donde yo sé, es la implementación BLAS portátil de código abierto más utilizada. Está escrito en C / ASM. Los fabricantes de CPU de alto rendimiento, como Intel, también ofrecen implementaciones BLAS especialmente optimizadas para sus chips. Garantizo que las partes de bajo nivel de la biblioteca de Intels están escritas en (duuh) ensamblaje x86, y estoy bastante seguro de que las partes de nivel medio se escribirían en C o C ++.
Andrew Tomazos
9
@KyleKanos: Estás confundido. Netlib BLAS es la implementación de referencia. La implementación de referencia es mucho más lenta que las implementaciones optimizadas (consulte la comparación de rendimiento ). Cuando alguien dice que está usando netlib BLAS en un clúster, no significa que realmente esté usando la implementación de referencia de netlib. Eso sería una tontería. Simplemente significa que están usando una lib con la misma interfaz que netlib blas.
Andrew Tomazos