La moltiplicazione di matrice 3x3 di Laderman con solo 23 moltiplicazioni, ne vale la pena?

Question

Prendere il prodotto di due matrici 3x3 A*B=C. Ingenitamente ciò richiede 27 moltiplicazioni usando lo standard algorithm. Se uno era intelligente, puoi farlo usando solo 23 moltiplicazioni, a result found in 1973 by Laderman. La tecnica consiste nel salvare i passaggi intermedi e combinarli nel modo giusto.

Ora consente di correggere una lingua e un tipo, ad esempio C++ con elementi di double. Se l'algoritmo Laderman fosse hard-coded rispetto al semplice doppio ciclo, potremmo aspettarci che le prestazioni di un compilatore moderno facciano da contrappeso alle differenze degli algoritmi?

Note su questa domanda: Si tratta di un programmazione sito, e viene posta la domanda nel contesto delle migliori pratiche per un ciclo interno time-critical; ottimizzazione prematura non lo è. I suggerimenti sull'implementazione sono molto apprezzati come commenti.

Answer 1

Test Timing:

ho eseguito i tempi io stesso test ed i risultati mi ha sorpreso (da qui il motivo per cui ho posto la domanda in primo luogo). Il corto di esso è, sotto una compilazione standard il laderman è ~ 225% più veloce, ma con il flag di ottimizzazione -03 è 50% più lento! Ho dovuto aggiungere un elemento casuale alla matrice ogni volta durante il flag -O3 o il compilatore completamente ottimizzato via la semplice moltiplicazione, prendendo un tempo di zero entro la precisione dell'orologio. Dal momento che l'algoritmo laderman era un problema da controllare/ricontrollare, posterò il codice completo qui sotto per i posteri.

Specifiche: Ubuntu 12.04, Dell Prevision T1600, gcc.differenza percentuale in tempi:

• 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]

codice comparativa insieme attuazione 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; 
} 

Answer 2

La chiave è il controllo del set di istruzioni sulla piattaforma. Dipende dalla tua piattaforma. Esistono diverse tecniche e quando si tende ad aver bisogno delle massime prestazioni possibili, il compilatore verrà fornito con strumenti di profilazione, alcuni dei quali dotati di suggerimenti per l'ottimizzazione integrati. Per le operazioni con grana più fine, guarda l'output dell'assemblatore e controlla se ci sono miglioramenti anche a quel livello.

I comandi di più dati di istruzione simultanea eseguono la stessa operazione su più operandi in parallelo. In modo che si può prendere

double a,b,c,d; 
double w = d + a; 
double x = a + b; 
double y = b + c; 
double z = c + d; 

e sostituirlo con

double256 dabc = pack256(d, a, b, c); 
double256 abcd = pack256(a, b, c, d); 
double256 wxyz = dabc + abcd; 

in modo che quando i valori vengono caricati in registri, vengono caricati in un unico ampio registro a 256 bit per qualche piattaforma immaginario con 256 -bit registri ampi.

Il punto in virgola mobile è una considerazione importante, alcuni DSP possono essere operativi in ​​modo significativamente più veloce su numeri interi. Le GPU tendono ad essere grandiose in virgola mobile, anche se alcune sono 2x più veloci con una precisione unica. Il caso 3x3 di questo problema potrebbe rientrare in un singolo thread CUDA, quindi è possibile eseguire lo streaming di 256 di questi calcoli contemporaneamente.

Scegli la tua piattaforma, leggi la documentazione, attua diversi metodi e personalizzali.

Answer 3

Mi aspetto che il problema principale delle prestazioni sia la latenza della memoria. A double[9] è in genere 72 byte. Questa è già una quantità non banale e ne stai utilizzando tre.

Answer 4

Benché la questione di cui C++, ho implementato matrice 3x3 moltiplicazione C = A * B in C# (.NET 4.5) e ha eseguito alcuni test di temporizzazione di base sulla mia macchina Windows 7 a 64 bit con ottimizzazioni. 10.000.000 moltiplicazioni sono voluti circa

1. 0,556 secondi con un'implementazione ingenuo e
2. 0.874 secondi con il codice Laderman dall'altro risposta.

È interessante notare che il codice laderman era più lento del modo ingenuo. Non ho investigato con un profiler, ma immagino che le allocazioni extra siano più costose di alcune moltiplicazioni extra.

Sembra che i compilatori attuali siano abbastanza intelligenti da fare quelle ottimizzazioni per noi, il che è positivo. Ecco il codice ingenua che ho usato, per il suo interesse:

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

dove Matrix3D è una struct immutabile (in sola lettura doppie campi).

La cosa più difficile è quello di trovare una valida riferimento , dove si misura il codice e non, quello che il compilatore ha fatto con il codice (debugger con tonnellate di materiale extra, o ottimizzato senza il tuo codice vero e proprio dal momento che il risultato non è mai stato usato). Di solito cerco di "toccare" il risultato, in modo che il compilatore non possa rimuovere il codice sotto test (ad esempio, controllare gli elementi della matrice per l'uguaglianza con 89038.8989384 e lanciare, se uguale). Tuttavia, alla fine non sono nemmeno sicuro se il compilatore ha fatto fuori questo paragone :)