c++ - resueltos - La multiplicación de la matriz 3x3 de Laderman con solo 23 multiplicaciones, ¿vale la pena?
multiplicacion de matrices 3x3 ejercicios resueltos (4)
Tomar el producto de dos matrices 3x3 A*B=C Ingeniosamente esto requiere 27 multiplicaciones usando el algoritmo estándar . Si uno fuera inteligente, podría hacerlo usando solo 23 multiplicaciones, un resultado encontrado en 1973 por Laderman . La técnica consiste en guardar pasos intermedios y combinarlos de la manera correcta.
Ahora vamos a arreglar un idioma y un tipo, digamos C ++ con elementos de double . Si el algoritmo de Laderman fuera de código fijo en comparación con el doble bucle simple, ¿podríamos esperar que el rendimiento de un compilador moderno supere las diferencias de los algoritmos?
Notas sobre esta pregunta: Este es un sitio de programación , y la pregunta se formula en el contexto de las mejores prácticas para un bucle interno crítico en el tiempo; optimización prematura esto no es. Consejos para la implementación son muy bienvenidos como comentarios.
Pruebas de tiempo:
Yo mismo realicé las pruebas de sincronización y los resultados me sorprendieron (de ahí el motivo por el que hice la pregunta en primer lugar). En pocas palabras, bajo una compilación estándar, el laderman es ~ 225% más rápido, pero con el indicador de optimización -03 es 50% más lento . Tenía que agregar un elemento aleatorio a la matriz cada vez que durante el indicador -O3 o el compilador optimizaba por completo la simple multiplicación, tomando un tiempo de cero dentro de la precisión del reloj. Debido a que el algoritmo de laderman fue laderman de verificar / verificar, publicaré el código completo a continuación para la posteridad.
Especificaciones: Ubuntu 12.04, Dell Prevision T1600, gcc. Porcentaje de diferencia en los tiempos:
-
g++ [2.22, 2.23, 2.27] -
g++ -O3 [-0.48, -0.49, -0.48] -
g++ -funroll-loops -O3 [-0.48, -0.48, -0.47]
Código de referencia junto con la implementación de Laderman:
#include <iostream>
#include <ctime>
#include <cstdio>
#include <cstdlib>
using namespace std;
void simple_mul(const double a[3][3],
const double b[3][3],
double c[3][3]) {
int i,j,m,n;
for(i=0;i<3;i++) {
for(j=0;j<3;j++) {
c[i][j] = 0;
for(m=0;m<3;m++)
c[i][j] += a[i][m]*b[m][j];
}
}
}
void laderman_mul(const double a[3][3],
const double b[3][3],
double c[3][3]) {
double m[24]; // not off by one, just wanted to match the index from the paper
m[1 ]= (a[0][0]+a[0][1]+a[0][2]-a[1][0]-a[1][1]-a[2][1]-a[2][2])*b[1][1];
m[2 ]= (a[0][0]-a[1][0])*(-b[0][1]+b[1][1]);
m[3 ]= a[1][1]*(-b[0][0]+b[0][1]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][2]);
m[4 ]= (-a[0][0]+a[1][0]+a[1][1])*(b[0][0]-b[0][1]+b[1][1]);
m[5 ]= (a[1][0]+a[1][1])*(-b[0][0]+b[0][1]);
m[6 ]= a[0][0]*b[0][0];
m[7 ]= (-a[0][0]+a[2][0]+a[2][1])*(b[0][0]-b[0][2]+b[1][2]);
m[8 ]= (-a[0][0]+a[2][0])*(b[0][2]-b[1][2]);
m[9 ]= (a[2][0]+a[2][1])*(-b[0][0]+b[0][2]);
m[10]= (a[0][0]+a[0][1]+a[0][2]-a[1][1]-a[1][2]-a[2][0]-a[2][1])*b[1][2];
m[11]= a[2][1]*(-b[0][0]+b[0][2]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][1]);
m[12]= (-a[0][2]+a[2][1]+a[2][2])*(b[1][1]+b[2][0]-b[2][1]);
m[13]= (a[0][2]-a[2][2])*(b[1][1]-b[2][1]);
m[14]= a[0][2]*b[2][0];
m[15]= (a[2][1]+a[2][2])*(-b[2][0]+b[2][1]);
m[16]= (-a[0][2]+a[1][1]+a[1][2])*(b[1][2]+b[2][0]-b[2][2]);
m[17]= (a[0][2]-a[1][2])*(b[1][2]-b[2][2]);
m[18]= (a[1][1]+a[1][2])*(-b[2][0]+b[2][2]);
m[19]= a[0][1]*b[1][0];
m[20]= a[1][2]*b[2][1];
m[21]= a[1][0]*b[0][2];
m[22]= a[2][0]*b[0][1];
m[23]= a[2][2]*b[2][2];
c[0][0] = m[6]+m[14]+m[19];
c[0][1] = m[1]+m[4]+m[5]+m[6]+m[12]+m[14]+m[15];
c[0][2] = m[6]+m[7]+m[9]+m[10]+m[14]+m[16]+m[18];
c[1][0] = m[2]+m[3]+m[4]+m[6]+m[14]+m[16]+m[17];
c[1][1] = m[2]+m[4]+m[5]+m[6]+m[20];
c[1][2] = m[14]+m[16]+m[17]+m[18]+m[21];
c[2][0] = m[6]+m[7]+m[8]+m[11]+m[12]+m[13]+m[14];
c[2][1] = m[12]+m[13]+m[14]+m[15]+m[22];
c[2][2] = m[6]+m[7]+m[8]+m[9]+m[23];
}
int main() {
int N = 1000000000;
double A[3][3], C[3][3];
std::clock_t t0,t1;
timespec tm0, tm1;
A[0][0] = 3/5.; A[0][1] = 1/5.; A[0][2] = 2/5.;
A[1][0] = 3/7.; A[1][1] = 1/7.; A[1][2] = 3/7.;
A[2][0] = 1/3.; A[2][1] = 1/3.; A[2][2] = 1/3.;
t0 = std::clock();
for(int i=0;i<N;i++) {
// A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
simple_mul(A,A,C);
}
t1 = std::clock();
double tdiff_simple = (t1-t0)/1000.;
cout << C[0][0] << '' '' << C[0][1] << '' '' << C[0][2] << endl;
cout << C[1][0] << '' '' << C[1][1] << '' '' << C[1][2] << endl;
cout << C[2][0] << '' '' << C[2][1] << '' '' << C[2][2] << endl;
cout << tdiff_simple << endl;
cout << endl;
t0 = std::clock();
for(int i=0;i<N;i++) {
// A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3
laderman_mul(A,A,C);
}
t1 = std::clock();
double tdiff_laderman = (t1-t0)/1000.;
cout << C[0][0] << '' '' << C[0][1] << '' '' << C[0][2] << endl;
cout << C[1][0] << '' '' << C[1][1] << '' '' << C[1][2] << endl;
cout << C[2][0] << '' '' << C[2][1] << '' '' << C[2][2] << endl;
cout << tdiff_laderman << endl;
cout << endl;
double speedup = (tdiff_simple-tdiff_laderman)/tdiff_laderman;
cout << "Approximate speedup: " << speedup << endl;
return 0;
}
Aunque la pregunta mencionaba C ++, implementé la multiplicación de matriz 3x3 C = A * B en C # (.NET 4.5) y realicé algunas pruebas básicas de sincronización en mi máquina Windows 7 de 64 bits con optimizaciones. 10,000,000 multiplicaciones tomaron aproximadamente
- 0.556 segundos con una implementación ingenua y
- 0.874 segundos con el código laderman de la otra respuesta.
Curiosamente, el código de laderman era más lento que el ingenuo. No investigué con un generador de perfiles, pero creo que las asignaciones adicionales son más costosas que unas pocas multiplicaciones adicionales.
Parece que los compiladores actuales son lo suficientemente inteligentes como para hacer esas optimizaciones por nosotros, lo cual es bueno. Aquí está el código ingenuo que utilicé, para su interés:
public static Matrix3D operator *(Matrix3D a, Matrix3D b)
{
double c11 = a.M11 * b.M11 + a.M12 * b.M21 + a.M13 * b.M31;
double c12 = a.M11 * b.M12 + a.M12 * b.M22 + a.M13 * b.M32;
double c13 = a.M11 * b.M13 + a.M12 * b.M23 + a.M13 * b.M33;
double c21 = a.M21 * b.M11 + a.M22 * b.M21 + a.M23 * b.M31;
double c22 = a.M21 * b.M12 + a.M22 * b.M22 + a.M23 * b.M32;
double c23 = a.M21 * b.M13 + a.M22 * b.M23 + a.M23 * b.M33;
double c31 = a.M31 * b.M11 + a.M32 * b.M21 + a.M33 * b.M31;
double c32 = a.M31 * b.M12 + a.M32 * b.M22 + a.M33 * b.M32;
double c33 = a.M31 * b.M13 + a.M32 * b.M23 + a.M33 * b.M33;
return new Matrix3D(
c11, c12, c13,
c21, c22, c23,
c31, c32, c33);
}
donde Matrix3D es una estructura inmutable (campos dobles de solo lectura).
Lo complicado es encontrar un punto de referencia válido , donde mida su código y no, qué hizo el compilador con su código (depurador con toneladas de material adicional u optimizado sin su código real, ya que el resultado nunca se usó). Por lo general, trato de "tocar" el resultado, de manera que el compilador no puede eliminar el código que se está probando (por ejemplo, verifique que los elementos de la matriz sean iguales a 89038.8989384 y lancen, si son iguales). Sin embargo, al final ni siquiera estoy seguro de si el compilador piratea esta comparación a un lado :)
Espero que el principal problema de rendimiento sea la latencia de memoria. Un double[9] es típicamente de 72 bytes. Eso ya es una cantidad no trivial, y estás usando tres de ellos.
La clave es dominar el conjunto de instrucciones en su plataforma. Depende de tu plataforma. Existen varias técnicas, y cuando tiendes a necesitar el máximo rendimiento posible, tu compilador vendrá con herramientas de creación de perfiles, algunas de las cuales tienen sugerencias de optimización integradas. Para las operaciones más detalladas, mira la salida del ensamblador y observa si hay mejoras. en ese nivel también.
Los comandos de datos múltiples de instrucciones simultáneas realizan la misma operación en varios operandos en paralelo. Para que puedas tomar
double a,b,c,d;
double w = d + a;
double x = a + b;
double y = b + c;
double z = c + d;
y reemplazarlo con
double256 dabc = pack256(d, a, b, c);
double256 abcd = pack256(a, b, c, d);
double256 wxyz = dabc + abcd;
De modo que cuando los valores se cargan en registros, se cargan en un único registro de 256 bits de ancho para algunas plataformas ficticias con registros de 256 bits de ancho.
El punto flotante es una consideración importante, algunos DSP pueden ser significativamente más rápidos operando en enteros. Las GPU tienden a ser excelentes en punto flotante, aunque algunas son 2 veces más rápidas con una sola precisión. El caso 3x3 de este problema podría caber en un solo hilo CUDA, por lo que podría transmitir 256 de estos cálculos simultáneamente.
Elija su plataforma, lea la documentación, implemente varios métodos diferentes y perfílelos.