Multiplicación de matrices optimizada en C

Estoy tratando de comparar diferentes métodos para la multiplicación de matrices. El primero es el método normal:

do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            suma = 0;
            for (l = 0; l < i; l++)
                suma += MatrixA[j][l]*MatrixB[l][k];
                MatrixR[j][k] = suma;
            }
        }
    }
    c++;
} while (c<iteraciones);

La segunda consiste en transponer primero la matriz B y luego hacer la multiplicación por filas:

int f, co;
for (f = 0; f < i; f++) {
    for ( co = 0; co < i; co++) {
        MatrixB[f][co] = MatrixB[co][f];
    }
}

c = 0;
do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            suma = 0;
            for (l = 0; l < i; l++)
                suma += MatrixA[j][l]*MatrixB[k][l];
                MatrixR[j][k] = suma;
            }
        }
     }
     c++;
} while (c<iteraciones);

Se supone que el segundo método es mucho más rápido, porque estamos accediendo a ranuras de memoria contiguas, pero no obtengo una mejora significativa en el rendimiento. ¿Estoy haciendo algo mal?

Puedo publicar el código completo, pero creo que no es necesario.

Lo que todo programador debe saber sobre la memoria (enlace pdf) por Ulrich Drepper tiene muchas buenas ideas sobre la eficiencia de la memoria, pero en particular, usa la multiplicación de matrices como un ejemplo de cómo conocer la memoria y usar ese conocimiento puede acelerar este proceso. Consulte el apéndice A.1 de su artículo y lea la sección 6.2.1. La Tabla 6.2 del documento muestra que podría obtener un tiempo de ejecución del 10 % a partir del tiempo de una implementación ingenua para una matriz de 1000×1000.

De acuerdo, su código final es bastante peludo y usa muchas cosas específicas del sistema y ajustes en tiempo de compilación, pero aun así, si De Verdad necesita velocidad, leer ese documento y leer su implementación definitivamente vale la pena.

Hacer esto bien no es trivial. Se recomienda encarecidamente utilizar una biblioteca BLAS existente.

Si realmente te inclinas a rodar tu propia multiplicación de matrices, mosaico de bucle es una optimización que es de particular importancia para matrices grandes. El mosaico debe ajustarse al tamaño de la memoria caché para garantizar que la memoria caché no esté siendo golpeada continuamente, lo que ocurrirá con una implementación ingenua. Una vez medí una diferencia de rendimiento de 12x al combinar una matriz multiplicada con tamaños de matriz seleccionados para consumir múltiplos de mi caché (alrededor del ’97, por lo que el caché probablemente era pequeño).

Los algoritmos de mosaico en bucle suponen que se utiliza una matriz lineal contigua de elementos, en lugar de filas o columnas de punteros. Con una opción de almacenamiento de este tipo, su esquema de indexación determina qué dimensión cambia más rápido y usted es libre de decidir si el acceso a filas o columnas tendrá el mejor rendimiento de caché.

Hay una lote de literatura sobre el tema. Las siguientes referencias, especialmente los libros de Banerjee, pueden ser útiles:

[Ban93] Banerjee, Utpal, Loop Transformations for Restructuring Compilers: the Foundations, Kluwer Academic Publishers, Norwell, MA, 1993.

[Ban94] Banerjee, Utpal, Paralelización de bucles, Kluwer Academic Publishers, Norwell, MA, 1994.

[BGS93] Bacon, David F., Susan L. Graham y Oliver Sharp, Compiler Transformations for High-Performance Computing, Computer Science Division, University of California, Berkeley, Calif., Technical Report No UCB/CSD-93-781.

[LRW91] Lam, Monica S., Edward E. Rothberg y Michael E Wolf. The Cache Performance and Optimizations of Blocked Algorithms, en la 4ª Conferencia Internacional sobre Soporte Arquitectónico para Lenguajes de Programación, celebrada en Santa Clara, California, abril de 1991, 63-74.

[LW91] Lam, Monica S. y Michael E Wolf. Una teoría de transformación de bucles y un algoritmo para maximizar el paralelismo, en IEEE Transactions on Parallel and Distributed Systems, 1991, 2(4):452-471.

[PW86] Padua, David A. y Michael J. Wolfe, Optimizaciones avanzadas del compilador para supercomputadoras, en Comunicaciones de la ACM, 29(12):1184-1201, 1986.

[Wolfe89] Wolfe, Michael J. Optimización de supercompiladores para supercomputadoras, The MIT Press, Cambridge, MA, 1989.

[Wolfe96] Wolfe, Michael J., Compiladores de alto rendimiento para computación en paralelo, Addison-Wesley, CA, 1996.

ATENCIÓN: Tienes un BUG en tu segunda implementación

for (f = 0; f < i; f++) {
    for (co = 0; co < i; co++) {
        MatrixB[f][co] = MatrixB[co][f];
    }
}

Cuando haces f=0, c=1

        MatrixB[0][1] = MatrixB[1][0];

tu sobrescribes MatrixB[0][1] y perder ese valor! Cuando el ciclo llega a f=1, c=0

        MatrixB[1][0] = MatrixB[0][1];

el valor copiado es el mismo que ya estaba ahí.

Si la matriz no es lo suficientemente grande o no repites las operaciones un número elevado de veces no verás diferencias apreciables.

Si la matriz es, digamos, de 1.000×1.000 empezarás a ver mejoras, pero yo diría que si está por debajo de 100×100 no debes preocuparte.

Además, cualquier ‘mejora’ puede ser del orden de milisegundos, a menos que esté trabajando con matrices extremadamente grandes o repitiendo la operación miles de veces.

Finalmente, si cambias la computadora que estás usando por una más rápida, ¡las diferencias serán aún más pequeñas!

No debes escribir la multiplicación de matrices. Debe depender de bibliotecas externas. En particular, debe utilizar el GEMM rutina de la BLAS Biblioteca. GEMM a menudo proporciona las siguientes optimizaciones

Bloqueo

La multiplicación eficiente de matrices se basa en bloquear su matriz y realizar varias multiplicaciones bloqueadas más pequeñas. Idealmente, el tamaño de cada bloque se elige para que encaje bien en el caché. mejorando mucho el rendimiento.

Afinación

El tamaño de bloque ideal depende de la jerarquía de memoria subyacente (¿qué tan grande es el caché?). Como resultado, las bibliotecas deben ajustarse y compilarse para cada máquina específica. Esto lo hace, entre otros, el ATLAS implementación de BLAS.

Optimización del nivel de ensamblaje

La multiplicación de matrices es tan común que los desarrolladores la optimizarán manualmente. En particular, esto se hace en GotoBLAS.

Computación heterogénea (GPU)

Matrix Multiply es muy intensivo en FLOP/computación, lo que lo convierte en un candidato ideal para ejecutarse en GPU. cuBLAS y MAGMA son buenos candidatos para esto.

En resumen, el álgebra lineal densa es un tema bien estudiado. La gente dedica su vida a la mejora de estos algoritmos. Deberías usar su trabajo; Los hará felices.

Las grandes mejoras que obtenga dependerán de:

1. El tamaño del caché
2. El tamaño de una línea de caché
3. El grado de asociatividad de la memoria caché.

Para tamaños de matriz pequeños y procesadores modernos, es muy probable que los datos de ambos MatrixA y MatrixB se mantendrá casi por completo en el caché después de tocarlo por primera vez.

Solo algo para que pruebe (pero esto solo marcaría la diferencia para matrices grandes): separe su lógica de suma de la lógica de multiplicación en el ciclo interno de la siguiente manera:

for (k = 0; k < i; k++)
{
    int sums[i];//I know this size declaration is illegal in C. consider 
            //this pseudo-code.
    for (l = 0; l < i; l++)
        sums[l] = MatrixA[j][l]*MatrixB[k][l];

    int suma = 0;
    for(int s = 0; s < i; s++)
       suma += sums[s];
}

Esto se debe a que terminas estancando tu tubería cuando escribes a suma. De acuerdo, gran parte de esto se soluciona en el cambio de nombre de registro y cosas por el estilo, pero con mi comprensión limitada del hardware, si quisiera exprimir cada gramo de rendimiento del código, lo haría porque ahora no tiene que hacerlo. detenga la canalización para esperar una escritura en suma. Dado que la multiplicación es más costosa que la suma, desea dejar que la máquina la paralelice tanto como sea posible, por lo que guardar sus puestos para la suma significa que pasará menos tiempo esperando en el ciclo de suma que en el ciclo de multiplicación.

Esta es solo mi lógica. Otros con más conocimiento en el área pueden estar en desacuerdo.