Deze keer een lekker low level onderwerp, performance. Voor de efficiëntie van een programma is het belangrijk om lineair door het geheugen te gaan. Dit is de basis regel. Helaas is dat niet altijd mogelijk, bijvoorbeeld bij een matrix matrix vermenigvuldiging. Hiervoor zijn er allerlei technieken om de snelste implementatie te vinden. Zie bijvoorbeeld wikipedia Loop_nest_optimization. Er zijn vele soortgelijke pagina's te vinden. De voorbeeld pagina geeft een redelijk complexe uitleg waarom bepaalde technieken werken, zonder meetdata te geven. Nu wil ik precies het tegenovergestelde doen. Namelijk meten welke vorm van matrix matrix vermenigvuldiging het snelst is, zonder een verklaring te geven.

Eerst de meest eenvoudige implementatie:

void simple(const size_t size, const std::vector<double>& a, const std::vector<double>& b, std::vector<double>& c)
{
  for (size_t i = 0; i < size; ++i) {
    for (size_t j = 0; j < size; ++j) {
       c[i*size + j] = 0.;
       for (size_t k = 0; k < size; ++k) {
         c[i*size + j] += a[i*size + k] * b[k*size + j];
       }
    }
  }
}

In dit voorbeeld zit het probleem bij matrix 'b'. Deze matrix word niet in de juiste volgorde ingelezen. Deze implementatie geeft de volgende resultaten als het gecompileerd word met "g++ -O0".

Het aantal uitgevoerde berekeningen schaalt met NxNxN, vandaar dat de tijd hiermee gedeeld word. Wat hier opvalt zijn de pieken bij N=256 en N=512. Meestal zijn machten van twee gunstige getallen voor berekeningen, maar hier is dat duidelijk niet het geval. Natuurlijk is het niet echt verstandig om naar de performance te kijken van een programma dat met -O0 gecompileerd is. Dus nu de resultaten wanneer er gecompileerd word met "g++ -O3"

De berekening is ruim 4 keer zo snel, zodat in ieder geval bevestigd kan worden dat -O3 inderdaad sneller is als -O0. De pieken bij N=256 en N=512 zijn nog steeds zichtbaar. Wat echter opvalt is dat er een enorme ruis zichtbaar is. Als we nauwkeuriger kijken dan blijkt dit helemaal geen ruis te zijn. Voor even N blijkt de berekening significant sneller te zijn dan voor oneven N, bijvoorbeeld:

Waarbij de berekening 100 keer is gedaan. In dit geval lijkt het effect veroorzaakt door de compiler. Een ander effect dat zichtbaar is in beide grafieken is dat de tijd onevenredig veel toeneemt tussen N=250 en N=300. Dit kan verder uitgewerkt worden door het experiment te herhalen met floats in plaats van doubles:

De berekening is, zoals verwacht, sneller voor floats. De grafiek voor floats laat hetzelfde effect zien, maar ditmaal tussen N=350 en N=400. Dit is ongeveer √2 later. Na deze inleidende beschietingen keren we terug naar het oorspronkelijke voorbeeld artiekel: wikipedia Loop_nest_optimization

Laten we een van de besproken technieken toepassen:

void opt2x2(const size_t size, const std::vector<double>& a, const std::vector<double>& b, std::vector<double>& c)
{
  for (size_t i = 0; i < size-1; i+=2) {
    for (size_t j = 0; j < size-1; j+=2) {
       double acc00 = 0.;
       double acc01 = 0.;
       double acc10 = 0.;
       double acc11 = 0.;
       for (size_t k = 0; k < size; ++k) {
         acc00 += a[i*size + k] * b[k*size + j];
         acc10 += a[(i+1)*size + k] * b[k*size + j];
         acc01 += a[i*size + k] * b[k*size + j+1];
         acc11 += a[(i+1)*size + k] * b[k*size + j+1];
       }
       c[i*size + j] += acc00;
       c[i*size + j+1] += acc01;
       c[(i+1)*size + j] += acc10;
       c[(i+1)*size + j+1] += acc11;
    }
  }
  if (size%2)
  {
     size_t i = size-1;
     for (size_t j = 0; j < size; ++j) {
       c[i*size + j] = 0.;
       for (size_t k = 0; k < size; ++k) {
         c[i*size + j] += a[i*size + k] * b[k*size + j];
       }
     }
     size_t j = size-1;
     for (size_t i = 0; i < size; ++i) {
       c[i*size + j] = 0.;
       for (size_t k = 0; k < size; ++k) {
         c[i*size + j] += a[i*size + k] * b[k*size + j];
       }
     }
  }
}

De truc is om matrix 'c' per 2x2 sub-matrix uit te rekenen. Dit trucje kan ook voor grotere sub-matrices toegepast worden. De code zal ik hier niet laten zien, alleen de resultaten.

Dit is weer een bevestiging van meten is weten. Het wikipedia artikel suggereert dat de snelheid afneemt wanneer de sub-matrices te groot worden gemaakt. Op mijn machine lijkt dit niet het geval te zijn. Hoe meer de loop uitgeschreven word hoe sneller de berekening is. Bovendien lijkt het dat het bizarre effect bij N=256 en N=512 kleiner word als de loop deels uitgeschreven word. Dit valt nog meer op als er met g++ -O0 gecompileerd word:

Helaas geen duidelijke verklaring, maar misschien wel een inspiratie om performance te meten.