Pregunta ¿Cómo consigue BLAS un rendimiento tan extremo?


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 de 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:

  1. 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.
  2. ¿Por qué mi implementación es mucho más lenta?

75
2017-08-19 23:30


origen


Respuestas:


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 escala mucho mejor.


-24
2017-08-19 23:36



Un buen punto de partida es el gran libro La ciencia de los cálculos de la matriz de programación por 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 matriciales-vectoriales, p. 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 caché optimizado / caché amigable. 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 nuevo y pionero documento para este tema es el 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 solamente 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 el tiempo 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 en caché

P.S .: 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:

Experimentos en uBLAS.

Aquí un proyecto similar con RESPLANDOR:

Experimentos sobre BLAZE.


98
2017-07-10 20:23



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).


19
2018-06-26 14:10



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 lo hace sólo ocurrir.

12
2017-08-20 12:12



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 Algoritmo Strassen 


9
2017-08-20 00:15



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_.


4
2017-11-30 20:11



Esta es una aceleración realista. Para ver un ejemplo de lo que se puede hacer con ensamblador SIMD sobre código C ++, vea algún ejemplo funciones de matriz de iPhone - Estos fueron más de 8 veces más rápidos que la versión C, y ni siquiera están ensamblados "optimizados". Aún no hay revestimiento de tubería y no hay operaciones de pila innecesarias.

También tu código no es "restringir correcto"¿Cómo sabe el compilador que cuando modifica C, no modifica A y B?


3
2017-08-20 00:10



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.

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


1
2018-05-02 12:07