una - Multiplicacion de matrices optimizada en C

multiplicacion de vectores en c

Question

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

El segundo consiste en transponer la matriz B primero 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 estoy obteniendo una mejora significativa en el rendimiento. ¿Estoy haciendo algo mal?

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

Answer 1

¿Puede publicar algunos datos comparando sus 2 enfoques para un rango de tamaños de matriz? Es posible que sus expectativas no sean realistas y que su segunda versión sea más rápida pero aún no haya realizado las mediciones.

No olvide, al medir el tiempo de ejecución, incluir el tiempo para transponer la matriz B.

Otra cosa que quizás desee probar es comparar el rendimiento de su código con el de la operación equivalente de su biblioteca BLAS. Es posible que esto no responda directamente a su pregunta, pero le dará una mejor idea de lo que puede esperar de su código.

Answer 2

Cómo las grandes mejoras que obtengas dependerán de:

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

Para tamaños de matriz pequeños y procesadores modernos, es muy probable que los datos de MatrixA y MatrixB se mantengan casi en su totalidad en el caché después de tocarlos por primera vez.

Answer 3

En general, la transposición de B debería ser mucho más rápida que la implementación ingenua, pero a expensas de desperdiciar otro NxN de memoria. Acabo de pasar una semana explorando la optimización de la multiplicación de matrices, y hasta ahora el ganador absoluto es este:

for (int i = 0; i < N; i++)
    for (int k = 0; k < N; k++)
        for (int j = 0; j < N; j++)
            if (likely(k)) /* #define likely(x) __builtin_expect(!!(x), 1) */
                C[i][j] += A[i][k] * B[k][j];
            else
                C[i][j] = A[i][k] * B[k][j];

Esto es incluso mejor que el método de Drepper mencionado en un comentario anterior, ya que funciona de manera óptima independientemente de las propiedades de caché de la CPU subyacente. El truco está en reordenar los bucles para que se pueda acceder a las tres matrices en orden de fila mayor.

Answer 4

Hacer esto bien puede ser no trivial. Una optimización que es de particular importancia para matrices grandes, es agrupar la multiplicación para mantener las cosas en el caché. Una vez medí una diferencia de rendimiento de 12x al hacerlo, pero elegí específicamente un tamaño de matriz que consumía múltiplos de mi caché (alrededor del ''97, por lo que el caché era pequeño).

Hay mucha literatura sobre el tema. Un punto de partida es:

http://en.wikipedia.org/wiki/Loop_tiling

Para un estudio más profundo, las siguientes referencias, especialmente los libros de Banerjee, pueden ser útiles:

[Ban93] Banerjee, Utpal, Transformaciones en bucle para reestructurar compiladores: The Foundations, Kluwer Academic Publishers, Norwell, MA, 1993.

[Ban94] Banerjee, Utpal, Loop Parallelization, Kluwer Academic Publishers, Norwell, MA, 1994.

[BGS93] Bacon, David F., Susan L. Graham y Oliver Sharp, Transformaciones de compilación para computación de alto rendimiento, División de informática, Universidad de California, Berkeley, California, Informe técnico No UCB / CSD-93-781.

[LRW91] Lam, Monica S., Edward E. Rothberg y Michael E Wolf. El rendimiento de caché y las optimizaciones de los algoritmos bloqueados, en la 4ta 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 la transformación de bucle y un algoritmo para maximizar el paralelismo, en transacciones IEEE en sistemas paralelos y distribuidos, 1991, 2 (4): 452-471.

[PW86] Padua, David A. y Michael J. Wolfe, Optimizaciones avanzadas de compilación para supercomputadoras, En comunicaciones de la ACM, 29 (12): 1184-1201, 1986.

[Wolfe89] Wolfe, Michael J. Optimizando supercompiladores para supercomputadoras, The MIT Press, Cambridge, MA, 1989.

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

Answer 5

La complejidad de cálculo de la multiplicación de dos matrices N * N es O (N ^ 3). El rendimiento mejorará dramáticamente si utiliza el algoritmo O (N ^ 2.73) que probablemente haya sido adoptado por MATLAB. Si instaló un MATLAB, intente multiplicar dos matrices 1024 * 1024. En mi computadora, MATLAB lo completa en 0.7s, pero la implementación C / C ++ del algoritmo ingenuo como el tuyo toma 20s. Si realmente te importa el rendimiento, consulta los algoritmos de menor complejidad. Escuché que existe un algoritmo O (N ^ 2.4), sin embargo, se necesita una matriz muy grande para que otras manipulaciones puedan ser descuidadas.

Answer 6

Muy antigua pregunta, pero aquí está mi implementación actual para mis proyectos de opengl:

typedef float matN[N][N];

inline void matN_mul(matN dest, matN src1, matN src2)
{
    unsigned int i;
    for(i = 0; i < N^2; i++)
    {
        unsigned int row = (int) i / 4, col = i % 4;
        dest[row][col] = src1[row][0] * src2[0][col] +
                         src1[row][1] * src2[1][col] +
                         ....
                         src[row][N-1] * src3[N-1][col];
    }
}

Donde N se reemplaza con el tamaño de la matriz. Entonces, si estás multiplicando matrices 4x4, entonces usas:

typedef float mat4[4][4];    

inline void mat4_mul(mat4 dest, mat4 src1, mat4 src2)
{
    unsigned int i;
    for(i = 0; i < 16; i++)
    {
        unsigned int row = (int) i / 4, col = i % 4;
        dest[row][col] = src1[row][0] * src2[0][col] +
                         src1[row][1] * src2[1][col] +
                         src1[row][2] * src2[2][col] +
                         src1[row][3] * src2[3][col];
    }
}

Esta función minimiza principalmente los bucles, pero el módulo podría estar sobrecargando ... En mi computadora, esta función se ejecutó aproximadamente un 50% más rápido que un triple para la función de multiplicación de bucles.

Contras:

• Se necesita mucho código (por ejemplo, diferentes funciones para mat3 x mat3, mat5 x mat5 ...)

• Ajustes necesarios para la multiplicación irregular (por ejemplo, mat3 x mat4) .....

Answer 7

No debes escribir la multiplicación de matrices. Debes depender de bibliotecas externas. En particular, debe usar la rutina GEMM de la biblioteca BLAS . GEMM a menudo proporciona las siguientes optimizaciones.

Bloqueo

La multiplicación de matrices eficiente se basa en el bloqueo de su matriz y en la realización de varios multiplicados de bloqueos más pequeños. Idealmente, el tamaño de cada bloque se elige para encajar bien en el caché, lo que mejora el rendimiento .

Sintonización

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

Optimización de nivel de montaje

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

Computación heterogénea (GPU)

Matrix Multiply es muy intensivo en FLOP / cómputo, por lo que es 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. Las personas dedican sus vidas a la mejora de estos algoritmos. Deberías usar su trabajo; Los hará felices.

Answer 8

Si está trabajando en números pequeños, entonces la mejora que está mencionando es insignificante. Además, el rendimiento variará dependiendo del hardware en el que se esté ejecutando. Pero si estás trabajando con números en millones, entonces tendrá efecto. Llegando al programa, puedes pegar el programa que has escrito.

Answer 9

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

Si la matriz es, digamos, 1,000x1,000, comenzarás a ver mejoras, pero diría que si está por debajo de 100x100 no debes preocuparte por eso.

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

Finalmente, si cambia la computadora que está utilizando por una más rápida, las diferencias serán aún más estrechas.

Answer 10

Solo algo para que lo intentes (pero esto solo haría una diferencia para matrices grandes): separe su lógica de suma de la lógica de multiplicación en el bucle interno, así:

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 canalización cuando escribes en suma. Por supuesto, gran parte de esto está a cargo del cambio de nombre de registro y similares, pero con mi conocimiento limitado del hardware, si quisiera eliminar cada onza de rendimiento del código, lo haría porque ahora no tiene que hacerlo. paralizar la tubería para esperar una escritura a suma. Dado que la multiplicación es más costosa que la suma, usted debe permitir que la máquina la paralice lo más posible, por lo que guardar sus paradas para la suma significa que pasa menos tiempo esperando en el ciclo de suma que en el ciclo de la multiplicación.

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

Answer 11

no tan especial pero mejor

    c = 0;
do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            sum = 0; sum_ = 0;
            for (l = 0; l < i; l++) {
                MatrixB[j][k] = MatrixB[k][j];
                sum += MatrixA[j][l]*MatrixB[k][l];
                l++;
                MatrixB[j][k] = MatrixB[k][j];
                sum_ += MatrixA[j][l]*MatrixB[k][l];

                sum += sum_;
            }
            MatrixR[j][k] = sum;
        }
     }
     c++;
} while (c<iteraciones);

Answer 12

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

Por supuesto, su código final es bastante velloso y usa muchas cosas específicas del sistema y ajustes de compilación, pero aún así, si realmente necesita velocidad, leer ese documento y leer su implementación definitivamente vale la pena.

Answer 13

ATENCIÓN: tienes un error 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];

sobrescribe MatrixB[0][1] y pierde ese valor! Cuando el bucle llega a f = 1, c = 0

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

El valor copiado es el mismo que ya estaba allí.